Wie berechnet man die Wahrscheinlichkeit, die mit absurd großen Z-Scores verbunden ist?

14

Softwarepakete zur Erkennung von Netzwerkmotiven können enorm hohe Z-Scores liefern (der höchste Wert, den ich gesehen habe, ist 600.000+, aber Z-Scores von mehr als 100 sind durchaus üblich). Ich habe vor zu zeigen, dass diese Z-Scores falsch sind.

Riesige Z-Scores entsprechen extrem niedrigen zugehörigen Wahrscheinlichkeiten. Die Werte der zugehörigen Wahrscheinlichkeiten sind zB auf der Normalverteilungs-Wikipedia-Seite (und wahrscheinlich in jedem Statistik-Lehrbuch) für Z-Scores von bis zu 6 angegeben. Also ...

Frage : Wie berechnet man die Fehlerfunktion 1erf(n/2)Sagen wir für n bis zu 1.000.000?

Ich bin besonders nach einem bereits implementierten Paket dafür (wenn möglich). Das beste, das ich bisher gefunden habe, ist WolframAlpha, das es schafft, es für n = 150 ( hier ) zu berechnen .

Douglas S. Stones
quelle
6
Vielleicht ist dies nicht die richtige Frage. Diese Z-Scores sind falsch, weil sie davon ausgehen, dass die Normalverteilung eine weitaus bessere Annäherung oder ein besseres Modell darstellt als sie tatsächlich ist. Es ist ein bisschen so, als würde man annehmen, dass die Newtonsche Mechanik bis auf 600.000 Dezimalstellen gut ist. Wenn Sie sich in der Tat ausschließlich für die Berechnung von erf für Extremwerte von interessieren n, dann gehört diese Frage zu math.SE, nicht hier.
Whuber
6
Für "absurd" große Werte werden Sie nichts Besseres tun, als die Obergrenze Pr(Z>z)(z2π)1ez2/2 für doppeltgenaue Gleitkomma-. Diese Annäherung und andere werden an anderer Stelle in der Statistik diskutiert.
Kardinal
Danke Kardinal, diese Grenze scheint ziemlich genau zu sein. Warum machst du das nicht zu einer Antwort?
Douglas S. Stones
@Douglas: Wenn du immer noch interessiert bist, kann ich am nächsten Tag etwas zusammenstellen und es als vollständigere Antwort posten.
Kardinal
1
Nun ... ich denke, es lohnt sich, es als Antwort hinzuzufügen. Vielleicht ist die Grenze in prob + stats allgemein bekannt, aber ich wusste es nicht. Auch sind die Fragen und Antworten hier nicht nur für das OP bestimmt.
Douglas S. Stones

Antworten:

19

Die Frage betrifft die komplementäre Fehlerfunktion

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

für "große" Werte von x ( =n/2 in der ursprünglichen Frage) - dh zwischen 100 und 700.000 oder so. (In der Praxis sollte jeder Wert größer als ungefähr 6 als "groß" betrachtet werden, wie wir sehen werden.) Da dies zur Berechnung von p-Werten verwendet wird, ist es wenig sinnvoll, mehr als drei signifikante (Dezimal-) Stellen zu erhalten .

Betrachten Sie zunächst die von @Iterator vorgeschlagene Annäherung.

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

wo

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

Obwohl dies eine ausgezeichnete Annäherung an die Fehlerfunktion selbst ist, ist es eine schreckliche Annäherung an . Es gibt jedoch eine Möglichkeit, dies systematisch zu beheben.erfc

Für die p-Werte, die mit so großen Werten von assoziiert sind, interessieren wir uns für den relativen Fehler f ( x ) / erfc ( x ) - 1 : Wir hoffen, dass sein absoluter Wert für drei signifikante Stellen der Genauigkeit kleiner als 0,001 ist. Leider ist dieser Ausdruck aufgrund von Unterläufen bei der Berechnung mit doppelter Genauigkeit für große x schwer zu untersuchen . Hier ist ein Versuch, der den relativen Fehler gegen x für 0 x 5,8 darstellt :x f(x)/erfc(x)1xx0x5.8

Plot 1

