Interwał przewidywania dla losowej zmiennej dwumianowej

14

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.YBinom(n,p)yYn

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

Estymacja punktowa jest , gdzie p = ynp^ . 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?p^=ynp^Ypp^

Statseeker
źródło
1
Zobacz Jakie nie-bayesowskie metody wnioskowania predykcyjnego? . W tym przypadku metoda wykorzystująca elementy przestawne nie jest dostępna (nie sądzę), ale można użyć jednego z prawdopodobieństw predykcyjnych. Lub oczywiście podejście bayesowskie.
Scortchi - Przywróć Monikę
1
Cześć chłopaki, chciałbym poświęcić chwilę, aby rozwiązać zgłoszone obawy. - odnośnie zaufania dla p: nie jestem tym zainteresowany. - jeśli chodzi o przewidywania stanowiące 95% rozkładu: tak, dokładnie takie są przedziały prognozowania niezależnie od kontekstu (w regresji musisz założyć normalne błędy, gdzie jako przedziały ufności opierają się na CLT - tak, przykład przewidywania liczby głów w rzut monetą jest prawidłowy. Co sprawia, że ​​ten problem jest trudny, to że nie „p”, jutro mamy oszacowanie.
Statseeker
3
@Addison Przeczytaj książkę Interwały statystyczne G. Hahna i W. Meekera. Wyjaśniają różnicę między przedziałami ufności, przedziałami prognoz, przedziałami tolerancji i przedziałami wiarygodności Bayesa. 95% przedział predykcji nie zawiera 95% rozkładu. Robi to, co robią najczęściej interwały. Jeśli wielokrotnie próbkujesz z B (n, p) i używasz tej samej metody za każdym razem, aby uzyskać 95% przedział predykcji dla p, wówczas 95% przedziałów predykcji będzie zawierać prawdziwą wartość p. Jeśli chcesz pokryć 95% rozkładu, zbuduj przedział tolerancji.
Michael R. Chernick,
Przedziały tolerancji obejmują procent rozkładu. W przypadku 95% przedziału tolerancji dla 90% rozkładu powtórzysz proces wiele razy i używasz tej samej metody do generowania przedziału za każdym razem, a następnie w około 95% przypadków co najmniej 90% rozkładu przypada na przedział i 5% czasu mniej niż 90% rozkładu będzie zawarte w przedziale.
Michael R. Chernick,
3
Lawless i Fredette (2005), „Frequentist Prediction Intervals and Predictive Distribution”, Biometrika , 92 , 3 to kolejne dobre odniesienie, oprócz tych pod linkiem, który podałem.
Scortchi - Przywróć Monikę

Odpowiedzi:

24

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ścippBeta(α,β)α^=α+k,β^=β+nkkny , 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 jestpαβα=β=1

Pr(p|n,k)=Beta(α+k,β+nk)

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 kyy=y1,,ymmnj=0,,mjmksukcesy w poprzednich próbach; tylna predykcyjna funkcja masy:n

Pr(j|m,y)=Pr(j|m,n,k)=01Pr(j,p|m,n,k)dp=01Pr(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 prostuYpjm

f(j|m,p)=(jm)pj(1p)j

W ten sposób wyrażenie staje się

Pr(j|m,n,k)=01(jm)pj(1p)jPr(p|n,k)dp=01(jm)pj(1p)jBeta(α+k,β+nk)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)=m!j!(mj)!Γ(α+β+n)Γ(α+k)Γ(β+nk)Γ(α+k+j)Γ(β+n+mkj)Γ(α+β+n+m)

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

μ=m(α+k)(α+β+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 1j 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 = 0[j1,j2]Pr(j1jj2)=0.95mm i stąd znajdź j 1 , j 2 takie, żePr(j=0|m,n,k),Pr(j1|m,n,k),,Pr(jm1|m,n,k)j1,j2

Pr(j1jj2)=Pr(jj2|m,n,k)Pr(j<j1|m,n,k)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[j1,j2]

Pr(j=0|m,n,k)=p0,Pr(j1|m,n,k)=p1,,Pr(jm1|m,n,k)=pm1

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:p0,,pm1

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.αβ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ępnien=100k=70jm=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

j2Pr(jj2|m,n,k)0.975j1Pr(j<j1|m,n,k)=Pr(jj11|m,n,k)0.025

Pr(j1jj2|m,n,k)=Pr(jj2|m,n,k)Pr(j<j1|m,n,k)0.9750.025=0.95

j2=18j1=9Pr(j1jj2|m,n,k)0.95

Rozwiązanie dla częstych

YBinom(m,p)XBinom(n,p)12αYXI=[L(X;n,m,α),U(X;n,m,α)]

PrX,Y(YI)=PrX,Y(L(X;n,m,α)YU(X;n,m,α)]12α

12αXX+Y=k+j=ssnn+m

Pr(X=k|X+Y=s,n,n+m)=(nk)(msk)(m+ns)

XX+Y=s

Pr(Xk|s,n,n+m)=H(k;s,n,n+m)=i=0k(ni)(msi)(m+ns)

pk1αL

Pr(Xk|k+L,n,n+m)=1H(k1;k+L,n,n+m)>α

1α

Pr(Xk|k+U,n,n+m)=H(k;k+U,n,n+m)>α

[L,U]Y12αpnm12α

Praktyczny przykład z rozwiązaniem Frequentist

αβ

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

p^=knm

frequentist_point_estimate <- m * k/n #14

W przypadku przedziału przewidywania procedura jest nieco inna. Szukamy największego takiego, że P r ( X UPr(Xk|k+U,n,n+m)=H(k;k+U,n,n+m)>α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ą, żeLPr(Xk|k+L,n,n+m)=1H(k1;k+L,n,n+m)>α

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

[L,U]=[8,18]

DeltaIV
źródło