Jak mogę (liczbowo) aproksymować wartości dla rozkładu beta z dużymi wartościami alfa i beta

12

Czy istnieje stabilny numerycznie sposób obliczania wartości rozkładu beta dla dużej liczby całkowitej alfa, beta (np. Alfa, beta> 1000000)?

Właściwie potrzebuję tylko 99% przedziału ufności wokół trybu, jeśli to w jakiś sposób ułatwi problem.

Dodaj : Przepraszam, moje pytanie nie było tak jasno określone, jak myślałem. Chcę to zrobić: mam maszynę, która sprawdza produkty na przenośniku taśmowym. Część frakcji tych produktów jest odrzucana przez maszynę. Teraz, jeśli operator maszyny zmieni jakieś ustawienia inspekcji, chcę mu pokazać szacunkową częstotliwość odrzucania i podpowiedź na temat tego, jak wiarygodne jest bieżące oszacowanie.

Pomyślałem więc, że traktuję rzeczywistą częstotliwość odrzucania jako zmienną losową X i obliczam rozkład prawdopodobieństwa dla tej zmiennej losowej na podstawie liczby odrzuconych obiektów N i zaakceptowanych obiektów M. Jeśli założę jednolity wcześniejszy rozkład dla X, jest to rozkład beta w zależności od N i M. Mogę albo wyświetlić ten rozkład użytkownikowi bezpośrednio, albo znaleźć przedział [l, r], aby rzeczywista częstość odrzucania była w tym przedziale przy p> = 0,99 (używając terminologii shabbychef) i wyświetlić to interwał. Dla małych M, N (tj. Bezpośrednio po zmianie parametru) mogę obliczyć rozkład bezpośrednio i przybliżać przedział [l, r]. Ale w przypadku dużych M, N takie naiwne podejście prowadzi do błędów niedopełnienia, ponieważ x ^ N * (1-x) ^ M jest zbyt małe, aby można je było przedstawić jako zmiennoprzecinkowe podwójnej precyzji.

Myślę, że moim najlepszym wyborem jest użycie mojej naiwnej dystrybucji beta dla małych M, N i przejście do normalnej dystrybucji z tą samą średnią i wariancją, gdy tylko M, N przekroczy pewien próg. Czy to ma sens?

nikie
źródło
1
Czy chcesz poznać matematykę, czy po prostu rozwiązanie kodu w R lub coś takiego?
John
Muszę zaimplementować to w języku C #, aby matematyka była dobra. Próbka kodu też byłaby w porządku, jeśli nie opiera się na jakiejś wbudowanej funkcji R / Matlab / Mathematica, której nie mogę przetłumaczyć na C #.
nikie
PDF, CDF czy odwrotny CDF?
JM nie jest statystykiem
Jeśli nie nalegasz na wersję beta, możesz użyć dystrybucji Kumaraswamy, która jest bardzo podobna i ma znacznie prostszą formę algebraiczną: en.wikipedia.org/wiki/Kumaraswamy_distribution
Tim

Odpowiedzi:

13

Normalne przybliżenie działa wyjątkowo dobrze, szczególnie w ogonach. Użyj średniej z i wariantu . Na przykład bezwzględny błąd względny w prawdopodobieństwie ogona w trudnej sytuacji (gdzie skośność może budzić obawy), taki jak osiąga około i jest mniejszy niż gdy jesteś więcej niż 1 SD od średniej. ( Nie dzieje się tak, ponieważ beta jest tak duża: przy bezwzględne błędy względne są ograniczone przezα/(α+β)αβ(α+β)2(1+α+β)α=106,β=1080.000260.00006α=β=1060.0000001.) To przybliżenie jest zatem doskonałe do zasadniczo dowolnego celu obejmującego przedziały 99%.

