Comment calculer la probabilité associée à des scores Z absurdement élevés?

14

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 1erf(n/2)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 ).

Douglas S. Stones
la source
6
Ce n'est peut-être pas la bonne question à poser. Ces z-scores sont faux car ils supposent que la distribution normale est une bien meilleure approximation, ou modèle, qu'elle ne l'est en réalité. C'est un peu comme si la mécanique newtonienne était bonne à 600 000 décimales. Si vous êtes en effet uniquement intéressé par le calcul de erf pour les valeurs extrêmes de n , alors cette question appartient à math.SE, pas ici.
whuber
6
Pour des valeurs "absurdement" grandes, vous ne ferez pas mieux que d'utiliser la borne supérieure Pr(Z>z)(z2π)1ez2/2 pour virgule flottantedouble précision. Cette approximation et d'autres sont discutées ailleurs sur stats.SE.
Cardinal
Merci cardinal, cette limite semble être tout à fait exacte. Pourquoi n'en faites-vous pas une réponse?
Douglas S. Stones
@Douglas: Si vous êtes toujours intéressé, je peux mettre quelque chose en place le lendemain ou alors et le poster comme une réponse plus complète.
cardinal
1
Eh bien ... je pense que ça vaudrait la peine de l'ajouter comme réponse. Peut-être que la limite est connue dans les statistiques prob +, mais je ne le savais pas. De plus, les Q et A ici ne sont pas uniquement destinés à l'OP.
Douglas S. Stones

Réponses:

19

La question concerne la fonction d'erreur complémentaire

erfc(x)=2πxexp(t2)dt

pour les "grandes" valeurs de x ( =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,

f(x)=11exp(x2(4+ax2π+ax2)),

a=8(π3)3(4π)0.439862.

Bien que ce soit une excellente approximation de la fonction d'erreur elle-même, c'est une terrible approximation de erfc . Cependant, il existe un moyen de corriger systématiquement cela.

Pour les valeurs de p associées à de si grandes valeurs de x , 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 0x5.8 :

Plot 1

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.xexp(5.82)1014.6x

L'exécution du calcul en arithmétique étendue (avec Mathematica ) améliore notre image de ce qui se passe:

Plot 2

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!xx=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 )xerfcfx xxerfc(x)/f(x), car l'espoir est que cela devrait avoir une valeur limite constante. Voici son graphique:

Plot 3

Notre supposition semble être confirmée: ce ratio semble approcher d'une limite autour de 8 environ. Lorsque demandé, Mathematica le fournira:

a1 = Limit[x (Erfc[x]/f[x]), x -> \[Infinity]]

La valeur est . Cela nous permet d'améliorer l'estimation:nous prenonsa1=2πe3(4+π)28(3+π)7.94325

f1(x)=f(x)a1x

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 2x5.320001erfc(x)/f1(x)1/x2xx2 et trouvez la prochaine limite:

a2 = Limit[x^2 (a1 - x (Erfc[x]/f[x])), x -> \[Infinity]] 

La valeur est

a2=132πe3(4+π)28(3+π)(329(4+π)3π(3+π)2)114.687.

Ce processus peut se poursuivre aussi longtemps que nous le souhaitons. Je l'ai sorti encore une fois, trouvant

a3 = Limit[x^2 (a2 - x^2 (a1 - x (Erfc[x]/f[x]))), x -> \[Infinity]] 

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

f3(x)=f(x)(a1a2/x2+a3/x4)/x.

L'erreur est proportionnelle à . D'importation est la constante de proportionnalité, nous traçons donc x 6 ( 1 - erfc ( x ) / f 3 ( x ) ) :x6x6(1erfc(x)/f3(x))

Plot 4

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 :f3erfc(x)2661/x6x>0xxx1020

 x  Erfc    Approximation      
