Intervalo de predicción para el número de parcialidad en lanzar una moneda para obtener 2 consecutivos cabezas

Question

Supongamos que tengo un posiblemente sesgada de la moneda, lo que convierte a la cabeza de una desconocida fracción $p$ del tiempo.

Al principio me tira la moneda al aire 10 veces (por ejemplo, la generación de datos = TTHTTHHTHH).

Siguiente, tengo que voltear la moneda repetidamente hasta que yo llegue a 2 cabezas en una fila. Deje que el número de lanzamientos realizar antes de llegar 2 cabezas en una fila ser $F_{HH}$.

Mi Pregunta: ¿Cómo puedo calcular frequentest intervalos de predicción para $F_{HH}$ (suponiendo que ya he realizado el 10 de coin flips)?

NOTA: a mí me parece el equivalente Bayesiano respuesta a esta pregunta es bastante sencilla (??) - pero no estoy seguro de la mejor Frequentest enfoque. Para mayor claridad, que incluyen la aproximación Bayesiana a continuación.

Enfoque Bayesiano:

Tenemos: $$P(F_{HH} | data) = \int_{0}^{1} P(F_{HH}|p)P(p|data)~dp$$. Computing $P(F_{HH}|p)$ debería ser directa (por ejemplo, utilizando una aproximación de fuerza bruta como esta: el Tiempo tomado para golpear un patrón de cabeza y la cola en una serie de monedas arrojadas )

Podemos calcular $P(p|data)$ usando el teorema de Bayes: $$P(p|data) = L(data|p)P(p)/P(data)$$ donde podemos especificar el antes de $P(p)$ (por ejemplo, un uniforme antes), calcular la probabilidad de $L(data|p)$ a partir de la distribución binomial, y $P(data)$ es sólo una normalización de la constante.

Teniendo en cuenta esto, podemos calcular la integral para obtener la distribución de $P(F_{HH}|data)$, de los cuales varios punto/intervalo de estimaciones pueden ser construidos.

Pero, ¿qué acerca de un frequentest solución, donde evitamos especificando $P(p)$?

Answer 1

Esta respuesta pone en práctica un enfoque a lo largo de las líneas de whuber del comentario, donde $p$ se estima que, ingenuamente, a partir de los 10 lanzamientos de inicio y, a continuación, 'enchufado' a $P(F_{HH}|p)$ para obtener el intervalo de predicción.

El enfoque no toma en cuenta explícitamente la incertidumbre en $p$, lo que conduce a un bajo rendimiento en algunos casos, como se muestra a continuación. Si tuviéramos muchos más de 10 lanzamientos disponibles para la estimación de $p$, luego de que este enfoque podría funcionar bien. Sería interesante saber de otros enfoques que pueden dar cuenta de la incertidumbre en $p$.

Todo el código de esta respuesta está en R.

Paso 1: Código para calcular $P(F_{HH}|p)$

En primer lugar, debemos ser capaces de calcular $P(F_{HH}|p)$. El siguiente código hace que analíticamente (desde el enfoque de simulación es muy ineficiente para las pequeñas $p$):

pmf_FHH<-function(p, Nout){
#############################################################
#
# Analytically compute the probability mass function for F_HH
# F_HH = number of coin flips required to give 2 consecutive heads
#
#
# p = probability of heads (length 1 vector)
# Nout = integer vector of values for which we want the pmf 
#

# Quick exit
if(p==0) return(Nout*0)
if(p==1) return((Nout==2))
if(max(Nout)==1) return(0*Nout)

# Recursively compute the pmf
N=max(Nout)

# Storage
PrN_T=rep(NA,N) # Probability that we got to the N'th flip without 2 consecutive heads, AND the N'th flip was a tail
PrN_2H=rep(NA,N) # Probability that the N'th flip gives 2 consecutive heads (for the first time)

# First flip
PrN_T[1]=(1-p) # Probability that we got to the first flip and it was a tail
PrN_2H[1]=0 # Can't have 2 heads on 1st flip

# Second flip
PrN_T[2] =(1- p) # Probability we get to the second flip and it was a tail
PrN_2H[2]=p*p # Probability that we get 2 heads after 2 flips

# Third flip and above
for(i in 3:length(PrN_2H)){
    # 'Probability that we got to the i'th flip, and it was a tail
    # = [1-(probability that we have terminated by i-1) ]*(1-p)
    PrN_T[i] = (1-sum(PrN_2H[1:(i-1)]))*(1-p)
    # Probability that flip i-2 was a tail, and i-1 and i were heads
    PrN_2H[i]=PrN_T[i-2]*p*p 
}

return(PrN_2H[Nout]) 
}

Para probar la función anterior y posterior uso de pruebas de los intervalos de predicción, se escribe una función para simular el lanzamiento de la moneda proceso.

sim_FHH_p<-function(p,n=round(1e+04/p**3), pattern='11'){
# Simulate many coin toss sequences, ending in the first occurrence of pattern
#
# p = probability of 1 (1=heads)
# n = number of individual tosses to sample (split into sequences ending in pattern)
# pattern = pattern to split on (1=heads)
#
# returns vector with the length of each toss sequence

# Make a data string of many coin flips e.g. '011010011'.
random_data=paste(sample(c(0,1),n,replace=T,prob=c(1-p,p)),
                  collapse="")

# Split up by occurrence of pattern, count characters, and add the number of characters in pattern.
# Each element of random_FHH gives a number of coin-tosses to get pattern
random_FHH=nchar(unlist(strsplit(random_data,pattern)))+nchar(pattern)

# The last string may not have ended in pattern. Remove it. 
random_FHH=random_FHH[-length(random_FHH)]

return(random_FHH)
}

