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?
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 ).
probability
normal-distribution
p-value
approximation
z-statistic
Douglas S. Kamienie
źródło
źródło
Odpowiedzi:
Pytanie dotyczy uzupełniającej funkcji błędu
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,
gdzie
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)−1 x x 0≤x≤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 xx exp(−5.82)≈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 erfc f x x x⋅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:
Wartość wynosi . . To pozwala nam poprawić szacunek: bierzemya1=2π√e3(−4+π)28(−3+π)≈7.94325
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 2x 5.3 2000 1−erfc(x)/f1(x) 1/x2 x x2
Wartość wynosi
Proces ten może trwać tak długo, jak chcemy. Zrobiłem to jeszcze jeden krok, znajdując
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
Błąd jest proporcjonalny do . Import jest stałą proporcjonalności, więc wykreślamy : x 6 ( 1 - erfc ( x ) / f 3 ( x ) )x- 6 x6( 1 - erfc ( x ) / f3)( 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 ( xfa3) 2661 / x 6 x > 0 x x x 10 20erfc (x) 2661 / x6 x > 0 x x x 10 20
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,fa x
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
Wykładnicze wydajności
Zastosowanie poprawki (w ) daje wynikf3
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/x≈1% 1,860038⋅10 - 434298exp(−x2)/(xπ−−√) 1.860038⋅10−434298
źródło
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
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.
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 .
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.
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)
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ź.
źródło
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)1−exp(−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 .
źródło