Jak obliczyć prawdopodobieństwo związane z absurdalnie dużymi wynikami Z?

Pakiety oprogramowania do wykrywania motywów sieciowych mogą zwracać niezwykle wysokie wyniki Z (najwyższy, jaki widziałem, to 600 000+, ale wyniki Z powyżej 100 są dość powszechne). Planuję pokazać, że te wyniki Z są fałszywe.

Ogromne wyniki Z odpowiadają bardzo niskim związanym prawdopodobieństwom. Wartości powiązanych prawdopodobieństw podano np. Na stronie wikipedii o normalnym rozkładzie (i prawdopodobnie w każdym podręczniku statystyk) dla Z-score do 6. Więc ...

Pytanie : Jak obliczyć funkcję błędu dla n do 1 000 000, powiedzmy?$1-\mathrm{erf}(n/\sqrt{2})$

W szczególności szukam już zaimplementowanego pakietu (jeśli to możliwe). Najlepsze, jakie do tej pory znalazłem, to WolframAlpha, któremu udało się go obliczyć dla n = 150 ( tutaj ).

Pytanie dotyczy uzupełniającej funkcji błędu

$$\textrm{erfc}(x) = \frac{2}{\sqrt{\pi}}\int_{x}^{\infty}\exp(-t^2) d t$$

dla „dużych” wartości ( w pierwotnym pytaniu) - to znaczy między 100 a 700 000 lub mniej więcej. (W praktyce jakakolwiek wartość większa niż około 6 powinna być uważana za „dużą”). Zauważ, że ponieważ zostanie ona wykorzystana do obliczenia wartości p, uzyskanie wartości większej niż trzy cyfry znaczące (dziesiętne) jest niewielkie. .= n / $x$$=n/\sqrt{2}$

Na początek rozważ przybliżenie sugerowane przez @Iterator,

$$f(x) = 1 - \sqrt{1 - \exp \left(-x^2 \left(\frac{4 + ax^2}{\pi + ax^2}\right)\right)},$$

gdzie

$$a = \frac{8(\pi-3)}{3(4-\pi)} \approx 0.439862.$$

Chociaż jest to doskonałe przybliżenie samej funkcji błędu, jest to straszne przybliżenie do . Istnieje jednak sposób, aby to systematycznie naprawić.$\textrm{erfc}$

Dla wartości p związanych z tak dużymi wartościami interesuje nas błąd względny : mamy nadzieję, że jego wartość bezwzględna byłaby mniejsza niż 0,001 dla trzech znaczących cyfry precyzji. Niestety tego wyrażenia trudno jest badać dla dużych ze względu na niedomiar w obliczeniach o podwójnej precyzji. Oto jedna próba, która wykreśla błąd względny względem dla :f ( x ) / erfc ( x ) - 1 x x 0 ≤ x ≤ $x$ $f(x)/\textrm{erfc}(x) - 1$$x$$x$$0 \le x \le 5.8$

Obliczenia stają się niestabilne, gdy przekroczy 5,3 lub mniej więcej i nie może dostarczyć jednej znaczącej cyfry po 5.8. Nie jest to zaskoczeniem: przesuwa granice arytmetyki podwójnej precyzji. Ponieważ nie ma dowodów na to, że błąd względny będzie akceptowalnie mały dla większego , musimy to zrobić lepiej.exp ( - 5,8 2 ) ≈ 10 - 14,6$x$$\exp(-5.8^2) \approx 10^{-14.6}$$x$

Wykonywanie obliczeń w rozszerzonej arytmetyki (z Mathematica ) poprawia nasz obraz tego, co się dzieje:

Błąd rośnie gwałtownie za pomocą i nie wykazuje żadnych oznak wyrównywania. Po lub więcej, to przybliżenie nie zapewnia nawet jednej wiarygodnej cyfry informacji!$x$$x=10$

