Les progiciels pour la détection de motifs de réseau peuvent renvoyer des scores Z extrêmement élevés (le plus élevé que j'ai vu est 600 000+, mais des scores Z de plus de 100 sont assez courants). Je prévois de montrer que ces scores Z sont faux.
D'énormes scores Z correspondent à des probabilités associées extrêmement faibles. Les valeurs des probabilités associées sont données par exemple sur la page wikipedia de distribution normale (et probablement tous les manuels de statistiques) pour des scores Z allant jusqu'à 6. Donc ...
Question : Comment calcule-t-on la fonction d'erreur pour n jusqu'à 1 000 000, disons?
Je suis particulièrement à la recherche d'un package déjà implémenté pour cela (si possible). Le meilleur que j'ai trouvé jusqu'à présent est WolframAlpha, qui parvient à le calculer pour n = 150 ( ici ).
la source
Réponses:
La question concerne la fonction d'erreur complémentaire
pour les "grandes" valeurs dex ( =n/2–√ dans la question initiale) - c'est-à-dire entre 100 et 700 000 environ. (Dans la pratique, toute valeur supérieure à environ 6 doit être considérée comme "grande", comme nous le verrons.) Notez que parce que cela sera utilisé pour calculer les valeurs de p, il y a peu de valeur à obtenir plus de trois chiffres significatifs (décimaux) .
Pour commencer, considérez l'approximation suggérée par @Iterator,
où
Bien que ce soit une excellente approximation de la fonction d'erreur elle-même, c'est une terrible approximation deerfc . Cependant, il existe un moyen de corriger systématiquement cela.
Pour les valeurs de p associées à de si grandes valeurs dex , nous nous intéressons à l' erreur relative f(x)/erfc(x)−1 : nous espérons que sa valeur absolue serait inférieure à 0,001 pour trois chiffres significatifs de précision. Malheureusement, cette expression est difficile à étudier pour les grands x raison de sous-écoulements dans le calcul en double précision. Voici une tentative, qui trace l'erreur relative par rapport à x pour 0≤x≤5.8 :
Le calcul devient instable une fois que dépasse 5,3 environ et ne peut pas délivrer un chiffre significatif au-delà de 5,8. Ce n'est pas une surprise: exp ( - 5,8 2 ) ≈ 10 - 14,6 repousse les limites de l'arithmétique double précision. Parce qu'il n'y a aucune preuve que l'erreur relative va être suffisamment petite pour un x plus grand , nous devons faire mieux.x exp(−5.82)≈10−14.6 x
L'exécution du calcul en arithmétique étendue (avec Mathematica ) améliore notre image de ce qui se passe:
L'erreur augmente rapidement avec et ne montre aucun signe de stabilisation. Au-delà de x = 10 environ, cette approximation ne fournit même pas un chiffre fiable d'informations!x x=10
Cependant, l'intrigue commence à sembler linéaire. On pourrait deviner que l'erreur relative est directement proportionnelle à . (Cela a du sens sur le plan théorique: erfc est manifestement une fonction impaire et f est manifestement pair, donc leur rapport devrait être une fonction impaire. Ainsi, nous nous attendrions à ce que l'erreur relative, si elle augmente, se comporte comme une puissance impaire de x .) Cela nous amène à étudier l'erreur relative divisée par x . De manière équivalente, je choisis d'examiner x ⋅ erfc ( x ) / f ( x )x erfc f x x x⋅erfc(x)/f(x) , car l'espoir est que cela devrait avoir une valeur limite constante. Voici son graphique:
Notre supposition semble être confirmée: ce ratio semble approcher d'une limite autour de 8 environ. Lorsque demandé, Mathematica le fournira:
La valeur est . Cela nous permet d'améliorer l'estimation:nous prenonsa1=2π√e3(−4+π)28(−3+π)≈7.94325
comme le premier raffinement de l'approximation. Lorsque est vraiment grand - supérieur à quelques milliers - cette approximation est très bien. Parce que cela ne sera toujours pas assez bon pour une gamme intéressante d'arguments entre 5.3 et 2000 , itérons la procédure. Cette fois, l'erreur relative inverse - en particulier, l'expression 1 - erfc ( x ) / f 1 ( x ) - devrait se comporter comme 1 / x 2 pour les grands x (en vertu des considérations de parité précédentes). En conséquence, nous multiplions par x 2x 5.3 2000 1−erfc(x)/f1(x) 1/x2 x x2 et trouvez la prochaine limite:
La valeur est
Ce processus peut se poursuivre aussi longtemps que nous le souhaitons. Je l'ai sorti encore une fois, trouvant
avec une valeur d'environ 1623,67. (L'expression complète implique une fonction rationnelle de degré huit de et est trop longue pour être utile ici.)π
Le déroulement de ces opérations donne notre approximation finale
L'erreur est proportionnelle à . D'importation est la constante de proportionnalité, nous traçons donc x 6 ( 1 - erfc ( x ) / f 3 ( x ) ) :x−6 x6(1−erfc(x)/f3(x))
Il s'approche rapidement d'une valeur limite autour de 2660,59. En utilisant l'approximation , nous obtenons des estimations de erfc ( x ) dont la précision relative est meilleure que 2661 / x 6 pour tout x > 0 . Une fois que x dépasse 20 ou plus, nous avons nos trois chiffres significatifs (ou bien plus, car x devient plus grand). À titre de vérification, voici un tableau comparant les valeurs correctes à l'approximation de x entre 10 et 20 :f3 erfc(x) 2661/x6 x>0 x x x 10 20
En fait, cette approximation fournit au moins deux chiffres de précision significatifs pour , ce qui est à peu près là où les calculs des piétons (tels que la fonction d'Excel ) s'épuisent.x=8
NormSDist
Enfin, on pourrait s'inquiéter de notre capacité à calculer l'approximation initiale . Cependant, ce n'est pas difficile: lorsque x est suffisamment grand pour provoquer des débordements dans l'exponentielle, la racine carrée est bien approximée de la moitié de l'exponentielle,f x
Le calcul du logarithme de celui-ci (en base 10) est simple et donne facilement le résultat souhaité. Par exemple, soit . Le logarithme commun de cette approximation estx=1000
Exposer les rendements
L'application de la correction (en ) produitf3
Notez que la correction réduit l'approximation d'origine de plus de 99% (et en effet, .) (Cette approximation ne diffère de la valeur correcte que dans le dernier chiffre. Une autre approximation bien connue, exp ( - x 2 ) / ( x √a1/x≈1% , égal à1,860038⋅10 - 434298 , errant dans le sixième chiffre significatif. Je suis sûr que nous pourrions aussi améliorer celui-ci, si nous le voulions, en utilisant les mêmes techniques.)exp(−x2)/(xπ−−√) 1.860038⋅10−434298
la source
Une borne supérieure simple
Then, a very simple, elementary upper bound is
There are several nice complementary lower bounds as well. One of the handiest and easiest to derive is the bound
A picture
Below is a plot of the two bounds (in grey) along with the actual functionS(z) .
How good is it?
From the plot, it seems that the bounds become quite tight even for moderately largez . We might ask ourselves how tight they are and
what sort of quantitative statement in that regard can be made.
One useful measure of tightness is the absolute relative error
Now, note that, since all of the involved functions are nonnegative, by using the bounding properties ofS^u(z) and S^ℓ(z) , we get
In fact, the simple form of the bounds provides a good check on other "approximations". If, in the numerical calculation of more complicated approximations, we get a value outside these bounds, we can simply "correct" it to take the value of, e.g., the upper bound provided here.
There are many refinements of these bounds. The Laplace bounds mentioned here provide a nice sequence of upper and lower bounds onS(z) of the form R(z)φ(z) where R(z) is a rational function.
Finally, here is another somewhat-related question and answer.
la source
You can approximate it with much simpler functions - see this Wikipedia section for more information. The basic approximation is thaterf(x)≈sgn(x)1−exp(−x24/π+ax21+ax2−−−−−−−−−−−−−−−−√)
The article has an incorrect link for that section. The PDF referenced can be found in Sergei Winitzki's files - or at this link.
la source