10  2.088*10^-45    2.094*10^-45
11  1.441*10^-54    1.443*10^-54
12  1.356*10^-64    1.357*10^-64
13  1.740*10^-75    1.741*10^-75
14  3.037*10^-87    3.038*10^-87
15  7.213*10^-100   7.215*10^-100
16  2.328*10^-113   2.329*10^-113
17  1.021*10^-127   1.021*10^-127
18  6.082*10^-143   6.083*10^-143
19  4.918*10^-159   4.918*10^-159
20  5.396*10^-176   5.396*10^-176

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=8NormSDist

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,fx

f(x)12exp(x2(4+ax2π+ax2)).

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

log10(f(x))(10002(4+a10002π+a10002)log(2))/log(10)434295.63047.

Exposer les rendements

f(1000)2.3416910434296.

L'application de la correction (en ) produitf3

erfc(1000)1.86003 70486 3232810434298.

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/x1%, égal à1,86003810 - 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.86003810434298

whuber
la source
1
+1 Ceci est une excellente réponse, d'une manière ou d'une autre, je n'ai jamais rencontré ce fil avant.
amibe dit Réintégrer Monica
15

Une borne supérieure simple

z>0

S(z):=P(Z>z)=zφ(z)dz,
where φ(z)=(2π)1/2ez2/2 is the standard normal pdf. I've used the notation S(z) in deference to the standard notation in survival analysis. In engineering contexts, they call this function the Q-function and denote it by Q(z).

Then, a very simple, elementary upper bound is

S(z)φ(z)z=:S^u(z),
where the notation on the right-hand side indicates this is an upper-bound estimate. This answer gives a proof of the bound.

There are several nice complementary lower bounds as well. One of the handiest and easiest to derive is the bound

S(z)zz2+1φ(z)=:S^(z).
There are at least three separate methods for deriving this bound. A rough sketch of one such method can be found in this answer to a related question.

A picture

Below is a plot of the two bounds (in grey) along with the actual function S(z).

Upper-tail of normal and bounds

How good is it?

From the plot, it seems that the bounds become quite tight even for moderately large z. 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

E(z)=|S^u(z)S(z)S(z)|.
This gives you the proportional error of the estimate.

Now, note that, since all of the involved functions are nonnegative, by using the bounding properties of S^u(z) and S^(z), we get

E(z)=S^u(z)S(z)S(z)S^u(z)S^(z)S^(z)=z2,
and so this provides a proof that for z10 the upper-bound is correct to within 1%, for z28 it is correct to within 0.1% and for z100 it is correct to within 0.01%.

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 on S(z) of the form R(z)φ(z) where R(z) is a rational function.

Finally, here is another somewhat-related question and answer.

cardinal
la source
1
Apologies for all the "self-citations". Once, several years ago, I took an intense, two-week-long interest in related questions and tried to learn as much as I could about this topic.
cardinal
+1 Agree with whuber. Very nice, and I appreciate the links to other answers.
Iterator
5

You can approximate it with much simpler functions - see this Wikipedia section for more information. The basic approximation is that erf(x)sgn(x)1exp(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.

Iterator
la source
1
Some amplification of this would be welcome, for two reasons. First, it's best when answers can stand alone. Second, that article writes ambiguously about the quality of the approximation "in a neighborhood of infinity": just how accurate is "very accurate"? (You implicitly have a good sense of this, but it's a lot to expect of all interested readers.) The stated value of ".00035" is useless here.
whuber
Thanks. I didn't notice that there was Javascript-based support for using TeX, which made the difference in writing that out.
Iterator
1
Incidentally, the Wikipedia reference to that approximation is broken. Mathematica finds, though, that the relative error (1 - approx(x)/erf(x)) behaves like the reciprocal of 2exp(x2+3(π4)2/(8(π3))).
whuber
@whuber, can you post the Mathematica code for that? :) I haven't seen Mathematica in 15+ years, and never for this kind of purpose.
Iterator
I posted it in a separate reply.
whuber