Ok, vamos a probar esto. Te voy a dar dos respuestas: el Bayesiano, que es, en mi opinión, simple y natural, y una de las posibles frecuentista.
Bayesiano solución
Asumimos una Beta antes de la en $p$, i,e., $p \sim Beta(\alpha,\beta)$, debido a que el Beta-modelo Binomial es conjugada, lo que significa que la distribución posterior es también una distribución Beta con parámetros de $\hat{\alpha}=\alpha+k,\hat{\beta}=\beta+n-k$, (estoy usando $k$ para denotar el número de éxitos en $n$ ensayos, en lugar de $y$). Por lo tanto, la inferencia se simplifica enormemente. Ahora, si usted tiene algún conocimiento previo sobre los posibles valores de $p$, que se podría utilizar para establecer los valores de $\alpha$$\beta$, es decir, para definir su Beta antes, de lo contrario se podría asumir un uniforme (noninformative) antes, con $\alpha=\beta=1$ u otras noninformative priores (ver por ejemplo aquí). En cualquier caso, su posterior es
$Pr(p|n,k)=Beta(\alpha+k,\beta+n-k)$
En la inferencia Bayesiana, todo lo que importa es la probabilidad posterior, lo que significa que una vez que usted sabe que usted puede hacer inferencias para todas las otras cantidades en el modelo. Usted desea hacer inferencia sobre las características observables $y$: en particular, en un vector de nuevos resultados $\mathbf{y}=y_1,\dots,y_m$ donde $m$ no es necesariamente igual a $n$. Específicamente, para cada una de las $j=0,\dots,m$, queremos calcular la probabilidad de tener exactamente $j$ éxitos en la próxima $m$ ensayos, dado que tenemos a $k$ éxitos en la anterior $n$ ensayos; la posterior predictivo de la función de masa:
$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$
Sin embargo, nuestro modelo Binomial para $Y$ significa que, condicionados en $p$ tener un cierto valor, la probabilidad de tener $j$ éxitos en $m$ ensayos no depende de los resultados del pasado: es simplemente
$f(j|m,p)=\binom{j}{m} p^j(1-p)^j$
Así, la expresión se convierte en
$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$
El resultado de esta integral es un conocido distribución llamada Beta-Binomial distribución: omitiendo los pasajes, tenemos la horrible expresión
$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)}$
Nuestra estimación para $j$, dada la pérdida cuadrática, es, por supuesto, la media de esta distribución, es decir,
$\mu=\frac{m(\alpha+k)}{(\alpha+\beta+n)}$
Ahora, vamos a buscar un intervalo de predicción. Dado que esta es una distribución discreta, no tenemos una forma cerrada de expresión para $[j_1,j_2]$, de tal manera que $Pr(j_1\leq j \leq j_2)= 0.95$. La razón es que, dependiendo de cómo se defina un cuantil, para una distribución discreta de la función cuantil es o no una función o es una función discontinua. Pero esto no es un gran problema: para las pequeñas $m$, usted puede simplemente escribir el $m$ probabilidades de $Pr(j=0|m,n,k),Pr(j\leq 1|m,n,k),\dots,Pr(j \leq m-1|m,n,k)$ y a partir de aquí encontrar $j_1,j_2$ tal que
$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$
Por supuesto, usted podría encontrar más de una pareja, por lo que lo ideal sería la mirada de los más pequeños a $[j_1,j_2]$ tales que el de arriba es satisfecho. Tenga en cuenta que
$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}$
son sólo los valores de la CMF (Acumulativas de Función de Masa) de la Beta-distribución Binomial, y como tal, hay una forma cerrada de expresión, pero esto es en términos de la función hipergeométrica generalizada y por lo tanto es bastante complicado. Prefiero solo tienes que instalar el paquete de R extraDistr y llame a pbbinom para el cálculo de la CMF de la Beta-distribución Binomial. Específicamente, si se desea calcular todas las probabilidades $p_0,\dots,p_{m-1}$ en uno, sólo tienes que escribir:
library(extraDistr)
jvec <- seq(0, m-1, by = 1)
probs <- pbbinom(jvec, m, alpha = alpha + k, beta = beta + n - k)
donde alpha y beta son los valores de los parámetros de su Beta antes, es decir, $\alpha$ $\beta$ (por lo que 1 si está usando un uniforme antes de $p$). Por supuesto, todo sería mucho más sencillo si R proporciona un cuantil de la función de la Beta-Binomial distribución, pero por desgracia no es así.
Ejemplo práctico con la solución Bayesiana
Vamos $n=100$, $k=70$ (así, hemos observado inicialmente 70 éxitos en 100 ensayos). Queremos una estimación de punto y un 95%-intervalo de predicción para el número de éxitos $j$ en los próximos $m=20$ ensayos. Entonces
n <- 100
k <- 70
m <- 20
alpha <- 1
beta <- 1
donde me supone un uniforme antes de $p$: dependiendo del conocimiento previo para su aplicación específica, esto puede o puede no ser una buena previa. Así
bayesian_point_estimate <- m * (alpha + k)/(alpha + beta + n) #13.92157
Claramente un no-entero estimación de $j$ no tiene sentido, así que simplemente redondear al entero más cercano (14). Entonces, para la predicción del intervalo:
jvec <- seq(0, m-1, by = 1)
library(extraDistr)
probabilities <- pbbinom(jvec, m, alpha = alpha + k, beta = beta + n - k)
Las probabilidades son
> 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
Para una igualdad de cola probabilidades de intervalo, queremos que los más pequeños $j_2$ tal que $Pr(j\leq j_2|m,n,k)\ge 0.975$ y el mayor $j_1$ tal que $Pr(j < j_1|m,n,k)=Pr(j \le j_1-1|m,n,k)\le 0.025$. De esta manera, tendremos
$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$
Por lo tanto, mirando a la anterior probabilidades, podemos ver que $j_2=18$$j_1=9$. La probabilidad de que este Bayesiano intervalo de predicción es 0.9778494, que es mayor que 0.95. Podríamos encontrar más cortos intervalos tales que $Pr(j_1\leq j \leq j_2|m,n,k)\ge 0.95$, pero en caso de que al menos una de las dos desigualdades de la cola probabilidades no estaría satisfecho.
Frecuentista solución
Yo voy a seguir el tratamiento de Krishnamoorthy y Peng, 2011. Deje $Y\sim Binom(m,p)$ $Y\sim Binom(n,p)$ independientemente Binominally distribuido. Queremos una $1-2\alpha-$intervalo de predicción para $Y$, basado en una observación de $X$. En otras palabras, buscamos $I=[L(X;n,m,\alpha),U(X;n,m,\alpha)]$ tal forma que:
$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$
El "$\geq 1-2\alpha$" es debido al hecho de que estamos tratando con un discreto de la variable aleatoria, y por lo tanto no podemos esperar para obtener la cobertura exacta...pero podemos buscar un intervalo que tiene siempre al menos la cobertura nominal, por lo tanto un conservador de intervalo. Ahora, se puede demostrar que la distribución condicional de $X$ $X+Y=k+j=s$ es hipergeométrica con el tamaño de la muestra $s$, el número de éxitos en la población $n$ y el tamaño de la población $n+m$. Por lo tanto el condicional pmf es
$Pr(X=k|X+Y=s,n,n+m)=\frac{\binom{n}{k}\binom{m}{s-k}}{\binom{m+n}{s}}$
El condicional CDF de $X$ $X+Y=s$ es así
$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}}$
El primer gran cosa acerca de este CDF es que no dependen de $p$, lo que no sabemos. El segundo gran cosa es que permite encontrar fácilmente nuestro PI: como una cuestión de hecho, si observamos un valor de $k$ de X, entonces la $1-\alpha$ inferior predicción límite es el entero más pequeño $L$ tal que
$Pr(X\geq k|k+L,n,n+m)=1-H(k-1;k+L,n,n+m)>\alpha$
en consecuencia, la de la $1-\alpha$ superior predicción límite es el mayor entero tal que
$Pr(X\leq k|k+U,n,n+m)=H(k;k+U,n,n+m)>\alpha$
Por lo tanto, $[L,U]$ es un intervalo de predicción para $Y$ de cobertura de al menos $1-2\alpha$. Tenga en cuenta que cuando se $p$ está cerca de 0 o 1, de este intervalo de tiempo es conservador, incluso para las grandes $n$, $m$, es decir, su cobertura es bastante más grande que la de $1-2\alpha$.
Ejemplo práctico con el Frecuentista solución
Misma configuración que antes, pero no necesitamos especificar $\alpha$ $\beta$ (no hay antecedentes en la Frecuentista marco):
n <- 100
k <- 70
m <- 20
La estimación puntual es ahora obtenidos mediante el MLE estimación de la probabilidad de éxitos, $\hat{p}=\frac{k}{n}$, que a su vez lleva a la siguiente estimación para el número de éxitos en $m$ ensayos:
frequentist_point_estimate <- m * k/n #14
Para la predicción de intervalo, el procedimiento es un poco diferente. Buscamos el mayor $U$ tal que $Pr(X\leq k|k+U,n,n+m)=H(k;k+U,n,n+m)>\alpha$, por lo tanto vamos a calcular la expresión anterior para todos los $U$$[0,m]$:
jvec <- seq(0, m, by = 1)
probabilities <- phyper(k,n,m,k+jvec)
Podemos ver que el mayor $U$ tal que la probabilidad es aún mayor que 0.025 es
jvec[which.min(probabilities > 0.025) - 1] # 18
Mismo como para el enfoque Bayesiano. La parte inferior de la predicción enlazado $L$ es el entero más pequeño tal que $Pr(X\geq k|k+L,n,n+m)=1-H(k-1;k+L,n,n+m)>\alpha$, lo que
probabilities <- 1-phyper(k-1,n,m,k+jvec)
jvec[which.max(probabilities > 0.025) - 1] # 8
Por lo tanto nuestra frecuentista "exacto" intervalo de predicción es $[L,U]=[8,18]$.
