Rozkład prawdopodobieństwa dla różnych prawdopodobieństw

Gdybym chciał uzyskać prawdopodobieństwo 9 sukcesów w 16 próbach z prawdopodobieństwem 0,6 każdej próby, mógłbym zastosować rozkład dwumianowy. Czego mogę użyć, jeśli każda z 16 prób ma inne prawdopodobieństwo sukcesu?

Jest to suma 16 (prawdopodobnie niezależnych) prób dwumianowych. Założenie niezależności pozwala nam pomnożyć prawdopodobieństwa. Stąd po dwóch próbach z prawdopodobieństwem sukcesu i szansa na sukces w obu próbach wynosi , szansa na brak sukcesu wynosi , a szansa na jeden sukces wynosi . To ostatnie wyrażenie jest ważne ze względu na fakt, że dwa sposoby osiągnięcia dokładnie jednego sukcesu wykluczają się wzajemnie: co najwyżej jeden może się zdarzyć. Oznacza to, że dodają ich prawdopodobieństwa .p 2 p 1 p 2 ( 1 - p 1 ) ( 1 - p 2 ) p 1 ( 1 - p 2 ) + ( 1 - p 1 ) p $p_1$$p_2$$p_1 p_2$$(1-p_1)(1-p_2)$$p_1(1-p_2) + (1-p_1)p_2$

Za pomocą tych dwóch zasad - mnożą się niezależne prawdopodobieństwa i wzajemnie się wykluczają - możesz wypracować odpowiedzi dla, powiedzmy, 16 prób z prawdopodobieństwami . Aby to zrobić, musisz wziąć pod uwagę wszystkie sposoby uzyskania każdej podanej liczby sukcesów (np. 9). Istnieje sposobów na osiągnięcie 9 sukcesów. Jeden z nich występuje na przykład, gdy próby 1, 2, 4, 5, 6, 11, 12, 14 i 15 są sukcesami, a inne porażkami. Sukcesy miały prawdopodobieństwa i a awarie miały prawdopodobieństwa . Pomnożenie tych 16 liczb daje szansę$p_1, \ldots, p_{16}$$\binom{16}{9} = 11440$p 15 1 - p 3 , 1 - p 7 , … , 1 - p 13 , 1 - p $p_1, p_2, p_4, p_5, p_6, p_{11}, p_{12}, p_{14},$$p_{15}$$1-p_3, 1-p_7, \ldots, 1-p_{13}, 1-p_{16}$tej konkretnej sekwencji wyników. Sumowanie tej liczby wraz z 11 439 pozostałymi takimi liczbami daje odpowiedź.

Oczywiście użyłbyś komputera.

Przy wielu ponad 16 próbach konieczne jest przybliżenie dystrybucji. Pod warunkiem, że żadne z prawdopodobieństw i stają się zbyt małe, normalne przybliżenie zwykle działa dobrze. Za pomocą tej metody można zauważyć, że oczekiwanie na sumę prób wynosi i (ponieważ próby są niezależne) wariancja wynosi . Udajesz, że rozkład sum jest normalny ze średnią i odchyleniem standardowym . Odpowiedzi wydają się być dobre do obliczania prawdopodobieństw odpowiadających odsetkowi sukcesów, który różni się 1 - p i n μ = p 1 + p 2 + ⋯ + p n σ 2 = p 1 ( 1 - p 1 ) + p 2 ( 1 - p 2 ) + ⋯ + p n ( 1 - p n ) μ σ μ σ n σ $p_i$$1-p_i$$n$$\mu = p_1 + p_2 + \cdots + p_n$$\sigma^2 = p_1(1-p_1) + p_2(1-p_2) + \cdots + p_n(1-p_n)$$\mu$$\sigma$$\mu$ przez nie więcej niż kilka wielokrotności . Gdy rośnie, to przybliżenie staje się coraz bardziej dokładne i działa dla jeszcze większych wielokrotności z dala od .$\sigma$$n$$\sigma$$\mu$

Jedną alternatywą dla normalnego przybliżenia @ whubera jest użycie prawdopodobieństwa „mieszania” lub modelu hierarchicznego. Miałoby to zastosowanie, gdy są w pewien sposób podobne i można to modelować za pomocą rozkładu prawdopodobieństwa z funkcją gęstości indeksowaną przez jakiś parametr . otrzymujesz równanie całkowe:p i ∼ D i s t ( θ ) g ( p | θ ) $p_i$$p_i\sim Dist(\theta)$$g(p|\theta)$$\theta$

