¿Por qué mi intervalo bootstrap tiene cobertura terrible?

Question

Yo quería hacer una clase de demostración donde puedo comparar un t-intervalo de un bootstrap de intervalo y calcular la probabilidad de cobertura de ambos. Yo quería que los datos provienen de una distribución desigual de manera que elegí para generar los datos como exp(rnorm(10, 0, 2)) + 1, una muestra de tamaño 10, a partir de un desplazado lognormal. Escribí un guión para dibujar 1000 muestras y, para cada muestra, calcular un 95% de intervalo y un 95% de bootstrap percentil intervalo basado en 1000 repeticiones.

Cuando ejecuto el script, ambos métodos dan muy similar intervalos y ambos tienen probabilidad de cobertura de 50-60%. Me sorprendió porque pensé que el bootstrap intervalo sería mejor.

Mi pregunta es, he

• cometió un error en el código?
• cometió un error en el cálculo de los intervalos?
• cometió un error esperar que el bootstrap intervalo para tener una mejor cobertura de las propiedades?

También, hay una manera de construir una forma más confiable de CI en esta situación?

 tCI.total <- 0
 bootCI.total <- 0
 m <- 10 # sample size
 true.mean <- exp(2) + 1

for (i in 1:1000){
 samp <- exp(rnorm(m,0,2)) + 1
 tCI <- mean(samp) + c(1,-1)*qt(0.025,df=9)*sd(samp)/sqrt(10)

 boot.means <- rep(0,1000)
 for (j in 1:1000) boot.means[j] <- mean(sample(samp,m,replace=T))
 bootCI <- sort(boot.means)[c(0.025*length(boot.means), 0.975*length(boot.means))]

 if (true.mean > min(tCI) & true.mean < max(tCI)) tCI.total <- tCI.total + 1
 if (true.mean > min(bootCI) & true.mean < max(bootCI)) bootCI.total <- bootCI.total + 1 
}
tCI.total/1000     # estimate of t interval coverage probability
bootCI.total/1000  # estimate of bootstrap interval coverage probability

Answer 1

Bootstrap diagnósticos y remedios por Canto, Davison, Hinkley Y Ventura (2006) parece ser un punto lógico de partida. Se discute múltiples formas en que el bootstrap se puede romper y - lo que es más importante aquí - ofrecer diagnósticos y posibles soluciones:

1. Los valores atípicos
2. Incorrecto remuestreo modelo
3. Nonpivotality
4. La contradicción del método bootstrap

No veo un problema con el 1, 2 y 4 en esta situación. Echemos un vistazo a las 3. Como @Ben Ogorek notas (aunque estoy de acuerdo con @Glen_b que la normalidad de discusión puede ser un arenque rojo), la validez de la secuencia de arranque depende de la pivotality de la estadística que nos interesa.

La sección 4 en Canty et al. sugiere remuestreo-dentro-de-cambia la resolución para obtener una medida de sesgo y la varianza de la estimación del parámetro dentro de cada uno de remuestreo bootstrap. Aquí está el código para replicar las fórmulas de p. 15 del artículo:

library(boot)
m <- 10 # sample size
n.boot <- 1000
inner.boot <- 1000

set.seed(1)
samp.mean <- bias <- vars <- rep(NA,n.boot)
for ( ii in 1:n.boot ) {
    samp <- exp(rnorm(m,0,2)) + 1
    samp.mean[ii] <- mean(samp)
    foo <- boot(samp,statistic=function(xx,index)mean(xx[index]),R=inner.boot)
    bias[ii] <- mean(foo$t[,1])-foo$t0
    vars[ii] <- var(foo$t[,1])
}

opar <- par(mfrow=c(1,2))
    plot(samp.mean,bias,xlab="Sample means",ylab="Bias",
        main="Bias against sample means",pch=19,log="x")
    abline(h=0)
    plot(samp.mean,vars,xlab="Sample means",ylab="Variance",
        main="Variance against sample means",pch=19,log="xy")
par(opar)

Nota: el registro de las escalas - sin registros, esto es aún más evidente. Vemos muy bien cómo la variación de la secuencia de arranque significa estimar va con la media de la secuencia de arranque de la muestra. Esto me parece bastante de una pistola humeante adjuntar la culpa a nonpivotality como culpable por el bajo nivel de confianza del intervalo de cobertura.

Sin embargo, me voy feliz admitir que uno podría seguir hasta en un montón de maneras. Por ejemplo, podríamos mirar cómo si el intervalo de confianza de un determinado bootstrap replicar incluye la media real depende de la media de la particular replicar.

Como para los remedios, Canty et al. discutir las transformaciones y los logaritmos vienen a la mente (p. ej., bootstrap y construir intervalos de confianza no para la media, pero para la media de los datos registrados), pero realmente no podía hacer el trabajo.

Canty et al. continuar para discutir cómo se puede reducir el número de interno se levanta y el resto de ruido por la importancia de muestreo y suavizado, así como añadir bandas de confianza para la dinámica de las parcelas.

Esto podría ser un divertido proyecto de tesis para un inteligente estudiante. Agradecería cualquier punteros a donde me salió mal, así como a cualquier otra literatura. Y me voy a tomar la libertad de agregar el diagnostic etiqueta a esta pregunta.

Answer 2