Die Berechnung wird instabil, sobald 5,3 oder mehr überschreitet, und kann keine signifikante Ziffer mehr nach 5,8 liefern. Dies ist keine Überraschung: exp ( - 5.8 2 ) 10 - 14.6 stößt an die Grenzen der Arithmetik mit doppelter Genauigkeit. Da es keinen Beweis dafür gibt, dass der relative Fehler für ein größeres x akzeptabel klein sein wird , müssen wir es besser machen.xexp(5.82)1014.6x

Das Durchführen der Berechnung in erweiterter Arithmetik (mit Mathematica ) verbessert unser Bild von dem, was vor sich geht:

Plot 2

Der Fehler steigt mit rapide an und zeigt keine Anzeichen einer Nivellierung. Nach x = 10 liefert diese Näherung nicht einmal eine zuverlässige Ziffer an Informationen!xx=10

Die Darstellung beginnt jedoch linear zu wirken. Wir könnten vermuten, dass der relative Fehler direkt proportional zu . (Dies ist aus theoretischen Gründen sinnvoll: erfc ist offensichtlich eine ungerade Funktion und f ist offensichtlich gerade, daher sollte ihr Verhältnis eine ungerade Funktion sein. Daher würden wir erwarten, dass sich der relative Fehler, wenn er zunimmt, wie eine ungerade Potenz von x verhält .) Dies führt uns dazu, den relativen Fehler geteilt durch x zu untersuchen . Gleichermaßen entscheide ich mich, x erfc ( x ) / f ( x ) zu untersuchen.xerfcfx xxerfc(x)/f(x), denn die Hoffnung ist, dass dies einen konstanten Grenzwert haben sollte. Hier ist seine Grafik:

Plot 3

Unsere Vermutung scheint sich zu bestätigen: Dieses Verhältnis scheint sich einer Grenze von ungefähr 8 zu nähern. Auf Anfrage liefert Mathematica Folgendes:

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

Der Wert ist . Dies ermöglicht es uns, die Schätzung zu verbessern:Wir nehmena1=2πe3(4+π)28(3+π)7.94325

f1(x)=f(x)a1x

als erste Verfeinerung der Annäherung. Wenn wirklich groß ist - größer als ein paar Tausend - ist diese Annäherung in Ordnung. Da es für eine Reihe interessanter Argumente zwischen 5.3 und 2000 immer noch nicht gut genug ist , lassen Sie uns die Prozedur wiederholen. Diesmal sollte sich der inverse relative Fehler - insbesondere der Ausdruck 1 - erfc ( x ) / f 1 ( x ) - wie 1 / x 2 für großes x verhalten (aufgrund der vorherigen Paritätsüberlegungen). Dementsprechend multiplizieren wir mit x 2x5.320001erfc(x)/f1(x)1/x2xx2 und finde das nächste Limit:

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

The value is

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

This process can proceed as long as we like. I took it out one more step, finding

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

with value approximately 1623.67. (The full expression involves a degree-eight rational function of π and is too long to be useful here.)

Unwinding these operations yields our final approximation

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

x6x6(1erfc(x)/f3(x)):

Plot 4

f3erfc(x)2661/x6x>0x exceeds 20 or so, we have our three significant digits (or far more, as x gets larger). As a check, here is a table comparing the correct values to the approximation for x between 10 and 20:

 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

In fact, this approximation delivers at least two significant figures of precision for x=8 on, which is just about where pedestrian calculations (such as Excel's NormSDist function) peter out.

Finally, one might worry about our ability to compute the initial approximation f. However, that's not hard: when x is large enough to cause underflows in the exponential, the square root is well approximated by half the exponential,

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

Computing the logarithm of this (in base 10) is simple, and readily gives the desired result. For example, let x=1000. The common logarithm of this approximation is

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

Exponentiating yields

f(1000)2.3416910434296.

Applying the correction (in f3) produces

erfc(1000)1.86003 70486 3232810434298.

Note that the correction reduces the original approximation by over 99% (and indeed, a1/x1%.) (This approximation differs from the correct value only in the last digit. Another well-known approximation, exp(x2)/(xπ), equals 1.86003810434298, erring in the sixth significant digit. I'm sure we could improve that one, too, if we wanted, using the same techniques.)

whuber
quelle
1
+1 This is a great answer, somehow I have never come across this thread before.
amoeba says Reinstate Monica
15

A simple upper bound

For very large values of the argument in the calculation of upper tail probability of a normal, excellent bounds exist that are probably as good as one will get using any other methods with double-precision floating point. For z>0, let

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
quelle
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
quelle
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