Intervalo de predicción para variable aleatoria binomial

14

¿Cuál es la fórmula (aproximada o exacta) para un intervalo de predicción para una variable aleatoria binomial?

Suponga que YBinom(n,p) , y observamos y (extraído de Y ). El n es conocido.

Nuestro objetivo es obtener un intervalo de predicción de 95% para un nuevo sorteo de Y .

La estimación puntual es np^ , dondep^=yn . Un intervalo de confianza para p es sencillo, pero no puede encontrar una fórmula para un intervalo de predicción para Y . Si supiéramos p (en lugar de p ), entonces un intervalo de predicción del 95% simplemente consiste en encontrar los cuantiles de un binomio. ¿Hay algo obvio que estoy pasando por alto?p^Ypp^

Statseeker
fuente
1
Vea ¿Qué métodos no bayesianos existen para la inferencia predictiva? . En este caso, el método que utiliza pivotes no está disponible (no lo creo), pero podría usar una de las probabilidades predictivas. O, por supuesto, un enfoque bayesiano.
Scortchi - Restablece a Monica
1
Hola chicos, me gustaría tomarme un momento para abordar las inquietudes que surgieron. - con respecto a la confianza para p: no estoy interesado en eso. - considerando que las predicciones son del 95% de la distribución: sí, eso es exactamente lo que son los intervalos de predicción independientemente del contexto (en la regresión, debe asumir errores normales, mientras que los intervalos de confianza dependen de CLT - sí, el ejemplo de predicción del número de cabezas en un lanzamiento de moneda es correcto. Lo que hace que este problema sea difícil es que ahora no "p", solo tenemos una estimación.
Statseeker
3
@Addison Lea el libro Intervalos estadísticos de G. Hahn y W. Meeker. Explican la diferencia entre los intervalos de confianza, los intervalos de predicción, los intervalos de tolerancia y los intervalos bayesianos creíbles. Un intervalo de predicción del 95% no contiene el 95% de la distribución. Hace lo que hacen los intervalos más frecuentes. Si toma muestras repetidamente de B (n, p) y usa el mismo método cada vez para producir un intervalo de predicción del 95% para p, entonces el 95% de los intervalos de predicción contendrá el valor verdadero de p. Si desea cubrir el 95% de la distribución, construya un intervalo de tolerancia.
Michael R. Chernick
Los intervalos de tolerancia cubren un porcentaje de la distribución. Para un intervalo de tolerancia del 95% para el 90% de la distribución, repita el proceso muchas veces y use el mismo método para generar el intervalo cada vez, luego en aproximadamente el 95% de los casos, al menos el 90% de la distribución caerá en el intervalo y 5% del tiempo menos del 90% de la distribución estará contenida en el intervalo.
Michael R. Chernick
3
Lawless & Fredette (2005), "Intervalos de predicción frecuente y distribuciones predictivas", Biometrika , 92 , 3 es otra buena referencia, además de las del enlace que di.
Scortchi - Restablece a Monica

Respuestas:

24

Ok, intentemos esto. Daré dos respuestas: la bayesiana, que en mi opinión es simple y natural, y una de las posibles frecuentas.

Solución bayesiana

Suponemos una previa Beta en , i, e., P ~ B e t un ( α , β ) , porque el modelo Beta-Binomial es conjugado, lo que significa que la distribución posterior es también una distribución Beta con parámetros α = α + k , β = β + 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 tiene algún conocimiento previo sobre los valores probables deppBeta(α,β)α^=α+k,β^=β+nkkny , podría usarlo para establecer los valores de α y β , es decir, para definir su beta anterior, de lo contrario podría asumir un previo uniforme (no informativo), con α = β = 1 , u otros previos no informativos (ver por ejemploaquí) En cualquier caso, su posterior espαβα=β=1

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

En la inferencia bayesiana, todo lo que importa es la probabilidad posterior, lo que significa que una vez que lo sepa, puede hacer inferencias para todas las demás cantidades en su modelo. Desea hacer inferencia en los observables : en particular, en un vector de nuevos resultados y = y 1 , ... , y m , donde m no es necesariamente igual a n . Específicamente, para cada j = 0 , ... , m , queremos calcular la probabilidad de tener exactamente j éxitos en las siguientes m pruebas, dado que obtuvimos kyy=y1,,ymmnj=0,,mjmkéxitos en las pruebas anteriores ; La función de masa predictiva posterior: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

Sin embargo, nuestro modelo binomial para significa que, condicionalmente si p tiene un cierto valor, la probabilidad de tener j éxitos en m ensayos no depende de resultados pasados: es simplementeYpjm

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

Así la expresión se convierte

Pr(j|m,n,k)=01(jm)pj(1p)jPr(p|n,k)dp=01(jm)pj(1p)jBeta(α+k,β+nk)dp

El resultado de esta integral es una distribución bien conocida llamada distribución Beta-Binomial: omitiendo los pasajes, obtenemos la horrible expresión

Pr(j|m,n,k)=m!j!(mj)!Γ(α+β+n)Γ(α+k)Γ(β+nk)Γ(α+k+j)Γ(β+n+mkj)Γ(α+β+n+m)

Nuestra estimación puntual para , dada la pérdida cuadrática, es, por supuesto, la media de esta distribución, es decir,j

μ=m(α+k)(α+β+n)

Ahora, busquemos un intervalo de predicción. Como esta es una distribución discreta, no tenemos una expresión de forma cerrada para , de modo que P r ( j 1j j 2 ) = 0.95 . La razón es que, dependiendo de cómo defina un cuantil, para una distribución discreta, la función del cuantil no es una función o es una función discontinua. Pero este no es un gran problema: para m pequeño , simplemente puede escribir las probabilidades m P r ( j = 0[j1,j2]Pr(j1jj2)=0.95mm y desde aquí encontrar j 1 , j 2 tal quePr(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

Por supuesto, encontraría más de una pareja, por lo que idealmente buscaría la más pequeña modo que se satisfaga lo anterior. Tenga en cuenta que[j1,j2]

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

son solo los valores de la CMF (función de masa acumulada) de la distribución beta-binomial, y como tal hay una expresión de forma cerrada , pero esto es en términos de la función hipergeométrica generalizada y, por lo tanto, es bastante complicado. Prefiero instalar el paquete R extraDistry llamar pbbinompara calcular el CMF de la distribución Beta-Binomial. Específicamente, si desea calcular todas las probabilidades de una vez, simplemente escriba:p0,,pm1

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

donde alphay betason los valores de los parámetros de su Beta anterior, es decir, y β (por lo tanto, 1 si está utilizando un previo uniforme sobre p ). Por supuesto, todo sería mucho más simple si R proporcionara una función cuantil para la distribución Beta-Binomial, pero desafortunadamente no lo hace.αβp

Ejemplo práctico con la solución bayesiana

Sea , k = 70 (por lo tanto, inicialmente observamos 70 éxitos en 100 ensayos). Queremos una estimación puntual y un intervalo de predicción del 95% para el número de éxitos j en los próximos m = 20 ensayos. Luegon=100k=70jm=20

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

donde asumí un uniforme previo en : dependiendo del conocimiento previo para su aplicación específica, este puede o no ser un buen previo. Asíp

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

Claramente, una estimación no entera para no tiene sentido, por lo que podríamos redondear al entero más cercano (14). Luego, para el intervalo de predicción:j

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 un intervalo de probabilidades de igual cola, queremos el más pequeño tal que P r ( j j 2 | m , n , k ) 0.975 y el j 1 más grande tal que P r ( j < j 1 | m , n , k ) = P r ( j j 1 - 1 | m , n , kj2Pr(jj2|m,n,k)0.975j1 . De esta manera, tendremosPr(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

Por lo tanto, al observar las probabilidades anteriores, vemos que y j 1 = 9 . La probabilidad de este intervalo de predicción bayesiano es 0.9778494, que es mayor que 0.95. Podríamos encontrar intervalos más cortos de modo que P r ( j 1j j 2 | m , n , k ) 0.95 , pero en ese caso al menos una de las dos desigualdades para las probabilidades de cola no se satisfaría.j2=18j1=9Pr(j1jj2|m,n,k)0.95

Solución frecuente

Seguiré el tratamiento de Krishnamoorthy y Peng, 2011 . Sea y X B i n o m ( n , p ) independientemente distribuidos binominalmente. Queremos un 1 - 2 α - intervalo de predicción para Y , en base a una observación de X . En otras palabras, buscamos I = [ L ( X ; nYBinom(m,p)XBinom(n,p)12αYX tal que:I=[L(X;n,m,α),U(X;n,m,α)]

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

The "12α" 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)=(nk)(msk)(m+ns)

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

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

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α lower prediction limit is the smallest integer L such that

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

correspondingly, the the 1α upper prediction limit is the largest integer such that

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

Thus, [L,U] is a prediction interval for Y of coverage at least 12α. 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 12α.

Practical example with the Frequentist solution

Same setting as before, but we don't need to specify α and β (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, p^=kn, 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(Xk|k+U,n,n+m)=H(k;k+U,n,n+m)>α, 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(Xk|k+L,n,n+m)=1H(k1;k+L,n,n+m)>α, 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].

DeltaIV
fuente