Vorhersageintervall für binomische Zufallsvariable

Wie lautet die (ungefähre oder genaue) Formel für ein Vorhersageintervall für eine Binomial-Zufallsvariable?

Es sei angenommen , $Y \sim \mathsf{Binom}(n, p)$ , und wir beobachten , $y$ (aus gezogenem $Y$ ). Das $n$ ist bekannt.

Unser Ziel ist es, ein Vorhersageintervall von 95% für eine neue Ziehung von $Y$ .

Die Punktschätzung ist $n\hat{p}$ , wobei p = y$\hat{p}=\frac{y}{n}$ . Ein Konfidenzintervall für p ist einfach, aber ich kann nicht eine Formel für ein Prognoseintervall für findenY. Wenn wir wusstendassp(statt p ), dann ein 95% Prognoseintervall beinhaltet nur die Quantile eines binomischenfinden. Gibt es etwas, das ich übersehen habe?$\hat{p}$$Y$$p$$\hat{p}$

Ok, lass es uns versuchen. Ich werde zwei Antworten geben - die Bayesianische, die meiner Meinung nach einfach und natürlich ist, und eine der möglichen häufigeren.

Bayes'sche Lösung

Wir gehen davon aus einem Beta vor auf , i, e., P ~ B e t einem ( α , β ) , weil die Beta-Binomialverteilung Modell - Konjugat, was bedeutet , dass die a posteriori Verteilung ist auch eine Beta - Verteilung mit Parametern α = α + k , β = β + n - k , (ich verwende k die Anzahl der Erfolge in bezeichnen n Versuchen anstelle von y ). Somit wird die Folgerung stark vereinfacht. Nun, wenn Sie einige Vorkenntnisse über die wahrscheinlichen Werte von haben$p$$p \sim Beta(\alpha,\beta)$$\hat{\alpha}=\alpha+k,\hat{\beta}=\beta+n-k$$k$$n$$y$ , Sie könnten es verwenden, um die Werte von α und β festzulegen, dh Ihre Beta-Priorität zu definieren, andernfalls könnten Sie eine einheitliche (nicht informative) Priorität mit α = β = 1 oder andere nicht informativePrioritäten annehmen(siehe zum Beispielhier). In jedem Fall ist Ihr posterior$p$$\alpha$$\beta$$\alpha=\beta=1$

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

Bei der Bayes'schen Folgerung kommt es nur auf die hintere Wahrscheinlichkeit an. Wenn Sie dies wissen, können Sie auf alle anderen Größen in Ihrem Modell schließen. Sie möchten Rückschlüsse auf die Observablen : insbesondere auf einen Vektor neuer Ergebnisse y = y 1 , … , y m , wobei m nicht unbedingt gleich n ist . Speziell wollen wir für jedes j = 0 , … , m die Wahrscheinlichkeit berechnen, genau j Erfolge in den nächsten m Versuchen zu haben, vorausgesetzt, wir haben $y$$\mathbf{y}=y_1,\dots,y_m$$m$$n$$j=0,\dots,m$$j$$m$$k$Erfolge in den vorhergehenden Versuchen; Die posteriore prädiktive Massenfunktion:$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$

Doch unsere Binomial - Modell für bedeutet , dass, bedingt auf p einen bestimmten Wert aufweist, wobei die Wahrscheinlichkeit, j Erfolge in m Studien nicht auf die letzten Ergebnisse sind abhängig: es ist einfach$Y$$p$$j$$m$

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

So wird der Ausdruck

$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$

Das Ergebnis dieses Integrals ist eine bekannte Distribution, die als Beta-Binomial-Distribution bezeichnet wird: Wenn wir die Passagen überspringen, erhalten wir den schrecklichen Ausdruck

$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)}$

Unsere Punktschätzung für bei quadratischem Verlust ist natürlich der Mittelwert dieser Verteilung, dh$j$

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

Lassen Sie uns nun nach einem Vorhersageintervall suchen. Da es sich um eine diskrete Verteilung handelt, haben wir keinen Ausdruck in geschlossener Form für , so dass P r ( j 1 ≤ j ≤ j 2 ) = 0,95 ist . Der Grund dafür ist, dass je nachdem, wie Sie ein Quantil definieren, die Quantilfunktion für eine diskrete Verteilung entweder keine Funktion oder eine diskontinuierliche Funktion ist. Dies ist jedoch kein großes Problem: Für kleines m können Sie einfach die m - Wahrscheinlichkeiten P r ( j = 0) aufschreiben$[j_1,j_2]$$Pr(j_1\leq j \leq j_2)= 0.95$$m$$m$ und von hier finden j 1 , j 2 , so dass$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$

Natürlich würden Sie mehr als ein Paar finden, also würden Sie idealerweise nach dem kleinsten suchen , so dass das oben Genannte erfüllt ist. Beachten Sie, dass$[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}$

