Jak mogę skutecznie modelować sumę zmiennych losowych Bernoulliego?

38

Modeluję zmienną losową ( ), która jest sumą około 15-40k niezależnych zmiennych losowych Bernoulliego ( ), z których każda ma inne prawdopodobieństwo powodzenia ( ). Formalnie gdzie i . $Y$ $X_i$ $p_i$ $Y=\sum X_i$ $\Pr(X_i=1)=p_i$ $\Pr(X_i=0)=1-p_i$

Interesuje mnie szybkie odpowiadanie na zapytania, takie jak $\Pr(Y<=k)$ (gdzie podano $k$ ).

Obecnie używam losowych symulacji, aby odpowiedzieć na takie pytania. Losowo rysuję każdy $X_i$ zgodnie z jego $p_i$ , a następnie wszystkie $X_i$ wartości , aby uzyskać $Y'$ . Powtarzam ten proces kilka tysięcy razy i zwracam ułamek razy $\Pr(Y'\leq k)$ .

Oczywiście nie jest to całkowicie dokładne (chociaż dokładność znacznie wzrasta wraz ze wzrostem liczby symulacji). Wydaje się również, że mam wystarczającą ilość danych na temat dystrybucji, aby uniknąć symulacji użytkowania. Czy możesz wymyślić rozsądny sposób, aby uzyskać dokładne prawdopodobieństwo $\Pr(Y\leq k)$ ?

ps

Używam Perla i R.

EDYTOWAĆ

Po odpowiedziach pomyślałem, że mogą być potrzebne pewne wyjaśnienia. Wkrótce opiszę ustawienie mojego problemu. Podano genom kołowy z obwodem ci nmapowanym do niego zestawem zakresów. Na przykład c=3*10^9i ranges={[100,200],[50,1000],[3*10^9-1,1000],...}. Uwaga: wszystkie zakresy są zamknięte (oba końce są włącznie). Pamiętaj również, że mamy do czynienia tylko z liczbami całkowitymi (całe jednostki).

Szukam regionów na okręgu, które są niedokryte przez podane nzakresy mapowane. Aby sprawdzić, czy dany zakres długości xna okręgu jest ukryty, testuję hipotezę, że nzakresy są mapowane losowo. Prawdopodobieństwo, że mapowany zakres długości w q>xpełni obejmie podany zakres długości, xwynosi (q-x)/c. Prawdopodobieństwo to staje się dość małe, gdy cjest duże i / lub qjest małe. Interesuje mnie liczba zakresów (poza n), które obejmują x. Tak Ypowstaje.

Testuję moją hipotezę zerową w porównaniu z jednostronną alternatywą (tajność). Zauważ też, że testuję wiele hipotez (różne xdługości) i na pewno to poprawię.

— David B.
źródło

Czy twoje p_i jest stałe w trakcie ćwiczenia modelowania, czy może zmieniać się z jednego obliczenia do następnego?

— whuber

Te p_isą naprawione.

— David B,

Czy w świetle bieżących odpowiedzi można podzielić szacunki (a) sumy punktów p i (b) sumy ich kwadratów? Te wartości określają opcje.

— whuber

@ whuber: różnią się one znacznie w zależności od przypadku. To nie jest jednorazowy moduł, który tworzę (niestety).

— David B,

@David Ale czy nie możesz podać wskazówek, takich jak typowe zakresy? Na przykład, jeśli suma wartości p waha się między 1 a 100, to jest użyteczna informacja i sugeruje pewne wydajne rozwiązania, ale jeśli można uzyskać do 10 000, może to wykluczyć niektóre podejścia.

— whuber

24

Jeśli często przypomina Poissona , czy próbowałeś przybliżyć go do Poissona za pomocą parametru ? $\lambda = \sum p_i$

EDYCJA : Znalazłem teoretyczny wynik, który to uzasadnia, a także nazwę dla rozkładu : nazywa się to rozkładem dwumianowym Poissona . Nierówność Le Cam mówi ci, jak blisko jej rozkład jest aproksymowany rozkładem Poissona z parametrem . Mówi ci, że jakość tego ok jest regulowana przez sumę kwadratów , parafrazując Steele (1994) . Więc jeśli wszystkie twoje są dość małe, tak jak się teraz wydaje, powinny być całkiem dobrym przybliżeniem. $Y$ $\lambda = \sum p_i$ $p_i$ $p_i$