$$Pr(s=9|n=16,\theta)={16 \choose 9}\int_{0}^{1} p^{9}(1-p)^{7}g(p|\theta)dp$$

Prawdopodobieństwo dwumianowe pochodzi z ustawienia , normalne przybliżenie pochodzi z (myślę) ustawienia (z i zgodnie z definicją w odpowiedzi @ whubera), a następnie odnotowując „ ogony ”tego pliku PDF gwałtownie opadają wokół szczytu.g ( p | θ ) = g ( p | μ , σ ) = $g(p|\theta)=\delta(p-\theta)$μ$g(p|\theta)=g(p|\mu,\sigma)=\frac{1}{\sigma}\phi\left(\frac{p-\mu}{\sigma}\right)$$\mu$$\sigma$

Można również użyć rozkładu wersji beta, który prowadziłby do prostej formy analitycznej i który nie musi cierpieć z powodu problemu „małego p”, jaki wykonuje normalne przybliżenie - ponieważ beta jest dość elastyczny. Przy użyciu rozkładu z ustawionymi przez rozwiązania następujących równań (jest to oszacowanie „minimalnej dywergencji KL”):α , $beta(\alpha,\beta)$$\alpha,\beta$

ψ(β)-ψ(α+β)=

$$\psi(\alpha)-\psi(\alpha+\beta)=\frac{1}{n}\sum_{i=1}^{n}log[p_{i}]$$ $$\psi(\beta)-\psi(\alpha+\beta)=\frac{1}{n}\sum_{i=1}^{n}log[1-p_{i}]$$

Gdzie To funkcja digamma - ściśle związana z szeregami harmonicznymi.$\psi(.)$

Otrzymujemy rozkład związku „dwumianowego”:

$${16 \choose 9}\frac{1}{B(\alpha,\beta)}\int_{0}^{1} p^{9+\alpha-1}(1-p)^{7+\beta-1}dp ={16 \choose 9}\frac{B(\alpha+9,\beta+7)}{B(\alpha,\beta)}$$

Rozkład ten jest zbliżony do rozkładu normalnego w przypadku, gdy @whuber wskazuje - ale powinien dać rozsądne odpowiedzi dla małego i pochylonego - ale nie dla multimodalnego , ponieważ rozkład beta ma tylko jeden pik. Ale możesz to łatwo naprawić, po prostu używając dystrybucji beta dla trybówDzielimy całkę z na części, aby każdy element miał unikalny tryb (i wystarczającą ilość danych do oszacowania parametrów) i dopasował rozkład beta w każdym kawałku. następnie zsumuj wyniki, zauważając, że dokonanie zmiany zmiennych dlap i p i M M 0 < p < 1 M p = x - $n$$p_i$$p_i$$M$$M$$0<p<1$$M$ L<x<$p=\frac{x-L}{U-L}$$L<x<U$ całka beta przekształca się w:

$$B(\alpha,\beta)=\int_{L}^{U}\frac{(x-L)^{\alpha-1}(U-x)^{\beta-1}}{(U-L)^{\alpha+\beta-1}}dx$$

