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

14

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

Douglas S. Kamienie
źródło
6
Może to nie jest właściwe pytanie. Te oceny Z są fałszywe, ponieważ zakładają, że rozkład normalny jest znacznie lepszym przybliżeniem lub modelem, niż jest w rzeczywistości. To trochę tak, jakby zakładać, że mechanika Newtona jest dobra do 600 000 miejsc po przecinku. Jeśli rzeczywiście jesteś zainteresowany wyłącznie obliczeniem erf dla ekstremalnych wartości , to pytanie należy do matematyki.SE, a nie tutaj. n
whuber
6
W przypadku „absurdalnie” dużych wartości nie da się lepiej niż górna granica dla zmiennoprzecinkowego podwójnej precyzji. To przybliżenie i inne są omówione w innym miejscu na stats.SE. Pr(Z>z)(z2π)1ez2/2
kardynał
Dzięki kardynałowi ta granica wydaje się dość dokładna. Dlaczego nie uczynisz tego odpowiedzią?
Douglas S. Stones
@Douglas: Jeśli nadal jesteś zainteresowany, mogę złożyć coś następnego dnia i opublikować jako pełniejszą odpowiedź.
kardynał
1
Cóż ... Myślę, że warto byłoby dodać to jako odpowiedź. Może granica jest powszechną wiedzą w statystykach prob +, ale nie wiedziałem o tym. Ponadto, Q i A tutaj nie są wyłącznie dla PO.
Douglas S. Stones

Odpowiedzi:

19

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

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

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/2

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

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

gdzie

a=8(π3)3(4π)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ć.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 5,8x f(x)/erfc(x)1xx0x5.8

Wykres 1

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

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

Wykres 2

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!xx=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 )xerfcfx xxerfc(x)/f(x)

Wykres 3

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: bierzemya1=2πe3(4+π)28(3+π)7.94325

f1(x)=f(x)a1x

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 2x5.320001erfc(x)/f1(x)1/x2xx2

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

Wartość wynosi

a2=132πe3(4+π)28(3+π)(329(4+π)3π(3+π)2)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.)π

Odkręcenie tych operacji daje nam ostateczne przybliżenie

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

Błąd jest proporcjonalny do . Import jest stałą proporcjonalności, więc wykreślamy : x 6 ( 1 - erfc ( x ) / f 3 ( x ) )x6x6(1erfc(x)/f3(x))

Wykres 4

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 ( xf32661 / x 6 x > 0 x x x 10 20erfc(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

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

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

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

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

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

Wykładnicze wydajności

f(1000)2.3416910434296.

Zastosowanie poprawki (w ) daje wynikf3

erfc(1000)1.86003 70486 3232810434298.

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).a1/x1%1,86003810 - 434298exp(x2)/(xπ)1.86003810434298

Whuber
źródło
1
+1 To świetna odpowiedź, jakoś nigdy wcześniej nie spotkałem tego wątku.
ameba mówi Przywróć Monikę
15

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

S(z):=P(Z>z)=zφ(z)dz,
φ(z)=(2π)1/2ez2/2S(z)QQ(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.

S(z)φ(z)z=:S^u(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)zz2+1φ(z)=:S^(z).

Obrazek

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

Górny ogon normy i granic

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.

E(z)=|S^u(z)S(z)S(z)|.

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%.S^u(z)S^(z)

E(z)=S^u(z)S(z)S(z)S^u(z)S^(z)S^(z)=z2,
z10z28z100

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)φ(z)R(z)

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

kardynał
źródło
1
Przepraszamy za wszystkie „cytaty własne”. Raz, kilka lat temu, intensywnie, dwutygodniowo interesowałem się pokrewnymi pytaniami i starałem się dowiedzieć jak najwięcej na ten temat.
kardynał
+1 Zgadzam się z jacuzzi. Bardzo miło i doceniam linki do innych odpowiedzi.
Iterator
5

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, żeerf(x)sgn(x)1exp(x24/π+ax21+ax2)

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

Iterator
źródło
1
Pewne wzmocnienie tego byłoby mile widziane z dwóch powodów. Po pierwsze, najlepiej, gdy odpowiedzi mogą być samodzielne. Po drugie, artykuł ten dwuznacznie pisze o jakości przybliżenia „w sąsiedztwie nieskończoności”: jak dokładna jest „bardzo dokładna”? (W domyśle masz to dobre pojęcie, ale od wszystkich zainteresowanych czytelników jest wiele.) Podana wartość „.00035” jest tutaj bezużyteczna.
whuber
Dzięki. Nie zauważyłem, że istnieje wsparcie dla Javascript w korzystaniu z TeXa, co zrobiło różnicę w napisaniu tego.
Iterator,
1
Nawiasem mówiąc, odwołanie Wikipedii do tego przybliżenia jest zepsute. Mathematica stwierdza jednak, że błąd względny (1 - około (x) / erf (x)) zachowuje się jak odwrotność . 2exp(x2+3(π4)2/(8(π3)))
whuber
@whuber, czy możesz w tym celu napisać kod Mathematica? :) Nie widziałem Mathematiki od ponad 15 lat i nigdy do tego celu.
Iterator
Opublikowałem to w osobnej odpowiedzi.
whuber