Aunque estoy de acuerdo con Stephan Kolassa del análisis y la conclusión, $$\hat{\mu} - \mu$$ con $\hat{\mu}$ la media de la muestra no es, definitivamente, un aproximado de pivote, permítanme hacer una observación complementaria. He investigado el uso de la $t$-estadística $$\sqrt{m} \frac{\hat{\mu} - \mu}{\hat{\sigma}}$$ junto con el arranque. El resultado fue una cobertura de alrededor de 0.8. No es la solución completa, pero es una mejora.

Luego pensé un poco más acerca de la configuración de conjunto. Con sólo 10 observaciones y extremadamente sesgada distribución, entonces es, básicamente, no imposible nonparametrically estimación de la media digamos construir intervalos de confianza con el derecho de la cobertura?

La log-normal de distribución ha considerado la media de $e^2 + 1 = 8.39$. Dado que $P(X \leq 2) = 0.84$ cuando $X \sim \mathcal{N}(0,4)$ la media es de $0.84$-cuantil de la distribución! Esto significa que la probabilidad de que todos los 10 de las observaciones son menores que la media es $0.84^{10} = 0.178$. Así que en un poco menos del 18% de los casos, la observación más grande es menor que la media. Para obtener una cobertura mayor que 0.82 necesitamos de una construcción de un intervalo de confianza para la media de la extiende más allá de la observación más grande. Tengo un tiempo difícil imaginar cómo este tipo de construcción se puede hacer (y justificado) sin supuestos previos que la distribución es muy desigual. Pero doy la bienvenida a cualquier sugerencia.

Answer 3

Yo estaba confundido acerca de esto también, y me pasé un montón de tiempo en el 1996 DiCiccio y Efron papel Bootstrap Intervalos de Confianza, sin mucho que mostrar por ello.

Lo que realmente me llevó a pensar menos de la bootstrap como un método de propósito general. Yo solía pensar que es algo que podría sacar de un apuro cuando realmente atascado. Pero he aprendido de su pequeño y sucio secreto: bootstrap intervalos de confianza se basan en la normalidad de alguna manera o de otra. Permítanme explicar.

El bootstrap le da una estimación de la distribución muestral del estimador, que es todo lo que usted podría esperar, la derecha? Pero recordemos que el vínculo clásico entre la distribución de muestreo y el intervalo de confianza se basa en el hallazgo de una fundamental cantidad. Para cualquier persona que esté oxidado, considere el caso de que $$x \sim N(\mu \sigma^2)$$ y $\sigma$ es conocido. A continuación, la cantidad de $a$z = \frac{x - \mu}{\sigma} \sim N(0,1)$$ es fundamental, es decir, su distribución no depende de $\mu$. Por lo tanto, $\Pr(-1.96 \le \frac{x - \mu}{\sigma} \le 1.96) = 0.95$ y el resto es historia.

Cuando usted piensa acerca de lo que justifica los percentiles de la distribución normal están relacionados con los intervalos de confianza, se basa por completo en esta práctica fundamental cantidad. Para una distribución arbitraria, no hay ningún vínculo teórico entre los percentiles de la distribución de muestreo e intervalos de confianza, y teniendo en raw proporciones de la distribución de muestreo bootstrap no se corte.

Así Efron del BCa (sesgo corregido) intervalos de utilizar transformaciones para llegar a la aproximación de la normalidad y bootstrap-t métodos se basan en el resultado de t-estadísticas de aproximadamente fundamental. Ahora el bootstrap puede estimar el infierno de los momentos, y siempre se puede asumir la normalidad y utilizar el estándar de +/-2*SE. Pero teniendo en cuenta todo el trabajo que iba a ir no paramétricas con el bootstrap, no parece muy justo, ¿no es así?

Answer 4

Los cálculos eran correctos, me comparan con el conocido paquete de arranque. Además he añadido la BCa-intervalo (Efron), un sesgo-versión corregida del percentil bootstrap intervalo:

for (i in 1:1000) {
  samp <- exp(rnorm(m, 0, 2)) + 1

  boot.out <- boot(samp, function(d, i) sum(d[i]) / m, R=999)
  ci <- boot.ci(boot.out, 0.95, type="all")

  ##tCI <- mean(samp) + c(1,-1)*qt(0.025,df=9)*sd(samp)/sqrt(10)
  tCI <- ci$normal[2:3]
      percCI <- ci$perc[4:5]
  bcaCI <- ci$bca[4:5]
      boottCI <- ci$student[4:5]

  if (true.mean > min(tCI) && true.mean < max(tCI)) tCI.total <- tCI.total + 1
  if (true.mean > min(percCI) && true.mean < max(percCI)) percCI.total <- percCI.total + 1 
  if (true.mean > min(bcaCI) && true.mean < max(bcaCI)) bcaCI.total <- bcaCI.total + 1
}

tCI.total/1000     # estimate of t interval coverage probability
0.53
percCI.total/1000  # estimate of percentile interval coverage probability
0.55
bcaCI.total/1000  # estimate of BCa interval coverage probability
0.61

Supongo que los intervalos sería mucho mejor si el tamaño de la muestra es mayor que 10, por ejemplo, 20 o 50.

Además, el bootstrap-t método por lo general conduce a mejores resultados para sesgada estadísticas. Sin embargo, se necesita un bucle anidado y, por tanto, más de 20 veces más tiempo de cálculo.

Para la prueba de hipótesis también es muy importante que el 1 cara coberturas son buenas. Así que mirando sólo a 2 caras coberturas a menudo puede ser engañosa.