EDYCJA 2 : Jak mały jest „dość mały”? To zależy od tego, jak dobre jest to przybliżenie! Artykuł w Wikipedii na temat twierdzenia Le Cam podaje dokładną formę wyniku, o którym wspomniałem powyżej: suma bezwzględnych różnic między funkcją masy prawdopodobieństwa (pmf) i pmf powyższego rozkładu Poissona jest nie większa niż dwukrotność sumy kwadratów s. Kolejny wynik z Le Cam (1960) może być łatwiejszy w użyciu: suma ta jest również nie więcej niż 18 razy większa od . Jest jeszcze kilka takich wyników ... patrz Serfling (1978) dla jednej recenzji. $Y$ $p_i$ $p_i$

— jeden przystanek
źródło

1

+1 Niezły pomysł. Jest prawdopodobne, że niewielka mieszanka Poissonów wykonałaby dobrą robotę, w zależności od wyjaśnienia pytania.

— whuber

1

Zastanawiałem się nad zasugerowaniem ujemnego rozkładu dwumianowego, który powstaje jako mieszanina Gamma-Poissona, ale który ma wariancję większą niż jego średnia, podczas gdy ten problem ma wariancję mniejszą niż jego średnia. W związku z tym nie jestem pewien, czy jakakolwiek mieszanina Poissonów będzie działać, ponieważ na pewno każda taka mieszanina będzie miała wariancję większą niż jej średnia?

— onestop

@onestop Gdzie powiedziano, że wariancja jest mniejsza niż średnia? Brakowało mi tego oświadczenia.

— whuber

Przepraszam, co tam, to było trochę tajemnicze, ale te komentarze nie pozwalają na tak szczegółowe opracowanie. mpiktas męska jest wariancją, która jest mniejsza niż średnia, . Tylko nieznacznie mniej, jeśli są średnio bardzo małe, więc standardowy Poisson może być wystarczająco dobry. Może powinienem rozwinąć moją odpowiedź powyżej ... ale wtedy wątek konwersacyjny staje się mylący.

B_{n} = \sum p_{i} (1 - p_{i})

$B_n = \sum p_i(1-p_i)$

\sum p_{i}

$\sum p_i$

p_{i}

$p_i$

— onestop

Co rozumiesz przez ? Jak uzyskać wartości ?

\sum X_{i}

$\sum X_i$

X_{i}

$X_i$

— David B,

11

Natknąłem się na twoje pytanie, szukając rozwiązania tego bardzo problemu. Odpowiedzi tutaj nie były do końca satysfakcjonujące, ale wydaje mi się, że istnieje dość proste rozwiązanie, które zapewnia dokładną dystrybucję i jest całkiem łatwe do wykonania.

Rozkład sumy dwóch dyskretnych zmiennych losowych jest splotem ich gęstości. Więc jeśli masz gdzie znasz i , możesz obliczyć: $Z = X + Y$ $P(X)$ $P(Y)$

P (Z = z) = \sum_{k = - \infty}^{\infty} P (X = k) P (Y = z - k)

$P(Z=z) = \sum_{k=-\infty}^{\infty} P(X=k) \; P(Y=z-k)$

(Oczywiście dla Bernoulliego zmiennych losowych nie trzeba iść dość do nieskończoności).

Możesz użyć tego, aby znaleźć dokładny rozkład sumy RV. Najpierw zsumuj dwa RV razem, zwijając ich pliki PDF (np. [0,3, 0,7] * [0,6, 0,4] = [0,18, 0,54, 0,28]). Następnie przekonaj tę nową dystrybucję do następnego pliku PDF Bernoulli (np. [0,18, 0,54, 0,28] * [0,5, 0,5] = [0,09, 0,36, 0,41, 0,14]). Powtarzaj tę czynność, aż wszystkie RV zostaną dodane. I voila, wynikowy wektor jest dokładnym plikiem PDF sumy wszystkich twoich zmiennych.