Fabuła zaczyna jednak wyglądać liniowo. Możemy zgadywać, że błąd względny jest wprost proporcjonalny do . (Ma to uzasadnienie teoretyczne: jest oczywiście funkcją nieparzystą, a jest oczywiście parzystą, więc ich stosunek powinien być funkcją nieparzystą. W związku z tym spodziewalibyśmy się, że błąd względny, jeśli wzrośnie, będzie zachowywał się jak nieparzysta moc .) To prowadzi nas do zbadania błędu względnego podzielonego przez . postanowiłem zbadać , ponieważ istnieje nadzieja, że ​​powinna mieć stałą wartość graniczną. Oto jego wykres:erfc f x x x ⋅ erfc ( x ) / f ( x $x$$\textrm{erfc}$$f$$x$ $x$$x \cdot \textrm{erfc}(x)/f(x)$

Nasze przypuszczenia wydają się potwierdzone: wydaje się, że stosunek ten zbliża się do granicy około 8. Zapytany, Mathematica dostarczy go:

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

Wartość wynosi . . To pozwala nam poprawić szacunek: bierzemy$a_1 = \frac{2 }{\sqrt{\pi }}e^{\frac{3 (-4+\pi )^2}{8 (-3+\pi )}} \approx 7.94325$

$$f_1(x) = f(x) \frac{a_1}{x}$$

jako pierwsze udoskonalenie przybliżenia. Gdy jest naprawdę duży - większy niż kilka tysięcy - to przybliżenie jest w porządku. Ponieważ nadal nie będzie wystarczająco dobry dla interesującego zakresu argumentów między a lub mniej, powtórzmy procedurę. Tym razem odwrotny błąd względny - w szczególności wyrażenie powinien zachowywać się jak dla dużych (na podstawie poprzednich rozważań dotyczących parzystości) . W związku z tym mnożymy przez i znajdujemy następny limit:5,3 2000 1 - erfc ( x ) / f 1 ( x ) 1 / x 2 x x $x$$5.3$$2000$$1 - \textrm{erfc}(x)/f_1(x)$$1/x^2$$x$$x^2$

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

Wartość wynosi

$$a_2 = \frac{1}{32 \sqrt{\pi }} e^{\frac{3 (-4+\pi )^2}{8 (-3+\pi )}} \left(32-\frac{9 (-4+\pi )^3 \pi }{(-3+\pi )^2}\right) \approx 114.687.$$

Proces ten może trwać tak długo, jak chcemy. Zrobiłem to jeszcze jeden krok, znajdując

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

o wartości około 1623,67. (Pełne wyrażenie obejmuje racjonalną funkcję stopnia i jest zbyt długa, aby była tu użyteczna.)$\pi$

Odkręcenie tych operacji daje nam ostateczne przybliżenie

$$f_3(x) = f(x)\left(a_1 - a_2/x^2 + a_3/x^4\right)/x.$$

Błąd jest proporcjonalny do . Import jest stałą proporcjonalności, więc wykreślamy : x 6 ( 1 - erfc ( x ) / f 3 ( x ) $x^{-6}$$x^6(1 - \textrm{erfc}(x) / f_3(x))$

Szybko zbliża się do wartości granicznej około 2660,59. Korzystając z aproksymacji , otrzymujemy oszacowania którego względna dokładność jest lepsza niż dla wszystkich . Gdy przekroczy około 20, mamy trzy cyfry znaczące (lub znacznie więcej, ponieważ staje się większy). Jako sprawdzenie, oto tabela porównująca prawidłowe wartości z przybliżeniem dla między a :erfc ( $f_3$2661 / x 6 x > 0 x x x 10 $\textrm{erfc}(x)$$2661/x^6$$x \gt 0$$x$$x$$x$$10$$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

W rzeczywistości, to przybliżenie dostarcza co najmniej dwie znaczące liczby precyzji dla , co oznacza, że ​​obliczenia pieszych (takie jak funkcja Excela ) się kończą.$x=8$NormSDist

Na koniec można się martwić naszą zdolnością do obliczenia wstępnego przybliżenia . Nie jest to jednak trudne: gdy jest wystarczająco duży, aby spowodować niedomiar wykładniczy, pierwiastek kwadratowy jest dobrze przybliżony do połowy wykładniczej,$f$$x$

$$f(x) \approx \frac{1}{2} \exp \left(-x^2 \left(\frac{4 + ax^2}{\pi + ax^2}\right)\right).$$

Obliczenie tego logarytmu (w podstawie 10) jest proste i łatwo daje pożądany rezultat. Na przykład niech . Wspólnym logarytmem tego przybliżenia jest$x=1000$

$$\log_{10}(f(x)) \approx \left(-1000^2 \left(\frac{4 + a \cdot 1000^2}{\pi + a \cdot 1000^2}\right) - \log(2)\right) / \log(10) \sim -434295.63047.$$

Wykładnicze wydajności

$$f(1000) \approx 2.34169 \cdot 10^{-434296}.$$

Zastosowanie poprawki (w ) daje wynik$f_3$

$$\textrm{erfc}(1000) \approx 1.86003\ 70486\ 32328 \cdot 10^{-434298}.$$

Zwróć uwagę, że poprawka zmniejsza pierwotne przybliżenie o ponad 99% (i rzeczywiście ). (To przybliżenie różni się od poprawnej wartości tylko w ostatniej cyfrze. Kolejne dobrze znane przybliżenie, , równa się , z szóstej cyfry znaczącej. Jestem pewien, że moglibyśmy to poprawić, jeśli poszukiwany przy użyciu tych samych technik).$a_1/x \approx 1\text{%}$1,860038⋅10 - $\exp(-x^2)/(x\sqrt{\pi})$$1.860038 \cdot 10^{-434298}$

Prosta górna granica

Dla bardzo dużych wartości argumentu przy obliczaniu prawdopodobieństwa górnego ogona normalnej istnieją doskonałe granice, które prawdopodobnie są tak dobre, jak można je uzyskać przy użyciu innych metod z zmiennoprzecinkową podwójnej precyzji. Dla , niech gdzie jest standardowym normalnym plikiem pdf. Użyłem notacji w odniesieniu do standardowej notacji w analizie przeżycia. W kontekście inżynierii nazywają tę funkcję funkcją i oznaczają ją przez .S ( z ) $z > 0$Φ ( z ) = ( 2 π ) - 1 / 2 e -

$$\renewcommand{\Pr}{\mathbb{P}}\newcommand{\rd}{\mathrm{d}} S(z) := \Pr(Z > z) = \int_z^\infty \varphi(z) \rd z \>,$$ $\varphi(z) = (2\pi)^{-1/2} e^{-z^2/2}$$S(z)$$Q$$Q(z)$

Zatem bardzo prostą, podstawową górną granicą jest gdzie notacja po prawej stronie wskazuje, że jest to oszacowanie górnej granicy. Ta odpowiedź stanowi dowód na związanie.

$$\newcommand{\Su}{\hat{S}_u} \newcommand{\Sl}{\hat{S}_\ell} S(z) \leq \frac{\varphi(z)}{z} =: \Su(z) \> ,$$

Istnieje również kilka miłych, uzupełniających się dolnych granic. Jednym z najłatwiejszych i najłatwiejszych do uzyskania jest związany Istnieją co najmniej trzy oddzielne metody wyprowadzenia tego ograniczenia. W tej odpowiedzi na powiązane pytanie można znaleźć przybliżony szkic jednej z takich metod .

$$S(z) \geq \frac{z}{z^2+1} \varphi(z) =: \Sl(z) \> .$$

Obrazek

Poniżej znajduje się wykres dwóch granic (w kolorze szarym) wraz z rzeczywistą funkcją .$S(z)$

Jak to jest dobre

Z fabuły wydaje się, że granice stają się dość ciasne, nawet dla umiarkowanie dużego . Możemy zadać sobie pytanie, jak ciasne są i jakie ilościowe stwierdzenie w tym względzie można sformułować.$z$

Jedną przydatną miarą szczelności jest bezwzględny błąd względny To daje błąd proporcjonalny oszacowania.

$$\newcommand{\err}{\mathcal{E}} \err(z) = \left|\frac{\Su(z) - S(z)}{S(z)}\right| \>.$$

Teraz zauważ, że ponieważ wszystkie zaangażowane funkcje są nieujemne, przez użycie właściwości ograniczających i , otrzymujemy a zatem stanowi to dowód że dla górna granica jest poprawna do 1%, dla jest poprawna z dokładnością do 0,1%, a dla z dokładnością do 0,01%.$\Su(z)$$\Sl(z)$

$$\err(z) = \frac{\Su(z) - S(z)}{S(z)} \leq \frac{\Su(z) - \Sl(z)}{\Sl(z)} = z^{-2} \> ,$$ $z \geq 10$$z \geq 28$$z \geq 100$

W rzeczywistości prosta forma granic zapewnia dobrą kontrolę innych „przybliżeń”. Jeśli w obliczeniach numerycznych bardziej skomplikowanych aproksymacji otrzymamy wartość poza tymi granicami, możemy po prostu „poprawić” ją, przyjmując wartość np. Podanej tutaj górnej granicy.

Istnieje wiele udoskonaleń tych granic. Wspomniane tutaj granice Laplace'a zapewniają ładną sekwencję górnych i dolnych granic na postaci gdzie jest funkcją wymierną.R ( z ) φ ( z ) R ( z $S(z)$$R(z) \varphi(z)$$R(z)$

Wreszcie, oto kolejne nieco związane pytanie i odpowiedź.

Możesz to przybliżyć za pomocą znacznie prostszych funkcji - więcej informacji można znaleźć w tej sekcji Wikipedii . Podstawowym przybliżeniem jest to, że$\textrm{erf}(x) \approx \textrm{sgn}(x)\sqrt{1 - \exp(-x^2 \frac{4/\pi + ax^2}{1+ax^2}})$

Artykuł zawiera niepoprawny link do tej sekcji. Odnośny plik PDF można znaleźć w plikach Siergieja Winitzkiego - lub pod tym linkiem .