Interwał przewidywania dla losowej zmiennej dwumianowej

Jaka jest formuła (przybliżona lub dokładna) przedziału predykcji dla losowej zmiennej dwumianowej?

Załóżmy, że , i obserwujemy y (na podstawie Y ). N jest znana.$Y \sim \mathsf{Binom}(n, p)$$y$$Y$$n$

Naszym celem jest uzyskanie 95% przedział predykcji dla nowego czerpać z .$Y$

Estymacja punktowa jest , gdzie p = $n\hat{p}$ . Przedział ufności dla p jest proste, ale nie mogę znaleźć formułę dla przedziału predykcji dlaY. Gdybyśmy wiedzieli,p(zamiast p ), a następnie 95% przedział predykcji właśnie polega znalezieniu quantiles o dwumianowy. Czy coś oczywistego przeoczam?$\hat{p}=\frac{y}{n}$$\hat{p}$$Y$$p$$\hat{p}$

Ok, spróbujmy tego. Dam dwie odpowiedzi - bayesowską, która moim zdaniem jest prosta i naturalna, i jedną z możliwych częstych.

Rozwiązanie bayesowskie

Zakładamy beta przed na , I, np., P ~ B E T a ( alfa , beta ) , ponieważ model beta dwumianowego jest sprzężone, co oznacza, że rozkład tylny jest również beta dystrybucyjnym parametry α = α + K , β = β + n - K (ja pomocą k w celu określenia liczby sukcesów n badaniach zamiast y ). Zatem wnioskowanie jest znacznie uproszczone. Teraz, jeśli masz wcześniejszą wiedzę na temat prawdopodobnych wartości$p$$p \sim Beta(\alpha,\beta)$$\hat{\alpha}=\alpha+k,\hat{\beta}=\beta+n-k$$k$$n$$y$ , można go użyć do ustawienia wartości α i β , tj. do zdefiniowania wcześniejszej Beta, w przeciwnym razie można założyć jednolity (nieinformacyjny) wcześniej, z α = β = 1 lub innymi nieinformacyjnymi priorytetami (patrz na przykładtutaj). W każdym razie twój tylny jest$p$$\alpha$$\beta$$\alpha=\beta=1$

$Pr(p|n,k)=Beta(\alpha+k,\beta+n-k)$

W wnioskowaniu bayesowskim liczy się tylko prawdopodobieństwo późniejsze, co oznacza, że ​​kiedy się o tym dowiesz, możesz wyciągać wnioski dla wszystkich innych wielkości w swoim modelu. Chcesz wnioskować na podstawie obserwowalnych : w szczególności na wektorze nowych wyników y = y 1 , … , y m , gdzie m niekoniecznie jest równe n . W szczególności dla każdego j = 0 , … , m , chcemy obliczyć prawdopodobieństwo osiągnięcia dokładnie j sukcesów w następnych m próbach, biorąc pod uwagę, że otrzymaliśmy $y$$\mathbf{y}=y_1,\dots,y_m$$m$$n$$j=0,\dots,m$$j$$m$$k$sukcesy w poprzednich próbach; tylna predykcyjna funkcja masy:$n$

$Pr(j|m,y)=Pr(j|m,n,k)=\int_0^1 Pr(j,p|m,n,k)dp = \int_0^1 Pr(j|p,m,n,k)Pr(p|n,k)dp$

Jednak nasz dwumianowy model oznacza, że warunkowo na str mający pewną wartość, prawdopodobieństwo konieczności j sukcesów w m prób nie zależy od ostatnich wyników: to po prostu$Y$$p$$j$$m$

$f(j|m,p)=\binom{j}{m} p^j(1-p)^j$

W ten sposób wyrażenie staje się

$Pr(j|m,n,k)=\int_0^1 \binom{j}{m} p^j(1-p)^j Pr(p|n,k)dp=\int_0^1 \binom{j}{m} p^j(1-p)^j Beta(\alpha+k,\beta+n-k)dp$

Wynikiem tej całki jest dobrze znany rozkład zwany rozkładem dwumianowym: pomijając fragmenty, otrzymujemy okropny wyraz

$Pr(j|m,n,k)=\frac{m!}{j!(m-j)!}\frac{\Gamma(\alpha+\beta+n)}{\Gamma(\alpha+k)\Gamma(\beta+n-k)}\frac{\Gamma(\alpha+k+j)\Gamma(\beta+n+m-k-j)}{\Gamma(\alpha+\beta+n+m)}$

Nasz punkt oszacowania dla , przy uwzględnieniu straty kwadratowej, jest oczywiście średnią tego rozkładu, tj.$j$

$\mu=\frac{m(\alpha+k)}{(\alpha+\beta+n)}$