W świetle zmian wprowadzonych do pytania należy zauważyć, że nie oblicza się całek beta poprzez faktyczne całkowanie całki: oczywiście dostaniesz niedopełnienia (chociaż tak naprawdę nie mają znaczenia, ponieważ nie mają znaczącego wpływu na całkę) . Istnieje wiele sposobów obliczania całki lub przybliżania jej, jak udokumentowano w Johnson & Kotz (Distribution in Statistics). Kalkulator online znajduje się na stronie http://www.danielsoper.com/statcalc/calc37.aspx . W rzeczywistości potrzebujesz odwrotności tej całki. Niektóre metody obliczania odwrotności są udokumentowane na stronie Mathematica pod adresem http://functions.wolfram.com/GammaBetaErf/InverseBetaRegularized/. Kod znajduje się w przepisach numerycznych (www.nr.com). Naprawdę fajnym kalkulatorem online jest strona Wolfram Alpha (www.wolframalpha.com): wprowadź inverse beta regularized (.005, 1000000, 1000001)lewy punkt końcowy i inverse beta regularized (.995, 1000000, 1000001)prawy punkt końcowy ( , przedział 99%).α=1000000,β=1000001

Whuber
źródło
Idealny! Cały czas miałem na biurku książkę z NR, ale nigdy nie pomyślałem, żeby tam zajrzeć. Wielkie dzięki.
nikie
3

Szybki eksperyment graficzny sugeruje, że rozkład beta wygląda bardzo podobnie do rozkładu normalnego, gdy alfa i beta są bardzo duże. Po przejrzeniu „normalnego limitu dystrybucji beta” znalazłem http://nrich.maths.org/discus/messages/117730/143065.html?1200700623 , co daje „dowód” na machanie ręką.

Strona wikipedii dla dystrybucji beta podaje jej średnią, tryb (v zbliżony do średniej dla dużej alfa i beta) i wariancję, więc można użyć rozkładu normalnego z tą samą średnią i wariancją, aby uzyskać przybliżenie. To, czy jest to wystarczająco dobre przybliżenie dla twoich celów, zależy od twoich celów.

jeden przystanek
źródło
Głupie pytanie: jak przeprowadziłeś ten eksperyment graficzny? Próbowałem wykreślić rozkład dla wersji alfa / beta na około 100, ale nic nie widziałem z powodu błędów niedopełnienia.
nikie
Nie chcesz wykreślić całki: chcesz wykreślić całkę. Można jednak uzyskać integrand na wiele sposobów. Jednym z nich jest wprowadzenie „wykresu D (beta (x, 1000000, 2000000), x) / beta (1, 1000000, 2000000) od 0,3325 do 0,334” na stronie Wolfram Alpha. Sama całka jest widoczna w „Wykresie beta (x, 1000000, 2000000) / beta (1, 1000000, 2000000) od 0,3325 do 0,334”.
whuber
Narysowałem integrand, tj. Pdf dystrybucji beta, w Stata - ma wbudowaną funkcję dla pdf. W przypadku dużych wersji alfa i beta musisz ograniczyć zasięg wykresu, aby był bliski normalności. Gdybym sam go programował, obliczałbym jego logarytm, a następnie wykładnik potęgował się na końcu. To powinno pomóc w problemach niedomiaru. Funkcja beta w mianowniku jest zdefiniowana w kategoriach funkcji gamma, równoważnych silniom dla liczb całkowitych alfa i beta, a wiele pakietów / bibliotek zawiera lngamma () lub lnfactorial () zamiast / a także funkcje gamma () i silnia ().
onestop
2

Mam zamiar wywnioskować, że chcesz przedział taki, że prawdopodobieństwo losowego losowania z Beta RV jest w przedziale z prawdopodobieństwem 0,99, z dodatkowymi punktami dla i symetrycznymi wokół trybu. Na podstawie Nierówności Gaussa lub nierówności Vysochanskii-Petunin można konstruować przedziały, które zawierają przedział , i byłyby to całkiem przyzwoite przybliżenia. W przypadku dostatecznie dużych będziesz miał problemy z niedopełnieniem liczb, nawet reprezentując i jako odrębne liczby, więc ta droga może być wystarczająco dobra.[l,r]lr[l,r]α,β lr