Za pomocą symulacji zweryfikowałem, że daje to prawidłowe wyniki. Nie opiera się na żadnych asymptotycznych założeniach i nie ma wymagań, aby sondy Bernoulliego były małe.

Może być też jakiś sposób, aby to zrobić bardziej efektywnie niż powtarzane splot, ale nie myślałem o tym zbyt głęboko. Mam nadzieję, że to komuś pomoże!

— alex
źródło

2

Czy próbowałeś tego ze zmiennymi 40K? (Zastanawiam się, ile godzin lub dni obliczeń potrzeba ...)

— whuber

5

(+1) Znalazłem sposób, aby ten pomysł zadziałał. Wymaga dwóch technik: po pierwsze, użyj FFT do zwojów; po drugie, nie rób ich sekwencyjnie, ale dziel i rządź: rób je w rozłącznych parach, a następnie rób wyniki w rozłącznych parach itp. Algorytm skaluje się teraz jako zamiast dla prawdopodobieństw. Na przykład Mathematica może obliczyć cały rozkład dla 40 000 prawdopodobieństw w zaledwie 0,4 sekundy. (1 000 000 jest obliczanych w 10,5 sekundy). Podam kod w komentarzu uzupełniającym.

O (n \log n)

$O(n\log n)$

O (n^{2})

$O(n^2)$

n

$n$

— whuber

7

Oto kod Mathematica :

