Test de logiciels statistiques

Quelles techniques / approches sont utiles pour tester des logiciels statistiques? Je m'intéresse particulièrement aux programmes qui font une estimation paramétrique en utilisant le maximum de vraisemblance.

Comparer les résultats à ceux d'autres programmes ou de sources publiées n'est pas toujours possible car la plupart du temps lorsque j'écris un programme personnel, c'est parce que le calcul dont j'ai besoin n'est pas déjà implémenté dans un système existant.

Je n'insiste pas sur les méthodes qui peuvent garantir l'exactitude. Je serais ravi des techniques qui peuvent détecter une fraction des erreurs.

Une technique utile est le test de monte carlo. S'il y a deux algorithmes qui font la même chose, implémentez-les tous les deux, alimentez-les en données aléatoires et vérifiez que (à une petite tolérance près pour les fuzz numériques) ils produisent la même réponse. Je l'ai fait plusieurs fois auparavant:

• J'ai écrit une implémentation efficace mais difficile à implémenter du Tau B de Kendall. Pour le tester, j'ai écrit une implémentation de 50 lignes extrêmement simple qui s'exécutait en . O ( N 2$O(N\ log\ N)$$O(N^2)$

• J'ai écrit du code pour faire une régression de crête. Le meilleur algorithme pour ce faire dépend de si vous êtes dans le cas ou , donc j'avais besoin de deux algorithmes de toute façon. p > $n > p$$p > n$

Dans ces deux cas, j'implémentais des techniques relativement connues dans le langage de programmation D (pour lesquelles aucune implémentation n'existait), j'ai donc également vérifié quelques résultats par rapport à R. Néanmoins, le test de monte carlo a détecté des bogues que je n'aurais jamais détectés autrement .

Un autre bon test est l' affirmation . Vous ne savez peut-être pas exactement quels devraient être les résultats corrects de votre calcul, mais cela ne signifie pas que vous ne pouvez pas effectuer de contrôles d'intégrité à différentes étapes du calcul. En pratique, si vous en avez beaucoup dans votre code et qu'ils réussissent tous, le code est généralement correct.

Edit: Une troisième méthode consiste à alimenter les données de l'algorithme (synthétique ou réel) où vous savez au moins approximativement quelle est la bonne réponse, même si vous ne savez pas exactement, et voyez par inspection si la réponse est raisonnable. Par exemple, vous ne savez peut-être pas exactement quelles sont les estimations de vos paramètres, mais vous savez peut-être lesquelles sont censées être "grandes" et lesquelles sont supposées être "petites".

Je ne sais pas si c'est vraiment une réponse à votre question, mais c'est au moins tangentiellement lié.

Je maintiens le package Statistics dans Maple . Un exemple intéressant de code difficile à tester est la génération d'échantillons aléatoires selon différentes distributions; il est facile de tester qu'aucune erreur n'est générée, mais il est plus délicat de déterminer si les échantillons générés sont "assez bien" conformes à la distribution demandée. Étant donné que Maple possède des caractéristiques symboliques et numériques, vous pouvez utiliser certaines des caractéristiques symboliques pour tester la génération d'échantillon (purement numérique):

1. Nous avons mis en œuvre quelques types de tests d'hypothèses statistiques, dont l'un est le test de modèle adapté au chi carré - un test du chi carré du nombre d'échantillons dans des casiers déterminé à partir du CDF inverse de la distribution de probabilité donnée. Ainsi, par exemple, pour tester la génération d'échantillons de distribution de Cauchy, je lance quelque chose comme

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

   Parce que je peux générer un échantillon aussi grand que je le souhaite, je peux rendre assez petit.$\alpha$

2. Pour les distributions à moments finis, je calcule d'une part un certain nombre d'exemples de moments, et d'autre part, je calcule symboliquement les moments de distribution correspondants et leur erreur standard. Ainsi, par exemple pour la distribution bêta:

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

   Cela montre une liste décroissante de nombres, dont le dernier est 255.1085766. Ainsi, même pour le 10e moment, la valeur du moment est plus de 250 fois la valeur de l'erreur standard du moment d'échantillon pour un échantillon de taille . Cela signifie que je peux implémenter un test qui s'exécute plus ou moins comme suit:$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;
   

   Les chiffres distributionMomentset standardErrorsproviennent de la première manche ci-dessus. Maintenant, si la génération d'échantillon est correcte, les nombres d'écarts devraient être relativement petits. Je suppose qu'ils sont à peu près normalement distribués (ce qui n'est pas vraiment le cas, mais cela se rapproche suffisamment - rappelez-vous que ce sont des versions à l'échelle des exemples de moments, pas les échantillons eux-mêmes) et ainsi je peux, par exemple, signaler un cas où une déviation est supérieur à 4 - correspondant à un moment d'échantillon qui s'écarte de plus de quatre fois l'erreur standard du moment de distribution. Il est très peu probable que cela se produise au hasard si la génération d'échantillon est bonne. D'un autre côté, si les 10 premiers moments de l'échantillon correspondent aux moments de distribution à moins d'un demi pour cent, nous avons une assez bonne approximation de la distribution.

La raison pour laquelle ces deux méthodes fonctionnent est que l'exemple de code de génération et le code symbolique sont presque complètement disjoints. S'il y avait un chevauchement entre les deux, une erreur dans ce chevauchement pourrait se manifester à la fois dans la génération de l'échantillon et dans sa vérification, et donc ne pas être détectée.

Bruce McCullough avait un peu d'industrie artisanale dans l'évaluation des logiciels statistiques (au sens le plus large; il a également testé Microsoft Excel. Et il l'a trouvé insuffisant). Deux articles qui illustrent une partie de son approche sont ici et ici.

Beaucoup de détails sont donnés par le président de StataCorp, William Gould, dans cet article du Stata Journal. 1 Il s'agit d'un article très intéressant sur le contrôle qualité des logiciels statistiques.