Ahora puedo ejecutar una prueba para comprobar que la simulación y los resultados analíticos son "el mismo" (aumento de la 1e+07 conseguir un mejor acuerdo).

set.seed(1)
p=0.3
# Simulate coin-toss
qq=sim_FHH_p(p,n=1e+07, pattern='11')

Nmax=round(10/p**2) # Convenient upper limit where we check pmf_FHH
empirical_pmf=rep(NA,Nmax)
for(i in 1:Nmax) empirical_pmf[i] = (sum(qq==i)/length(qq))

png('test_analytical_relation.png',width=6,height=5,res=200,units='in')
plot(1:Nmax,empirical_pmf,main='Test of analytical relation',ylab='pmf')
points(1:Nmax,pmf_FHH(p, 1:Nmax),col='red',t='l')
legend('topright', c('Approximate empirical pmf', 'Analytical pmf'), pch=c(1,NA),lty=c(NA,1),col=c(1,2))
dev.off()

Se ve bien.

Paso 2: Código para calcular el intervalo de predicción, asumiendo $p$ es conocido.

Si p es conocida, entonces podemos utilizar directamente $P(F_{HH}|p)$ para obtener un intervalo de predicción para $F_{HH}$. Por un lado (1-$\alpha$) intervalo de predicción, sólo tenemos que pasar el (1-$\alpha$) cuantil de $P(F_{HH}|p)$. El código es:

ci_FHH<-function(p, alpha=0.1,Nmax=round(10/max(p,0.001)**2), two.sided=FALSE){
## Compute a prediction interval for FHH, assuming p
## is known exactly
##
## By default, compute 1-sided prediction interval to bound the upper values of FHH
if(p==0){
    return(c(Inf, Inf, NA, NA))
}else if(p==1){
    return(c(2, 2, 0, 1))
}else{
    cdf_FHH=cumsum(pmf_FHH(p, 1:Nmax))
    if(two.sided){ 
        lowerInd=max(which(cdf_FHH<(alpha/2)))+1 
        upperInd=min(which(cdf_FHH>(1-alpha/2)))
    }else{
        lowerInd=2 
        upperInd=min(which(cdf_FHH>(1-alpha)))
    }
    return(c(lowerInd,upperInd, cdf_FHH[lowerInd-1],cdf_FHH[upperInd]))
}

Paso 3: Probar la predicción de intervalo de cobertura Teóricamente podemos esperar que los intervalos de predicción desarrollados anteriormente a ser muy bueno si $p$ se calcula correctamente, pero tal vez muy mal si no lo es. Para probar la cobertura, la siguiente función supone el verdadero valor de $p$ es conocido, y luego repetidamente hace una estimación de $p$ basado en 10 coin flips (utilizando la fracción observada de la cabeza), y calcula un intervalo de predicción con el valor estimado de $p$.

test_ci_with_estimated_p<-function(true_p=0.5, theoretical_coverage=0.9, len_data=10, Nsim=100){

# Simulate many coin-toss experiments
simRuns=sim_FHH_p(true_p,n=1e+07)

# Simulate many prediction intervals with ESTIMATED p, and see what their
# coverage is like
store_est_p=rep(NA,Nsim)
store_coverage=rep(NA,Nsim) 
for(i in 1:Nsim){
    # Estimate p from a sample of size len_data
    mysim=rbinom(len_data,1,true_p)
    est_p = mean(mysim) # sample estimate of p
    myci=ci_FHH(est_p,alpha=(1-theoretical_coverage))

    #
    store_est_p[i]=est_p
    store_coverage[i] = sum(simRuns>=myci[1] & simRuns<=myci[2])/length(simRuns)
}
return(list(est_p=store_est_p,coverage=store_coverage, simRuns=simRuns))
}

Un par de pruebas confirman que la cobertura es casi lo correcto al $p$ se calcula correctamente, pero puede ser muy malo cuando no lo es. Las cifras muestran pruebas con el real $p$ =0,2 y 0,5 (líneas verticales), y una cobertura teórica de 0,9 (líneas horizontales). Es claro que si el estimado de $p$ es demasiado alto, entonces la predicción intervalos tienden a encubiertos, mientras que si la estimación de $p$ es demasiado bajo, sobre-cubierta, excepto si el estimado de $p$ es cero, en cuyo caso no podemos calcular cualquier intervalo de predicción (ya que con el plug-in de estimación, los jefes nunca debe ocurrir). Con sólo 10 muestras para la estimación de $p$, a menudo, la cobertura está lejos del nivel teórico.

t5=test_ci_with_estimated_p(0.5,theoretical_coverage=0.9)
t2=test_ci_with_estimated_p(0.2,theoretical_coverage=0.9)

png('test_CI.png',width=12,height=10,res=300,units='in')
par(mfrow=c(2,2))
plot(t2$est_p,t2$coverage,xlab='Estimated value of p', ylab='Coverage',cex=2,pch=19,main='CI performance when p=0.2')
abline(h=0.9)
abline(v=0.2)
plot(t5$est_p,t5$coverage,xlab='Estimated value of p', ylab='Coverage',cex=2,pch=19,main='CI performance when p=0.5')
abline(h=0.9)
abline(v=0.5)
#dev.off()

barplot(table(t2$est_p),main='Estimated p when p=0.2')
    barplot(table(t5$est_p),main='Estimated p when p=0.5')
dev.off()

En los ejemplos anteriores, la cobertura media estaba bastante cerca de la cobertura deseada cuando el verdadero $p$=0.5 (87% en comparación con el 90%), pero no tan bueno cuando es true $p$=0.2 (71% vs 90%).

# Compute mean coverage + other stats
summary(t2$coverage)
    summary(t5$coverage)