sind nur die Werte der CMF (Cumulative Mass Function) der Beta-Binomial-Verteilung, und als solche gibt es einen Ausdruck in geschlossener Form , der sich jedoch auf die verallgemeinerte hypergeometrische Funktion bezieht und daher ziemlich kompliziert ist. Ich möchte lieber nur das R-Paket installieren extraDistrund aufrufen pbbinom, um die CMF der Beta-Binomial-Distribution zu berechnen. Wenn Sie alle Wahrscheinlichkeiten auf einmal berechnen möchten , schreiben Sie einfach:$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)

wo alphaund betasind die Werte der Parameter Ihres Beta-Prior, dh und β (also 1, wenn Sie ein einheitliches Prior vor p verwenden ). Natürlich wäre alles viel einfacher, wenn R eine Quantilfunktion für die Beta-Binomial-Verteilung bereitstellen würde, aber leider nicht.$\alpha$$\beta$$p$

Praxisbeispiel mit der Bayes'schen Lösung

Sei , k = 70 (so haben wir zunächst 70 Erfolge in 100 Versuchen beobachtet). Wir wollen eine Punktschätzung und ein 95% -Vorhersageintervall für die Anzahl der Erfolge j in den nächsten m = 20 Versuchen. Dann$n=100$$k=70$$j$$m=20$

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

wo ich auf einen einheitlichen Stand vorausgesetzt habe : Abhängig von den Vorkenntnissen für Ihre spezifische Anwendung kann dies ein guter Stand sein oder nicht. Somit$p$

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

Es ist klar, dass eine nicht ganzzahlige Schätzung für keinen Sinn ergibt. Wir könnten also einfach auf die nächste ganze Zahl (14) runden. Dann für das Vorhersageintervall:$j$

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

Die Wahrscheinlichkeiten sind

> 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

Bei gleichem-tail Wahrscheinlichkeiten Intervall, wollen wir die kleinste derart , dass P r ( j ≤ j 2 | m , n , k ) ≥ 0,975 und die größte j 1 , so daß P r ( j < j 1 | m , n , k ) = P r ( j ≤ j 1 - 1 | m , n , $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$. This way, we will have

$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$

Thus, by looking at the above probabilities, we see that $j_2=18$ and $j_1=9$. The probability of this Bayesian prediction interval is 0.9778494, which is larger than 0.95. We could find shorter intervals such that $Pr(j_1\leq j \leq j_2|m,n,k)\ge 0.95$, but in that case at least one of the two inequalities for the tail probabilities wouldn't be satisfied.

Frequentist solution

I'll follow the treatment of Krishnamoorthy and Peng, 2011. Let $Y\sim Binom(m,p)$ and $X\sim Binom(n,p)$ be independently Binominally distributed. We want a $1-2\alpha-$prediction interval for $Y$, based on a observation of $X$. In other words we look for $I=[L(X;n,m,\alpha),U(X;n,m,\alpha)]$ such that:

$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$

The "$\geq 1-2\alpha$" is due to the fact that we are dealing with a discrete random variable, and thus we cannot expect to get exact coverage...but we can look for an interval which has always at least the nominal coverage, thus a conservative interval. Now, it can be proved that the conditional distribution of $X$ given $X+Y=k+j=s$ is hypergeometric with sample size $s$, number of successes in the population $n$ and population size $n+m$. Thus the conditional pmf is

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

The conditional CDF of $X$ given $X+Y=s$ is thus

$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}}$

The first great thing about this CDF is that it doesn't depend on $p$, which we don't know. The second great thing is that it allows to easily find our PI: as a matter of fact, if we observed a value $k$ of X, then the $1-\alpha$ lower prediction limit is the smallest integer $L$ such that

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

correspondingly, the the $1-\alpha$ upper prediction limit is the largest integer such that

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

Thus, $[L,U]$ is a prediction interval for $Y$ of coverage at least $1-2\alpha$. Note that when $p$ is close to 0 or 1, this interval is conservative even for large $n$, $m$, i.e., its coverage is quite larger than $1-2\alpha$.

Practical example with the Frequentist solution

Same setting as before, but we don't need to specify $\alpha$ and $\beta$ (there are no priors in the Frequentist framework):

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

The point estimate is now obtained using the MLE estimate for the probability of successes, $\hat{p}=\frac{k}{n}$, which in turns leads to the following estimate for the number of successes in $m$ trials:

frequentist_point_estimate <- m * k/n #14

For the prediction interval, the procedure is a bit different. We look for the largest $U$ such that $Pr(X\leq k|k+U,n,n+m)=H(k;k+U,n,n+m)>\alpha$, thus let's compute the above expression for all $U$ in $[0,m]$:

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

We can see that the largest $U$ such that the probability is still larger than 0.025 is

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

Same as for the Bayesian approach. The lower prediction bound $L$ is the smallest integer such that $Pr(X\geq k|k+L,n,n+m)=1-H(k-1;k+L,n,n+m)>\alpha$, thus

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

Thus our frequentist "exact" prediction interval is $[L,U]=[8,18]$.