Prueba de software estadístico

¿Qué técnicas / enfoques son útiles para probar el software estadístico? Estoy particularmente interesado en los programas que realizan estimaciones paramétricas utilizando la máxima probabilidad.

No siempre es posible comparar resultados con los de otros programas o fuentes publicadas, ya que la mayoría de las veces cuando escribo un programa propio es porque el cálculo que necesito no está implementado en un sistema existente.

No estoy insistiendo en métodos que puedan garantizar la corrección. Estaría contento con las técnicas que pueden detectar una fracción de los errores.

Una técnica útil es la prueba de Monte Carlo. Si hay dos algoritmos que hacen lo mismo, implemente ambos, alimente con datos aleatorios y verifique que (dentro de una pequeña tolerancia para el fuzz numérico) produzcan la misma respuesta. He hecho esto varias veces antes:

• Escribí una implementación de eficiente pero difícil de implementar de Kendall's Tau B. Para probarlo, escribí una implementación de 50 líneas muy simple que se ejecutó en . $O(N\ log\ N)$$O(N^2)$

• Escribí un código para hacer la regresión de cresta. El mejor algoritmo para hacer esto depende de si estás en el caso o , por lo que necesitaba dos algoritmos de todos modos. $n > p$$p > n$

En ambos casos, estaba implementando técnicas relativamente conocidas en el lenguaje de programación D (para el cual no existía ninguna implementación), así que también verifiqué algunos resultados contra R. Sin embargo, la prueba de Monte Carlo detectó errores que nunca habría detectado de otra manera. .

Otra buena prueba es afirma . Es posible que no sepa exactamente cuáles deberían ser los resultados correctos de su cálculo, pero eso no significa que no pueda realizar comprobaciones de cordura en varias etapas del cálculo. En la práctica, si tiene muchos de estos en su código y todos pasan, entonces el código generalmente es correcto.

Editar: Un tercer método es alimentar los datos del algoritmo (sintéticos o reales) donde usted sabe al menos aproximadamente cuál es la respuesta correcta, incluso si no sabe exactamente, y ver mediante inspección si la respuesta es razonable. Por ejemplo, es posible que no sepa exactamente cuáles son las estimaciones de sus parámetros, pero puede saber cuáles se supone que son "grandes" y cuáles se supone que son "pequeños".

No estoy seguro de si esta es realmente una respuesta a su pregunta, pero al menos está relacionada tangencialmente.

Mantengo el paquete de estadísticas en Maple . Un ejemplo interesante de código difícil de probar es la generación aleatoria de muestras de acuerdo con diferentes distribuciones; es fácil comprobar que no se generan errores, pero es más difícil determinar si las muestras que se generan cumplen con la distribución solicitada "suficientemente bien". Dado que Maple tiene características simbólicas y numéricas, puede usar algunas de las características simbólicas para probar la generación de muestras (puramente numéricas):

1. Hemos implementado algunos tipos de pruebas de hipótesis estadísticas, una de las cuales es la prueba del modelo adecuado de chi cuadrado : una prueba de chi cuadrado de los números de muestras en contenedores determinados a partir del CDF inverso de la distribución de probabilidad dada. Entonces, por ejemplo, para probar la generación de muestras de distribución de Cauchy, ejecuto algo como

   with(Statistics):
   infolevel[Statistics] := 1:
   distribution := CauchyDistribution(2, 3):
   sample := Sample(distribution, 10^6):
   ChiSquareSuitableModelTest(sample, distribution, 'bins' = 100, 'level' = 0.001);
   

   Debido a que puedo generar una muestra tan grande como quiera, puedo hacer que bastante pequeño.$\alpha$

2. Para distribuciones con momentos finitos, calculo por un lado una cantidad de momentos de muestra y, por otro lado, calculo simbólicamente los momentos de distribución correspondientes y su error estándar. Entonces, por ejemplo, la distribución beta:

   with(Statistics):
   distribution := BetaDistribution(2, 3):
   distributionMoments := Moment~(distribution, [seq(1 .. 10)]);
   standardErrors := StandardError[10^6]~(Moment, distribution, [seq(1..10)]);
   evalf(distributionMoments /~ standardErrors);
   

   Esto muestra una lista decreciente de números, el último de los cuales es 255.1085766. Entonces, incluso para el décimo momento, el valor del momento es más de 250 veces el valor del error estándar del momento muestral para una muestra de tamaño . Esto significa que puedo implementar una prueba que se ejecuta más o menos de la siguiente manera:$10^6$

   with(Statistics):
   sample := Sample(BetaDistribution(2, 3), 10^6):
   sampleMoments := map2(Moment, sample, [seq(1 .. 10)]);
   distributionMoments := [2/5, 1/5, 4/35, 1/14, 1/21, 1/30, 4/165, 1/55, 2/143, 1/91];
   standardErrors := 
     [1/5000, 1/70000*154^(1/2), 1/210000*894^(1/2), 1/770000*7755^(1/2), 
      1/54600*26^(1/2), 1/210000*266^(1/2), 7/5610000*2771^(1/2), 
      1/1567500*7809^(1/2), 3/5005000*6685^(1/2), 1/9209200*157366^(1/2)];
   deviations := abs~(sampleMoments - distributionMoments) /~ standardErrors;
   

   Los números entran distributionMomentsy standardErrorsprovienen de la primera ejecución anterior. Ahora, si la generación de la muestra es correcta, los números en las desviaciones deberían ser relativamente pequeños. Supongo que están aproximadamente distribuidos normalmente (que en realidad no lo son, pero se acerca lo suficiente; recuerde que estas son versiones escaladas de momentos de muestra, no las muestras en sí) y, por lo tanto, puedo marcar un caso en el que hay una desviación mayor que 4: correspondiente a un momento de muestra que se desvía más de cuatro veces el error estándar del momento de distribución. Es muy poco probable que esto ocurra al azar si la generación de la muestra es buena. Por otro lado, si los primeros 10 momentos de muestra coinciden con los momentos de distribución en menos de medio por ciento, tenemos una aproximación bastante buena de la distribución.

La clave de por qué ambos métodos funcionan es que el código de generación de muestra y el código simbólico son casi completamente disjuntos. Si hubiera una superposición entre los dos, entonces un error en esa superposición podría manifestarse tanto en la generación de la muestra como en su verificación, y por lo tanto no quedar atrapado.

Bruce McCullough tenía un poco de una industria casera en la evaluación de software estadístico (en el sentido más amplio; también probó Microsoft Excel. Y lo encontró deficiente). Dos documentos que ilustran parte de su enfoque están aquí y aquí.

El presidente de StataCorp, William Gould, da muchos detalles en este artículo de Stata Journal. 1 Es un artículo muy interesante sobre control de calidad de software estadístico.