Niech ~ z funkcją generującą prawdopodobieństwo (pgf): B e r n o u l l i ( p i$X_i$$Bernoulli(p_i)$

$$\text{pgf} = E[t^{X_i}] = 1 - p_i (1-t)$$

Niech oznacza sumę takich niezależnych zmiennych losowych. Następnie PGF dla sumy z takich zmiennych jest: n S n = $S = \sum_{i=1}^n X_i$$n$$S$$n=16$

$$\begin{align*}\displaystyle \text{pgfS} &= E[t^S] \\&= E[t^{X_1}] E[t^{X_2}] \dots E[t^{X_{16}}] \text{ (... by independence)} \\ &= \prod _{i=1}^{16} \left(1-p_i(1-t) \right)\end{align*}$$

Szukamy , czyli:$P(S=9)$

$$\frac{1}{9!}\frac{d^9 \text{pgfS}}{dt^9}|_{t=0}$$

WSZYSTKIE GOTOWE. Daje to dokładne rozwiązanie symboliczne jako funkcję . Odpowiedź jest dość długa, aby wydrukować ją na ekranie, ale jest całkowicie wykonalna i zajmuje mniej niż sekundy do oceny za pomocą Mathematiki na moim komputerze.$p_i$$\frac{1}{100}$

Przykłady

Jeśli , to: P(S=9)=$p_i = \frac{i}{17}, i= 1 \text{ to } 16$$P(S=9) = \frac{9647941854334808184}{48661191875666868481} = 0.198268 \dots$

Jeśli , to: P(S=9)=0,000228613$p_i = \frac{\sqrt{i}}{17}, i= 1 \text{ to } 16$$P(S=9) = 0.000228613 \dots$

Ponad 16 prób?

Przy ponad 16 próbach nie ma potrzeby przybliżania dystrybucji. Powyższa dokładna metoda działa równie łatwo dla przykładów z powiedzmy lub . Na przykład, gdy , zajmuje mniej niż sekundy do oceny całego pmf ( tj. Przy każdej wartości ) przy użyciu poniższego kodu.n = 100 n = 50 $n = 50$$n = 100$$n = 50$ s=0,1,…,$\frac{1}{10}$$s = 0, 1, \dots, 50$

Kod matematyczny

Biorąc pod uwagę wektor wartości , powiedz:$p_i$

n = 16;   pvals = Table[Subscript[p, i] -> i/(n+1), {i, n}];

... oto kod Mathematica, aby zrobić wszystko, co konieczne:

pgfS = Expand[ Product[1-(1-t)Subscript[p,i], {i, n}] /. pvals];
D[pgfS, {t, 9}]/9! /. t -> 0  // N

0,198268

Aby uzyskać cały pmf:

Table[D[pgfS, {t,s}]/s! /. t -> 0 // N, {s, 0, n}]

... lub użyj jeszcze ładniejszego i szybszego (dzięki sugestii Raya Koopmana poniżej):

CoefficientList[pgfS, t] // N

Na przykład przy , obliczenie zajmuje tylko 1 sekundę , a następnie 0,002 sekundy, aby uzyskać cały pmf przy użyciu , więc jest to niezwykle wydajne.$n = 1000$pgfSCoefficientList

Komentarz @wolfies, a moja próba odpowiedzi na to pytanie ujawniła ważny problem z moją drugą odpowiedzią, o której omówię później.

Przypadek szczególny (n = 16)

Istnieje dość skuteczny sposób na zakodowanie pełnej dystrybucji przy użyciu „sztuczki” polegającej na użyciu liczb podstawowych (binarnych) w obliczeniach. Wymaga tylko 4 wierszy kodu R, aby uzyskać pełny rozkład gdzie . Zasadniczo istnieją w sumie wyborów wektora które mogą przyjąć zmienne binarne . Teraz załóżmy, że każdy z osobnych wyborów numerujemy od do . To samo w sobie nie jest niczym specjalnym, ale teraz załóżmy, że reprezentujemy „liczbę wyboru” za pomocą arytmetyki podstawy 2. Teraz weź , abym mógł zapisać wszystkie opcje, aby były P r ( Z i = 1 ) = p i 2 n z = ( z 1 , … , z n ) Z i 1 2 n n = 3 2 3 = 8 1 , 2 , 3 , 4 , 5 , 6 , 7 , 8 $Y=\sum_{i=1}^{n} Z_i$$Pr(Z_i=1)=p_i$$2^n$$z=(z_1,\dots,z_n)$$Z_i$$1$$2^n$$n=3$$2^3=8$wybory Następnie w „liczbach zwykłych” staje się w „liczbach binarnych”. Załóżmy teraz, że zapisujemy je jako liczby czterocyfrowe, a następnie mamy . Teraz spójrz na ostatnie cyfry każdej liczby - można traktować jako itd. Liczenie w postaci binarnej zapewnia skuteczny sposób organizacji sumowania . Na szczęście istnieje funkcja R, która może wykonać dla nas tę konwersję binarną, wywołana i przekonwertujemy surową formę binarną na liczbę za pomocą , wtedy otrzymamy wektor z$1,2,3,4,5,6,7,8$0001 , 0010 , 0011 , 0100 , 0101 , 0110 , 0111 , 1000 3 001 ( Z 1 = 0 , Z 2 = 0 , Z 3 = 1 $1,10,11,100,101,110,111,1000$$0001,0010,0011,0100,0101,0110,0111,1000$$3$$001$32 y = $(Z_1=0,Z_2=0,Z_3=1)\implies Y=1$intToBits(x)as.numeric(intToBits(x))$32$elementy, przy czym każdy element jest cyfrą podstawowej wersji 2 naszego numeru (czytany od prawej do lewej, nie od lewej do prawej). Używając tej sztuczki w połączeniu z innymi wektoryzacjami R, możemy obliczyć prawdopodobieństwo, że w 4 liniach kodu R:$y=9$

exact_calc <- function(y,p){
    n       <- length(p)
    z       <- t(matrix(as.numeric(intToBits(1:2^n)),ncol=2^n))[,1:n] #don't need columns n+1,...,32 as these are always 0
    pz      <- z%*%log(p/(1-p))+sum(log(1-p))
    ydist   <- rowsum(exp(pz),rowSums(z))
    return(ydist[y+1])
}

Podłączenie jednolitego przypadku i przypadek root sqrt daje pełną dystrybucję dla y jako: p ( 2 ) i =$p_i^{(1)}=\frac{i}{17}$$p_i^{(2)}=\frac{\sqrt{i}}{17}$

$$\begin{array}{c|c}y & Pr(Y=y|p_i=\frac{i}{17}) & Pr(Y=y|p_i=\frac{\sqrt{i}}{17})\\ \hline 0 & 0.0000 & 0.0558 \\ 1 & 0.0000 & 0.1784 \\ 2 & 0.0003 & 0.2652 \\ 3 & 0.0026 & 0.2430 \\ 4 & 0.0139 & 0.1536 \\ 5 & 0.0491 & 0.0710 \\ 6 & 0.1181 & 0.0248 \\ 7 & 0.1983 & 0.0067 \\ 8 & 0.2353 & 0.0014 \\ 9 & 0.1983 & 0.0002 \\ 10 & 0.1181 & 0.0000 \\ 11 & 0.0491 & 0.0000 \\ 12 & 0.0139 & 0.0000 \\ 13 & 0.0026 & 0.0000 \\ 14 & 0.0003 & 0.0000 \\ 15 & 0.0000 & 0.0000 \\ 16 & 0.0000 & 0.0000 \\ \end{array}$$

Tak dla specyficznego problemu z sukcesów w próbach, dokładne obliczenia są prosto do przodu. Działa to również dla wielu prawdopodobieństw do około - poza tym prawdopodobnie zaczniesz mieć problemy z pamięcią i potrzebne są różne sztuczki obliczeniowe.16 n = $y$$16$$n=20$

Zauważ, że stosując mój sugerowany „rozkład beta” otrzymujemy oszacowania parametrów a to daje oszacowanie prawdopodobieństwa, które jest prawie równomierne w , dając przybliżoną wartość . Wydaje się to dziwne, biorąc pod uwagę, że gęstość rozkładu beta z ściśle zbliża się do histogramu wartości . Co poszło nie tak?y p r ( y = 9 ) = 0,06799 ≈ $\alpha=\beta=1.3206$$y$ α=β=1,3206p$pr(y=9)=0.06799\approx\frac{1}{17}$$\alpha=\beta=1.3206$$p_i$

Ogólna sprawa

Omówię teraz bardziej ogólny przypadek i dlaczego moje proste przybliżenie wersji beta nie powiodło się. Zasadniczo, pisząc a następnie mieszając nad z innym rozkładem faktycznie przyjmuje ważne założenie - że możemy oszacować rzeczywiste prawdopodobieństwo za pomocą jedno prawdopodobieństwo dwumianowe - jedynym problemem, który pozostaje, jest to, której wartości użyć. Jednym ze sposobów na to jest użycie gęstości mieszania, która jest dyskretnie jednorodna w stosunku do rzeczywistego . Zastępujemy więc rozkład beta gęstością dyskretnąp p ∼ f ( θ ) p p i p ∼ B e t a ( a , b ) p ∼ ∑ 16 i = 1 w i δ ( p - p i ) p i w i p $(y|n,p)\sim Binom(n,p)$$p$$p\sim f(\theta)$$p$$p_i$$p\sim Beta(a,b)$$p\sim \sum_{i=1}^{16}w_i\delta(p-p_i)$. Następnie za pomocą aproksymacji mieszania można wyrazić słowami, wybierając wartość z prawdopodobieństwem i zakładając, że wszystkie próby mają to prawdopodobieństwo$p_i$$w_i$ . Oczywiście, aby takie przybliżenie działało dobrze, większość wartości powinna być do siebie podobna. Zasadniczo oznacza to, że dla równomiernego rozkładu wartości @wolfies, powoduje bardzo złe przybliżenie przy użyciu rozkładu mieszania beta. To wyjaśnia również, dlaczego przybliżenie jest znacznie lepsze dla - są one mniej rozłożone.$p_i$ pi=$p_i=\frac{i}{17}$$p_i=\frac{\sqrt{i}}{17}$

Mieszanie wykorzystuje następnie zaobserwowane do uśrednienia wszystkich możliwych wyborów pojedynczego . Ponieważ „miksowanie” jest jak średnia ważona, nie może być nic lepszego niż użycie pojedynczego najlepszego . Jeśli więc są wystarczająco rozłożone, nie może istnieć pojedynczy który mógłby zapewnić dobre przybliżenie wszystkich . p p p i p p $p_i$ $p$$p$$p_i$$p$$p_i$

Jedną rzeczą, którą powiedziałem w mojej drugiej odpowiedzi było to, że może być lepiej użyć mieszanki rozkładów beta w ograniczonym zakresie - ale to nadal nie pomoże tutaj, ponieważ wciąż miesza się w jednym . Bardziej sensowny jest podział przedziału na kawałki i dwumian w każdym kawałku. Na przykład możemy wybrać jako nasze podziały i dopasować dziewięć dwumianów w każdym zakresie prawdopodobieństwa . Zasadniczo w ramach każdego podziału pasowalibyśmy do prostego przybliżenia, takiego jak użycie dwumianu z prawdopodobieństwem równym średniej( 0 , 1 ) ( 0 , 0,1 , 0,2 , … , 0,9 , 1 ) 0,1 p $p$$(0,1)$$(0,0.1,0.2,\dots,0.9,1)$$0.1$$p_i$w tym zakresie. Jeśli zmniejszymy odstępy, przybliżenie staje się arbitralnie dobre. Ale zauważ, że wszystko to sprawia, że ​​mamy do czynienia z sumą niezależnych prób dwumianowych z różnymi prawdopodobieństwami, zamiast prób Bernoulliego . Jednak poprzednia część tej odpowiedzi pokazała, że ​​możemy wykonać dokładne obliczenia pod warunkiem, że liczba dwumianów jest wystarczająco mała, powiedzmy 10-15 lub więcej.

Aby rozszerzyć odpowiedź opartą na bernoulli na odpowiedź dwumianową, po prostu „ponownie interpretujemy” zmienne . Po prostu stwierdzamy, że - redukuje to do pierwotnego opartego na ale teraz mówi, z których dwumianów odnoszą sukcesy. Tak więc przypadek oznacza teraz, że wszystkie „sukcesy” pochodzą z trzeciego dwumianu, a żaden z pierwszych dwóch.Z i = I ( X i > 0 ) Z i ( Z 1 = 0 , Z 2 = 0 , Z 3 = 1 $Z_i$$Z_i=I(X_i>0)$$Z_i$$(Z_1=0,Z_2=0,Z_3=1)$

Zauważ, że wciąż jest to „wykładnicze”, ponieważ liczba obliczeń jest podobna do gdzie jest liczbą dwumianów, a jest wielkością grupy - więc masz gdzie . Jest to jednak lepsze niż , z którymi miałbyś do czynienia przy użyciu losowych zmiennych bernoulli. Załóżmy na przykład, że podzieliliśmy prawdopodobieństw na grupy z prawdopodobieństwami w każdej grupie. Daje to obliczeń, w porównaniu do$k^g$$g$$k$$Y\approx\sum_{j=1}^{g}X_j$$X_j\sim Bin(k,p_j)$$2^{gk}$$n=16$$g=4$$k=4$$4^4=256$$2^{16}=65536$

Wybierając grup i zauważając, że limit wynosił około czyli około komórek, możemy skutecznie zastosować tę metodę, aby zwiększyć maksimum do .$g=10$$n=20$$10^7$$n$$n=50$

Jeśli dokonamy przybliżenia crudera, obniżając , zwiększymy „możliwy” rozmiar dla . oznacza, że można mieć efektywną około . Poza tym normalne przybliżenie powinno być niezwykle dokładne.$g$$n$$g=5$$n$$125$

(Na ogół trudny do wyodrębnienia) pmf to Kod R:

$$\Pr(S=k) = \sum_{\substack{A\subset\{1,\dots,n\}\\ |A|=k}} \left( \prod_{i\in A} p_i \right)\left(\prod_{j\in \{1,\dots,n\}\setminus A} (1-p_j) \right) \, .$$

p <- seq(1, 16) / 17
cat(p, "\n")
n <- length(p)
k <- 9
S <- seq(1, n)
A <- combn(S, k)
pr <- 0
for (i in 1:choose(n, k)) {
    pr <- pr + exp(sum(log(p[A[,i]])) + sum(log(1 - p[setdiff(S, A[,i])])))
}
cat("Pr(S = ", k, ") = ", pr, "\n", sep = "")

Dla używanych w odpowiedzi na wilki mamy:$p_i$

Pr(S = 9) = 0.1982677

Kiedy rośnie, użyj splotu .$n$