shabbychef
źródło
Gdy alfa i beta nie są zbyt daleko od siebie (tj. Alfa / beta są ograniczone powyżej i poniżej), SD beta (alfa, beta) jest proporcjonalne do 1 / Sqrt (alfa). Np. Dla alpha = beta = 10 ^ 6, SD jest bardzo bliskie 1 / Sqrt (8) / 1000. Myślę, że nie będzie problemu z reprezentacją lir, nawet jeśli używasz tylko pływaków pojedynczej precyzji .
whuber
co oznacza, że nie jest „wystarczająco duże”;)106
shabbychef
1
Tak, to szalona liczba dla wersji beta. BTW, te nierówności w ogóle nie spowodują dobrych odstępów, ponieważ są skrajne w stosunku do wszystkich rozkładów (spełniając pewne ograniczenia).
whuber
@whuber: Masz rację, to szalone liczby. Dzięki mojemu naiwnemu algorytmowi „rozsądne” liczby były łatwe i działały dobrze, ale nie wyobrażałem sobie, jak je obliczyć dla „szalonych” parametrów. Stąd pytanie.
nikie
2
OK, masz rację: gdy alfa + beta przekroczy 10 ^ 30 lub więcej, będziesz mieć trudności z podwójnymi :-). (Ale jeśli reprezentujesz l i r jako różnice w stosunku do średniej alfa / (alfa + beta), nic ci nie będzie, dopóki alfa lub beta nie przekroczą około 10 ^ 303.)
whuber
1

Jeśli jest zmienną rozproszoną w wersji beta, to logarytmiczna szansa (tj. jest w przybliżeniu normalnie rozłożona. Jest to prawdą nawet w przypadku mocno wypaczonych rozkładów beta wraz zp l o g ( p / ( 1 - p ) ) m i n ( α , β ) > 100pplog(p/(1p))min(α,β)>100

Na przykład

f <- function(n, a, b) {
    p <- rbeta(n, a, b)
    lor <- log(p/(1-p))
    ks.test(lor, 'pnorm', mean(lor), sd(lor))$p.value
}
summary(replicate(50, f(10000, 100, 1000000)))

zazwyczaj daje wynik podobny do

podsumowanie (powtórzenie (50, f (10000, 100, 1000000))) min. 1st Qu. Mediana oznacza 3. kwartę Max. 0,01205 0,108 0,08680 0,248 0,36670 0,68730

tzn. typowe wartości p wynoszą około 0,2.

Tak więc nawet przy 10000 próbkach test Kołmogorowa-Smirnowa nie ma mocy, aby odróżnić transformację ilorazu logarytmicznego silnie wypaczonej zmiennej rozproszonej beta o .α=100,β=100000

Jednak podobny test rozkładu samegop

f2 <- function(n, a, b) {
    p <- rbeta(n, a, b)
    ks.test(p, 'pnorm', mean(p), sd(p))$p.value
}
summary(replicate(50, f2(10000, 100, 1000000)))

produkuje coś podobnego

summary(replicate(50, f2(10000, 100, 1000000)))
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
2.462e-05 3.156e-03 7.614e-03 1.780e-02 1.699e-02 2.280e-01 

z typowymi wartościami p około 0,01

Funkcja R qqnormdaje również pomocną wizualizację, tworząc bardzo prosto wyglądający wykres rozkładu logarytmiczno-szansowego wskazujący przybliżoną normalność, rozkład zmiennej dsitribute beta daje charakterystyczną krzywą wskazującą na nienormalność

Dlatego rozsądne jest zastosowanie aproksymacji Gaussa w przestrzeni logarytmicznej, nawet dla mocno wypaczonych wartości o ile oba są większe niż 100.α,β

Daniel Mahler
źródło