multinomial[p_] := Module[{lc, condense},   lc = Function[{s}, ListConvolve[s[[1]], s[[2]], {1, -1}, 0]];   condense = Function[{s}, Map[lc, Partition[s, 2, 2, {1, 1}, {{1}}]]];   Flatten[NestWhile[condense, Transpose[{1 - p, p}], Length[#] > 1 &]]]

Aby go zastosować, zrób coś takiego p = RandomReal[{0, 1}, 40000]; pp = multinomial[p];. Stwarza to prawdopodobieństwa, pa następnie oblicza dokładny rozkład pp. Uwaga: Kiedy średnia pnie jest ekstremalna, rozkład jest bardzo zbliżony do normalnego: prowadzi to do znacznie szybszego algorytmu.

— whuber

9

@onestop zapewnia dobre referencje. Artykuł w Wikipedii o rozkładzie dwumianowym Poissona podaje rekurencyjną formułę obliczania dokładnego rozkładu prawdopodobieństwa; wymaga wysiłku . Niestety, jest to suma naprzemienna, więc będzie niestabilna numerycznie: obliczenia z arytmetyką zmiennoprzecinkową są beznadziejne. Na szczęście, gdy są małe, wystarczy obliczyć niewielką liczbę prawdopodobieństw, więc wysiłek jest naprawdę proporcjonalny do . Precyzja niezbędna do przeprowadzenia obliczeń za pomocą racjonalnej arytmetyki ( tj. Dokładnie tak , aby niestabilność numeryczna nie stanowiła problemu) rośnie na tyle wolno, że całkowity czas może nadal wynosić w przybliżeniu $O(n^2)$ $p_i$ $O(n \log(\sum_i{p_i}))$ $O(n^2)$ . To wykonalne.

Jako test stworzyłem tablicę prawdopodobieństw dla różnych wartości od do , co jest rozmiarem tego problemu. Dla małych wartości (do ) czas do dokładnego obliczenia prawdopodobieństwa był w sekundach i skalowany kwadratowo, więc zaryzykowałem obliczenie dla do trzech SD powyżej średnia (prawdopodobieństwa 0, 1, ..., 22 sukcesów). Zajęło to 80 minut (z Mathematica 8), zgodnie z przewidywanym czasem. (Wynikowe prawdopodobieństwa to ułamki, których liczniki i mianowniki mają około 75 000 cyfr na sztukę!) To pokazuje, że można wykonać obliczenia. $p_i = 1/(i+1)$ $n$ $n = 2^{16}$ $n$ $n = 2^{12}$ $n = 2^{16}$

Alternatywą jest przeprowadzenie długiej symulacji (należy wykonać milion prób). Należy to zrobić tylko raz, ponieważ się nie zmieniają. $p_i$

— Whuber
źródło

9

(Ponieważ takie podejście jest niezależne od innych opublikowanych rozwiązań, w tym jednego, które opublikowałem, oferuję je jako osobną odpowiedź).

Możesz obliczyć dokładny rozkład w sekundach (lub mniej), pod warunkiem, że suma wartości p jest niewielka.

Widzieliśmy już sugestie, że rozkład może być w przybliżeniu Gaussa (w niektórych scenariuszach) lub Poissona (w innych scenariuszach). W każdym razie wiemy, że jego średnia jest sumą a jej wariancja jest sumą . Dlatego rozkład będzie skoncentrowany w obrębie kilku standardowych odchyleń jego średniej, powiedzmy SD między 4 a 6 lub więcej. Dlatego musimy tylko obliczyć prawdopodobieństwo, że suma jest równa (liczba całkowita) dla do . Kiedy większość $\mu$ $p_i$ $\sigma^2$ $p_i(1-p_i)$ $z$ $z$ $X$ $k$ $k = \mu - z \sigma$ $k = \mu + z \sigma$ $p_i$ są małe, jest w przybliżeniu równa (ale nieco mniejsza niż) , więc aby zachować ostrożność, możemy wykonać obliczenia dla w przedziale . Na przykład, gdy suma wynosi i wybranie aby dobrze pokryć ogony, potrzebowalibyśmy obliczeń, aby pokryć w = , czyli tylko 28 wartości. $\sigma^2$ $\mu$ $k$ $[\mu - z \sqrt{\mu}, \mu + z \sqrt{\mu}]$ $p_i$ $9$ $z = 6$ $k$ $[9 - 6 \sqrt{9}, 9 + 6 \sqrt{9}]$ $[0, 27]$

Rozkład jest obliczany rekurencyjnie . Niech będzie rozkładem sumy pierwszego tych zmiennych Bernoulliego. Dla dowolnego od do suma pierwszych zmiennych może być równa na dwa wzajemnie wykluczające się sposoby: suma pierwszych zmiennych jest równa a wynosi w przeciwnym razie suma pierwszych zmiennych jest równa a wynosi . W związku z tym $f_i$ $i$ $j$ $0$ $i+1$ $i+1$ $j$ $i$ $j$ $i+1^\text{st}$ $0$ $i$ $j-1$ $i+1^\text{st}$ $1$

f_{i + 1} (j) = f_{i} (j) (1 - p_{i + 1}) + f_{i} (j - 1) p_{i + 1} .

$f_{i+1}(j) = f_i(j)(1 - p_{i+1}) + f_i(j-1) p_{i+1}.$

Musimy przeprowadzić to obliczenie dla całki w przedziale od do $j$ $\max(0, \mu - z \sqrt{\mu})$ $\mu + z \sqrt{\mu}.$

Kiedy większość jest niewielkich (ale są nadal odróżnialne od z rozsądną precyzją), podejście to nie jest nękane ogromną kumulacją błędów zaokrągleń zmiennoprzecinkowych używanych w rozwiązaniu, które wcześniej opublikowałem. Dlatego obliczenia o zwiększonej precyzji nie są wymagane. Na przykład obliczenia o podwójnej precyzji dla tablicy prawdopodobieństw ( , wymagające obliczeń dla prawdopodobieństw sum od do $p_i$ $1 - p_i$ $1$ $2^{16}$ $p_i = 1/(i+1)$ $\mu = 10.6676$ $0$ $31$ ) zajęło 0,1 sekundy w przypadku Mathematica 8 i 1-2 sekundy w programie Excel 2002 (oba uzyskały te same odpowiedzi). Powtarzanie go z poczwórną precyzją (w Mathematica) zajęło około 2 sekund, ale nie zmieniło żadnej odpowiedzi więcej niż . Zakończenie rozkładu przy SD w górny ogon straciło tylko całkowitego prawdopodobieństwa. $3 \times 10^{-15}$ $z = 6$ $3.6 \times 10^{-8}$

Kolejne obliczenia dla tablicy 40 000 losowych wartości podwójnej precyzji od 0 do 0,001 ( ) zajęło 0,08 sekundy dla Mathematica. $\mu = 19.9093$

Ten algorytm można zrównoleglać. Po prostu podziel zestaw na rozłączne podzbiory o mniej więcej równej wielkości, po jednym na procesor. Oblicz rozkład dla każdego podzbioru, a następnie zbierz wyniki (używając FFT, jeśli chcesz, chociaż to przyspieszenie prawdopodobnie nie jest konieczne), aby uzyskać pełną odpowiedź. To sprawia, że jest praktyczny w użyciu, nawet gdy staje się duży, gdy trzeba patrzeć daleko w ogony ( duże) i / lub jest duży. $p_i$ $\mu$ $z$ $n$

Czas dla tablicy zmiennych z procesorami jest skalowany jako . Prędkość Mathematiki jest rzędu miliona na sekundę. Na przykład, przy procesorze, zmienia się, całkowite prawdopodobieństwo wynosi , a wychodząc na standardowych odchyleń w górnym ogonie, miliona: oblicz kilka sekund czasu obliczeniowego. Jeśli to skompilujesz, możesz przyspieszyć działanie o dwa rzędy wielkości. $n$ $m$ $O(n(\mu + z \sqrt{\mu})/m)$ $m = 1$ $n = 20000$ $\mu = 100$ $z = 6$ $n(\mu + z \sqrt{\mu})/m = 3.2$

Nawiasem mówiąc, w tych przypadkach testowych wykresy rozkładu wyraźnie pokazały pewną dodatnią skośność: nie są one normalne.

Dla przypomnienia, oto rozwiązanie Mathematica:

pb[p_, z_] := Module[
  {\[Mu] = Total[p]},
  Fold[#1 - #2 Differences[Prepend[#1, 0]] &, 
   Prepend[ConstantArray[0, Ceiling[\[Mu] + Sqrt[\[Mu]] z]], 1], p]
  ]

( Uwaga: kodowanie kolorami stosowane przez tę stronę nie ma znaczenia dla kodu Mathematica. W szczególności szare elementy nie są komentarzami: tam cała praca jest wykonywana!)

Przykładem jego zastosowania jest

pb[RandomReal[{0, 0.001}, 40000], 8]

Edytować

RRozwiązaniem jest dziesięć razy wolniej niż Mathematica w tym przypadku testowego - być może nie zostały zakodowane go optymalnie - ale nadal wykonuje szybko (około jednej sekundy):

pb <- function(p, z) {
  mu <- sum(p)
  x <- c(1, rep(0, ceiling(mu + sqrt(mu) * z)))
  f <- function(v) {x <<- x - v * diff(c(0, x));}
  sapply(p, f); x  
}
y <- pb(runif(40000, 0, 0.001), 8)
plot(y)

Wykres PDF

— Whuber
źródło

8

Z innym najlepszym moim zdaniem jest normalne przybliżeniem. Niech . Następnie $p_i$ $B_n=\sum_{i=1}^np_i(1-p_i)$

a, pod warunkiem, że dla każdego

\begin{aligned} B_{n}^{- 1 / 2} (\sum_{i = 1}^{n} X_{i} - \sum_{i = 1}^{n} p_{i}) \to N (0, 1), \end{aligned}

$\begin{align*} B_n^{-1/2}\left(\sum_{i=1}^nX_i-\sum_{i=1}^np_i\right)\to N(0,1), \end{align*}$

n \to \infty

$n\to\infty$

ε > 0

$\varepsilon>0$

a, co zmiennych Bernoulliego zostanie przytrzymaj, jeśli

\begin{aligned} B_{n}^{- 1} \sum_{i = 1}^{n} E ((X_{i} - p_{i})^{2} 1 {| X_{i} - p_{i} | > ε B_{n}^{1 / 2}}) \to 0, \end{aligned}

$\begin{align*} B_n^{-1}\sum_{i=1}^nE\left((X_i-p_i)^2\mathbf{1}\{|X_i-p_i|>\varepsilon B_n^{1/2}\}\right)\to 0, \end{align*}$

n \to \infty

$n\to\infty$

B_{n} \to \infty

$B_n\to\infty$ . Jest to tak zwany warunek Lindeberga, który jest wystarczający i niezbędny do konwergencji do standardowej normy.

Aktualizacja: błąd aproksymacji można obliczyć na podstawie następujących nierówności:

w którym ijest ED skalowanego ześrodkowany suma.

\begin{aligned} sup_{x} | F_{n} (x) - Φ (x) | \leq A L_{n}, \end{aligned}

$\begin{align*} \sup_x|F_n(x)-\Phi(x)|\le AL_n, \end{align*}$

\begin{aligned} L_{n} = B_{n}^{- 3 / 2} \sum_{i = 1}^{n} E | X_{i} - p_{i} |^{3} \end{aligned}

$\begin{align*} L_n=B_n^{-3/2}\sum_{i=1}^nE|X_i-p_i|^3 \end{align*}$

F_{n}

$F_n$

X_{i}

$X_i$

$p_i$ $p_i=\frac{1}{1+i}$ $B_n\approx \ln n$ $L_n\approx (\ln n)^{-1/2}$ $n=2^{16}$

— mpiktas
źródło

3

Nie jest to prawdą, gdy p_i zbliża się do zera, gdy rośnie. W przeciwnym razie właśnie udowodniono, że rozkład Poissona jest normalny!

— whuber

1

B_{n} \to \infty

$B_n\to\infty$

p_{i}

$p_i$

1 / i

$1/i$

lim B_{n} < \infty

$\lim B_n<\infty$

@mpiktas ma rację. Tutaj analogia do rozkładu Poissona nie do końca pasuje.

Nawiasem mówiąc, tak naprawdę nie sprawdziłem tego monstrualnego stanu w drugim akapicie.

@SOL. Jay Kerns Zgadzam się, że analogia do Poissona jest niedoskonała, ale myślę, że daje dobre wskazówki. Wyobraź sobie sekwencję p, p_i = 10 ^ {- j}, gdzie j jest rzędem wielkości i (równym 1 dla i <= 10, 2 dla i <= 100 itd.). Gdy n = 10 ^ k, 90% p jest równe 10 ^ {- k}, a ich suma wygląda na Poissona z oczekiwaniem 0,9. Kolejne 9% równa się 10 ^ {1-k}, a ich suma wygląda na Poissona (z takimi samymi oczekiwaniami). Zatem rozkład wygląda w przybliżeniu jak suma zmienności k Poissona. To oczywiście nigdzie w pobliżu Normalnego. Stąd potrzeba „monstrualnego stanu”.

— whuber

4

$Y$ $\sum_i p_i$ $\sum_i p_{i}(1-p_{i})$ $Y$ $p_i$ $p_i$ $Y$ $\sum p_i$ $\sum_i p_i$ $y$

$p$ $p$ $p$

$p$ $p$ $Y$

set.seed(1)
M <- 5000
N <- 15000
p <- rbeta(N, shape1 = 1, shape2 = 10)
Y <- replicate(M, sum(rbinom(N, size = 1, prob = p)))

Teraz spójrz na wyniki.

hist(Y)
mean(Y)
sum(p)
var(Y)
sum(p*(1 - p))

Baw się dobrze; Oczywiście zrobiłem.

p

$p$

p

$p$

p

$p$

α < 1

$\alpha < 1$

2

Myślę, że inne odpowiedzi są świetne, ale nie widziałem żadnych bayesowskich metod szacowania prawdopodobieństwa. Odpowiedź nie ma jednoznacznej formy, ale prawdopodobieństwo można zasymulować za pomocą R.

Oto próba:

X_{i} | p_{i} \sim B e r (p_{i})

$X_i | p_i \sim Ber(p_i)$

p_{i} \sim B e t a (α, β)

$p_i \sim Beta(\alpha, \beta)$

$\hat{\alpha}$ $\hat{\beta}$

$i^{th}$ $p_i$ $Beta(\hat{\alpha},\hat{\beta})$ $X_i$ $Ber(p_i)$ $N$ $Y = \sum X_i$ $M$ $M$ Ys będzie oszacowaniem gęstości Y.

P r o b [Y \leq y] = \frac{# Y \leq y}{M}

$Prob[Y \leq y] = \frac {\#Y \leq y} {M}$

$p_i$

— suncoolsu
źródło

1

Dla niektórych purystów może to nie być Bayesian. W rzeczywistości jest to empiryczny Bayesian, ale jest to szybki sposób na symulację prawdopodobieństwa w R, bez uciekania się do hiper-wcześniejszego mumbo jumbo.

— suncoolsu,

1

Dlaczego potrzebujesz priors, gdy podane są p_i?

— whuber

p_{i}

$p_i$

p_{i} \to \frac{α}{α + β}

$p_i\to\frac{\alpha}{\alpha+\beta}$

(\binom{1}{x_{i}}) \frac{B (α + x_{i}, β + 1 - x_{i})}{B (α, β)} = \frac{α^{x_{i}} β^{1 - x_{i}}}{α + β}

${1\choose x_i}\frac{B(\alpha+x_i,\beta+1-x_i)}{B(\alpha,\beta)}=\frac{\alpha^{x_i}\beta^{1-x_i}}{\alpha+\beta}$

p_{i}

$p_i$

p_{1} = p_{2} = \dots = p_{n}

$p_1=p_2=\dots=p_n$

2

Jak wspomniano w innych odpowiedziach, rozkład prawdopodobieństwa, który opisujesz, jest rozkładem dwumianowym Poissona. Wydajną metodę obliczania CDF podano w Hong, Yili. Przy obliczaniu funkcji rozkładu dla rozkładu dwumianowego Poissona .

Podejście polega na wydajnym obliczeniu DFT (dyskretnej transformaty Fouriera) funkcji charakterystycznej.

$\phi(t) = \prod_j^n [(1-p_j)+p_je^{it}]$ $i=\sqrt{-1}$

Algorytm to:

$z_j(k) = 1-p_j+p_j \text{cos}(\omega k)+ i p_j \text{sin}(\omega k)$ $\omega=\frac{2\pi}{n+1}$
$x_k=\text{exp}\{\sum_j^n log(z_j(k))\}$ $x_0=1$
$x_k$ $k=1,\dots,[n/2]$ $\bar{x}_k=x_{n+1-k}$
$\frac{1}{n+1}<x_0,x_1,\dots,x_n>$
Weź łączną sumę wyniku, aby uzyskać CDF.

Algorytm jest dostępny w pakiecie poibin R.

Takie podejście daje znacznie lepsze wyniki niż preparaty rekurencyjne, ponieważ mają tendencję do braku stabilności numerycznej.

— Kyle
źródło

3

Mam dostęp tylko do streszczenia tego artykułu, ale wygląda na to, że implementuje metodę, której użyłem na stronie stats.stackexchange.com/questions/41247/... i omawia, jak działa w porównaniu z innymi metodami podanymi w tym wątku. Jeśli wiesz więcej o tym, co udało się osiągnąć w pracy, chętnie przeczytamy streszczenie.

— whuber

1

$Y$ $Z$ $\sum_i p_i$

sup_{A} | P (Y \in A) - P (Z \in A) | \leq min {1, \frac{1}{\sum_{i} p_{i}}} \sum_{i} p_{i}^{2} .

$\sup_A |{\bf P}(Y\in A) - {\bf P}(Z\in A)| \le \min \left\{ 1, \frac{1}{\sum_i p_i} \right\} \sum_i p_i^2.$

$|{\bf E}f(Y) - {\bf E}f(Z)|$ $f$ $A$

P (Y \in A) \leq \frac{1}{(1 - max_{i} p_{i})^{2}} P (Z \in A) .

${\bf P}(Y\in A) \le \frac{1}{(1-\max_i p_i)^2} {\bf P}(Z\in A).$

— Pavel Ruzankin
źródło

1

+1 Dziękujemy za przydatne informacje ilościowe na temat granic aproksymacji. Witamy na naszej stronie!

— whuber