Teraz spójrzmy na przedział przewidywania. Ponieważ jest to rozkład dyskretny, nie mamy wyrażenia w postaci zamkniętej dla , takiego, że P r ( j 1 ≤ j ≤ j 2 ) = 0,95 . Powodem jest to, że w zależności od tego, jak zdefiniujesz kwantyl, dla dyskretnego rozkładu funkcja kwantylu albo nie jest funkcją, albo jest funkcją nieciągłą. Ale to nie jest duży problem: dla małego m można po prostu zapisać m prawdopodobieństwa P r ( j = $[j_1,j_2]$$Pr(j_1\leq j \leq j_2)= 0.95$$m$$m$ i stąd znajdź j 1 , j 2 takie, że$Pr(j=0|m,n,k),Pr(j\leq 1|m,n,k),\dots,Pr(j \leq m-1|m,n,k)$$j_1,j_2$

$Pr(j_1\leq j \leq j_2)=Pr(j\leq j_2|m,n,k)-Pr(j < j_1|m,n,k)\geq 0.95$

Oczywiście można znaleźć więcej niż jedną parę, więc idealnie byłoby szukać najmniejszej takiej, aby powyższe było spełnione. Zauważ, że$[j_1,j_2]$

$Pr(j=0|m,n,k)=p_0,Pr(j\leq 1|m,n,k)=p_1,\dots,Pr(j \leq m-1|m,n,k)=p_{m-1}$

są tylko wartościami CMF (Cumulative Mass Function) rozkładu Beta-Binomial, i jako taki istnieje wyrażenie postaci zamkniętej , ale jest to pod względem uogólnionej funkcji hipergeometrycznej, a zatem jest dość skomplikowane. Wolałbym po prostu zainstalować pakiet R extraDistri wywołać, pbbinomaby obliczyć CMF dystrybucji Beta-Binomial. W szczególności, jeśli chcesz obliczyć wszystkie prawdopodobieństwa za jednym razem, po prostu napisz:$p_0,\dots,p_{m-1}$

library(extraDistr)  
jvec <- seq(0, m-1, by = 1) 
probs <- pbbinom(jvec, m, alpha = alpha + k, beta = beta + n - k)

gdzie alphai betasą wartościami parametrów Beta przed, tj. i β (a więc 1, jeśli używasz munduru przed p ). Oczywiście wszystko byłoby znacznie prostsze, gdyby R zapewniał funkcję kwantylową dla rozkładu Beta-Dwumianowego, ale niestety nie.$\alpha$$\beta$$p$

Praktyczny przykład z rozwiązaniem bayesowskim

Niech , k = 70 (dlatego początkowo zaobserwowaliśmy 70 sukcesów w 100 próbach). Chcemy oszacowania punktowego i przedziału 95% prognozy dla liczby sukcesów jw następnych m = 20 próbach. Następnie$n=100$$k=70$$j$$m=20$

n <- 100
k <- 70
m <- 20
alpha <- 1
beta  <- 1

$p$

bayesian_point_estimate <- m * (alpha + k)/(alpha + beta + n) #13.92157

$j$

jvec <- seq(0, m-1, by = 1)
library(extraDistr)
probabilities <- pbbinom(jvec, m, alpha = alpha + k, beta = beta + n - k)

Prawdopodobieństwa są

> probabilities
 [1] 1.335244e-09 3.925617e-08 5.686014e-07 5.398876e-06
 [5] 3.772061e-05 2.063557e-04 9.183707e-04 3.410423e-03
 [9] 1.075618e-02 2.917888e-02 6.872028e-02 1.415124e-01
[13] 2.563000e-01 4.105894e-01 5.857286e-01 7.511380e-01
[17] 8.781487e-01 9.546188e-01 9.886056e-01 9.985556e-01

$j_2$$Pr(j\leq j_2|m,n,k)\ge 0.975$$j_1$$Pr(j < j_1|m,n,k)=Pr(j \le j_1-1|m,n,k)\le 0.025$

$Pr(j_1\leq j \leq j_2|m,n,k)=Pr(j\leq j_2|m,n,k)-Pr(j < j_1|m,n,k)\ge 0.975-0.025=0.95$

$j_2=18$$j_1=9$$Pr(j_1\leq j \leq j_2|m,n,k)\ge 0.95$

Rozwiązanie dla częstych

$Y\sim Binom(m,p)$$X\sim Binom(n,p)$$1-2\alpha-$$Y$$X$$I=[L(X;n,m,\alpha),U(X;n,m,\alpha)]$

$Pr_{X,Y}(Y\in I)=Pr_{X,Y}(L(X;n,m,\alpha)\leq Y\leq U(X;n,m,\alpha)]\geq 1-2\alpha$

$\geq 1-2\alpha$$X$$X+Y=k+j=s$$s$$n$$n+m$

$Pr(X=k|X+Y=s,n,n+m)=\frac{\binom{n}{k}\binom{m}{s-k}}{\binom{m+n}{s}}$

$X$$X+Y=s$

$Pr(X\leq k|s,n,n+m)=H(k;s,n,n+m)=\sum_{i=0}^k\frac{\binom{n}{i}\binom{m}{s-i}}{\binom{m+n}{s}}$

$p$$k$$1-\alpha$$L$

$Pr(X\geq k|k+L,n,n+m)=1-H(k-1;k+L,n,n+m)>\alpha$

$1-\alpha$

$Pr(X\leq k|k+U,n,n+m)=H(k;k+U,n,n+m)>\alpha$

$[L,U]$$Y$$1-2\alpha$$p$$n$$m$$1-2\alpha$

Praktyczny przykład z rozwiązaniem Frequentist

$\alpha$$\beta$

n <- 100
k <- 70
m <- 20

$\hat{p}=\frac{k}{n}$$m$

frequentist_point_estimate <- m * k/n #14

W przypadku przedziału przewidywania procedura jest nieco inna. Szukamy największego takiego, że P r ( X $U$$Pr(X\leq k|k+U,n,n+m)=H(k;k+U,n,n+m)>\alpha$$U$$[0,m]$

jvec <- seq(0, m, by = 1)
probabilities <- phyper(k,n,m,k+jvec)

$U$

jvec[which.min(probabilities > 0.025) - 1] # 18

To samo, co w przypadku podejścia bayesowskiego. Dolna granica przewidywania jest najmniejszą liczbą całkowitą taką, że$L$$Pr(X\geq k|k+L,n,n+m)=1-H(k-1;k+L,n,n+m)>\alpha$

probabilities <- 1-phyper(k-1,n,m,k+jvec)
jvec[which.max(probabilities > 0.025) - 1] # 8

$[L,U]=[8,18]$