6  Bootstrap

6.1 Bootstrap y el principio de sustitución (plug-in)

El bootstrap es un método muy versátil para aproximar la distribución en el muestreo de un estadístico, y en particular calcular el error típico de un estimador. Bradley Efron lo introdujo en 1979. La potencia del método procede de combinar técnicas de simulación con el principio de sustitución o plug-in. Este principio es muy simple, consiste en reemplazar la verdadera distribución de los datos por la distribución empírica muestral para estimar cualquier cantidad que dependa de $F$. Por ejemplo, en la media $\mu$ de una distribución $F$ se puede escribir $\mu =h(F)$, donde $h(F)=\int xdF(x)$. Si no conocemos $F$ pero tenemos una muestra de v.a.i.i.d. de $F$, $X_1,\dots ,X_n$, la aplicación del principio de sustitución lleva entonces al estimador $$\hat{\mu }=\int xd{F}_n(x)=\frac{1}{n}\sum _{i=1}^n{X}_i=\overline{X}.$$ En la última ecuación hemos usado que la función de distribución empírica es la distribución discreta que asigna probabilidad $1/n$ a cada dato muestral $X_i$.

La palabra bootstrap alude a una de las aventuras del Barón de Münchausen, escritas en el siglo XVIII por R. E. Raspe, según la cual el Barón cayó a las aguas de un profundo lago y consiguió salir tirando de los cordones de sus botas (de donde procede la expresión en inglés to pull oneself up by one’s own bootstrap) o de su coleta en otras versiones.

Figura 6.1: El barón saliendo del lago. (Foto: Axel Hindemith / Licencia: Creative Commons CC-by-sa-3.0)

Vamos a aplicar ahora el principio de sustitución de forma más ambiciosa, con el objetivo de estimar la distribución de un estadístico de la forma $T(X_1,\dots ,X_n;F)=\sqrt{n}(\hat{\theta }-\theta )$. Matemáticamente, podemos denotar $H_n$ a la función de distribución correspondiente: $$H_n(x)={\text{P}}_F(T(X_1,\dots ,X_n;F)\le x),$$ donde la notación ${\text{P}}_F$ indica que la probabilidad se calcula suponiendo que las v.a. $X_i$ tienen distribución $F$.

Como hemos dicho, el bootstrap esencialmente consiste en sustituir $F$ por $F_n$, con lo que resulta el estimador bootstrap ideal. La aplicación de esta sustitución lleva al estimador $${\hat{H}}_n(x)={\text{P}}_{F_n}(T(X_1^{\ast },\dots ,X_n^{\ast };F_n)\le x).$$ En la expresión anterior la distribución de las variables con asterisco es la empírica, no la verdadera. La notación es bastante habitual.

En la práctica, es necesario poder calcular ${\hat{H}}_n(x)$ de forma efectiva. Obtener una expresión cerrada suele ser imposible en general, pero al sustituir $F$ por $F_n$ pasamos del duro mundo real, en el que disponemos de una única muestra de tamaño $n$, al mundo bootstrap en el que la distribución de la que proceden los datos es totalmente conocida (la función de distribución empírica), por lo que podemos obtener tantas réplicas de $T(X_1^{\ast },\dots ,X_n^{\ast };F_n)$ como nuestra capacidad de cálculo permita. Podemos por tanto obtener muestras bootstrap o remuestras $X_1^{\ast b},\dots ,X_n^{\ast b}$ de $F_n$, donde $b=1,\dots ,B$ y $B$ es un número grande pero factible ($B\approx 1000$ suele bastar). Para ello, basta con muestrear con reemplazamiento entre los datos originales. Para cada una de estas muestras artificiales podemos calcular el estadístico de interés $T^{\ast (b)}=T^{\ast }(X_1^{\ast b},\dots ,X_n^{\ast b};F_n)$. El valor de ${\hat{H}}_n(x)$ se puede entonces aproximar de la siguiente forma: $${\hat{H}}_n(x)\approx \frac{1}{B}\sum _{b=1}^B{\mathbb{I}}_{\{T^{\ast (b)}\le x\}}.$$

En resumen, el procedimiento a seguir es:

• Se estima $F$ mediante $F_n$. (Principio de sustitución o plug-in.)
• Se obtienen $B$ muestras bootstrap $X_1^{\ast b},\dots ,X_n^{\ast b}$ ($b=1,\dots ,B$) procedentes de la distribución $F_n$, sorteando con reemplazamiento entre los datos originales $X_1,\dots ,X_n$.
• Se calcula $T^{\ast (b)}=T^{\ast }(X_1^{\ast b},\dots ,X_n^{\ast b};F_n)$ ($b=1,\dots ,B$) para cada una de las muestras bootstrap.
• Para cada $x$ fijo se calcula la proporción ${\tilde{H}}_B(x)=B^{-1}\sum _{b=1}^B{\mathbb{I}}_{\{T^{\ast (b)}\le x\}}$.

El proceso anterior implica dos aproximaciones diferentes: $$\begin{array}{}\text{(6.1)}& H_n(x)\approx {\hat{H}}_n(x)\approx {\tilde{H}}_B(x).\end{array}$$ La ley fuerte de los grandes números (aplicada a observaciones de la distribución $F_n$) garantiza que si $B\to \mathrm{\infty }$ el último término converge al segundo. En esta aproximación todo ocurre dentro del mundo bootstrap. Sin embargo, la aproximación entre los dos primeros términos cuando $n\to \mathrm{\infty }$ requiere trabajo teórico adicional. Es lo que se llama establecer la consistencia del bootstrap. En términos más intuitivos, hay que demostrar que al aumentar el tamaño muestral el mundo bootstrap se parece lo suficiente al mundo real. Dado que el teorema de Glivenko-Cantelli garantiza la convergencia uniforme de $F_n$ a $F$ con probabilidad uno, hay esperanzas fundadadas de que el bootstrap sea consistente.

El mismo algoritmo sirve esencialmente para estimar cualquier aspecto de la distribución, en lugar de la función de distribución completa. Por ejemplo, muchas veces lo que interesa es estimar la varianza o la desviación típica de un estimador $\hat{\theta }$, es decir, determinar su error típico. En este caso el estimador bootstrap ideal de la varianza de $\hat{\theta }$ es ${\text{Var}}_{F_n}({\hat{\theta }}^{\ast })$. En la práctica usaremos la correspondiente aproximación basada en $B$ remuestras: $${\text{Var}}_{F_n}({\hat{\theta }}^{\ast })\approx \frac{1}{B-1}\sum _{j=1}^B({\hat{\theta }}_j^{\ast }-{\overline{\theta }}^{\ast }{)}^2,$$ donde ${\hat{\theta }}_j^{\ast }$ es el valor del estimador para la remuestra $j$, y ${\overline{\theta }}^{\ast }=B^{-1}\sum _{j=1}^B{\hat{\theta }}_j^{\ast }$ es el promedio de todas las versiones bootstrap.

Ejemplo: el error típico de la mediana para datos de la distribución de Cauchy

Supongamos que $X_1,\dots ,X_n$ son v.a.i.i.d. de una distribución de Cauchy centrada en $\theta$ y con parámetro de escala igual a uno. La función de densidad es $$f(x)=\frac{1}{\pi }\frac{1}{1+(x-\theta {)}^2},x\in \mathbb{R}.$$ La esperanza de esta distribución no existe por lo que para estimar $\theta$ es razonable el uso de la mediana. ¿Cuál es el error típico de esta mediana? Vamos a comparar distintos métodos de calcularlo en función de la información que suponemos conocida sobre la distribución. En el código siguiente generamos una muestra original de una distribución de Cauchy con $n=30$ datos, a la que vamos a aplicar los diferentes métodos.

set.seed(100)

# Parámetros
n <- 30
theta <- 1

# Generamos los datos
muestra_original <- rt(n, 1) + theta   # Cauchy con theta = 0 coincide con t Student con 1 gl
mediana_original <- median(muestra_original)

Imaginemos primero que la distribución es conocida completamente (lo que no es muy realista). En este caso la muestra no hace falta para nada. Para aproximar la desviación típica de la mediana bastaría obtener un número grande de muestras de tamaño 30 de la ditribución de Cauchy con $\theta =1$, calcular la mediana de cada una de ellas, y la desviación típica de todas las medianas.

B <- 1000
muestras <- matrix(rt(n*B, 1) + theta, nrow = n) # cada columna una muestra
medianas <- apply(muestras, 2, median) # calcula la mediana de cada columna
sd(medianas)
#> [1] 0.2799922

En realidad no podemos suponer que conocemos la distribución. En un contexto paramétrico sí podemos suponer que la distribución es de Cauchy, pero que el valor de $\theta$ es desconocido. En esta situación podríamos repetir el esquema anterior pero sustituyendo $\theta$ por su estimador (la mediana de la muestra original):

B <- 1000
theta_estimado <- median(muestra_original)
muestras <- matrix(rt(n*B, 1) + theta_estimado, nrow = n) # cada columna una muestra
medianas <- apply(muestras, 2, median) # calcula la mediana de cada columna
sd(medianas)
#> [1] 0.2893767

El resultado obtenido es muy similar. A este procedimiento en el que se asume un modelo paramétrico $F=F_{\theta }$ y se reemplaza $F$ por $F_{\hat{\theta }}$ se le suele llamar bootstrap paramétrico.

Cuando no estamos dispuestos a asumir un modelo paramétrico, aplicamos el principio de sustitución de la forma como hemos descrito anteriormente y en lugar de generar muestras de la distribución de Cauchy, las generamos de la empírica de la muestra original:

# Generamos las remuestras (matriz n x B, cada columna una remuestra)
muestras_bootstrap <- sample(muestra_original, n*B, rep = TRUE)
muestras_bootstrap <- matrix(muestras_bootstrap, nrow = n)

# Medianas de las remuestras
medianas_bootstrap <- apply(muestras_bootstrap, 2, median)

# Estimador bootstrap del error típico de la mediana
sd_mediana <- sd(medianas_bootstrap)
sd_mediana
#> [1] 0.3355713

Vemos que se obtiene 0.336, que no es muy diferente al obtenido en el caso en que se suponía conocida la distribución o la familia paramétrica, solo unas centésimas de diferencia. Esta es la magia del bootstrap. De hecho, una vez que tenemos todas las medianas artificiales podemos usarlas para aproximar la distribución completa de la mediana muestral, no solo su error típico:

# Histograma de las medianas bootstrap
df <- data.frame(medianas_bootstrap = medianas_bootstrap)
ggplot(df) +
  geom_histogram(aes(x = medianas_bootstrap, y = ..density..),
                 bins = 10, fill = 'olivedrab4', col = 'black') +
  geom_vline(xintercept = mediana_original, size = 1.1)

La alternativa clásica al bootstrap es la teoría asintótica. Recordemos que en la sección 2.6.4 demostramos la fórmula (2.1): $$\sqrt{n}(M_n-\theta )\stackrel{\text{d}}{\to }\text{N}(0,\frac{1}{4f(\theta {)}^2}).$$

Para la distribución de Cauchy (centrada en $\theta$ y con parámetro de escala igual a 1) tenemos $f(\theta )=1/\pi$. Aplicando la fórmula, podemos aproximar la desviación típica de la mediana muestral mediante $0.5\pi /\sqrt{n}\approx 0.29$, para $n=30$. Este valor tampoco está lejos del obtenido mediante bootstrap (0.336), y es prácticamente idéntico al del bootstrap paramétrico. Al igual que este último, también hace uso de la hipótesis paramétrica de que los datos proceden de una distribución de Cauchy pero no usa el valor de $\theta$. Como curiosidad, en este caso se conoce una expresión exacta de la desviación típica de la mediana y se puede aproximar numéricamente por 0.37 (DasGupta (2008) pag. 473).

Superponemos a continuación, la distribución normal aproximada al histograma de medianas bootstrap que habíamos obtenido:

df <- data.frame(medianas_bootstrap = medianas_bootstrap)
ggplot(df) +
  geom_histogram(aes(x = medianas_bootstrap, y = ..density..),
                 bins = 10, fill = 'olivedrab4', col = 'black') +
  geom_vline(xintercept = mediana_original, size = 1.1) +
  geom_function(fun = dnorm,
                args = list(mean = mediana_original, sd = .5*pi/sqrt(n)),
                linetype = 2, size = 1.1)

Este gráfico motiva una reflexión sobre el papel desempeñado por las matemáticas en la inferencia estadística. Su papel clásico era el de proporcionar aproximaciones (TCL) o fórmulas cerradas que hicieran posibles los cálculos. En la actualidad estas aproximaciones son menos necesarias ya que como vemos se pueden sustituir por cálculos intensivos. No obstante, el papel de las matemáticas sigue siendo fundamental para interpretar los resultados, descubrir nuevos métodos y establecer las condiciones bajo las que los métodos funcionan. Por ejemplo, ¿es el bootstrap consistente en este caso? ¿Podemos demostrar que, en algún sentido, la diferencia entre el histograma y la distribución asintótica converge a cero al aumentar el tamaño muestral?

6.2 Consistencia

6.2.1 ¿Qué significa que el bootstrap es consistente?

Como muestra el ejemplo anterior, el bootstrap nos ha permitido reemplazar un TCL para la mediana, que no es trivial de obtener, por cálculos llevados a cabo con el ordenador. Sin embargo, garantizar la validez asintótica del procedimiento no es siempre fácil.

Recordemos que en general queremos garantizar que la primera aproximación de la ecuación (6.1) se cumple para $n$ grande. Una manera de formalizar esta propiedad es considerar alguna distancia entre distribuciones.

Dada una distancia entre distribuciones de probabilidad $\rho$, se dice que el bootstrap es fuertemente consistente si $\rho (H_n,{\hat{H}}_n)\to 0$ con probabilidad 1, si $n\to \mathrm{\infty }$. Es débilmente consistente si $\rho (H_n,{\hat{H}}_n)\stackrel{\text{p}}{\to }0$. En los primeros artículos que incluyeron este tipo de resultados se consideraron la distancia de Kolmogorov $\rho (F,G)=\parallel F-G{\parallel }_{\mathrm{\infty }}=\underset{x}{sup}|F(x)-G(x)|$ y la llamada distancia de Mallows.

Uno de los primeros resultados de validez asintótica del bootstrap se refiere a la estimación de la distribución de la media muestral. Se puede encontrar en (Singh, 1981). y lo enunciamos aquí sin demostración:

Teorema 6.1 Supongamos ${\text{E}}_F(X^2)<\mathrm{\infty }$ y denotemos $\mu ={\text{E}}_F(X)$, $H_n(x)={\text{P}}_F(\sqrt{n}(\overline{X}-\mu )\le x)$ y ${\hat{H}}_n(x)={\text{P}}_{F_n}(\sqrt{n}({\overline{X}}^{\ast }-\overline{X})\le x)$. Entonces $\parallel H_n-{\hat{H}}_n{\parallel }_{\mathrm{\infty }}\to 0$, con probabilidad 1.

Resultados muy generales de validez del bootstrap para la mediana se pueden encontrar en Ghosh et al (1984). Enunciamos aquí uno de ellos:

Teorema 6.2 Sea $X_1\dots ,X_n$ una muestra de v.a.i.i.d. de una distribución $F$ tal que ${\text{E}}_F|X{|}^{\alpha }<\mathrm{\infty }$, para algún $\alpha >0$. Supongamos que $F$ tiene una única mediana poblacional $\theta$ y una función de densidad $f$ continua en un entorno de $\theta$ tal que $f(\theta )>0$. Entonces $${\text{Var}}_{F_n}(\sqrt{n}(M_n^{\ast }-M_n))\stackrel{\text{c.s.}}{\to }\frac{1}{4f(\theta {)}^2}.$$

Según el resultado anterior, el estimador bootstrap de la varianza converge con probabilidad uno al valor que obtendríamos usando la distribución asintótica, pero el bootstrap no requiere conocer la función de densidad $f$.

La distribución de Cauchy verifica la condición del teorema ya que si $\alpha \in (0,1)$ y $X$ tiene distribución de Cauchy con parámetro de escala uno, entonces $\text{E}(|X-\theta {|}^{\alpha })=1/\cos (\pi \alpha /2)<\mathrm{\infty }$. Por lo tanto, el teorema anterior justifica los resultados numéricos que obtuvimos en el ejemplo de la sección anterior sobre la distribución de Cauchy.

6.2.2 ¿Cuándo falla el bootstrap?

El método bootstrap no siempre es consistente. La situación típica en la que falla es aquella en la que el estadístico $T(X_1,\dots ,X_n;F)$ no admite un teorema central del límite. Por ejemplo,

• $T(X_1,\dots ,X_n;F)=\sqrt{n}(\overline{X}-\mu )$, pero $\text{Var}(X)=\mathrm{\infty }$, por lo que no podemos usar el TCL.
• $T(X_1,\dots ,X_n;F)=\sqrt{n}(g(\overline{X})-g(\mu ))$, pero $g$ no es derivable en $\mu$, con lo que no es aplicable el método delta.
• $T(X_1,\dots ,X_n;F)=\sqrt{n}(F_n^{-1}(p)-F^{-1}(p))$, pero $f(F^{-1}(p))=0$, y no se aplican los resultados asintóticos disponibles para estadísticos de orden.
• La distribución de los datos es $F_{\theta }$ y el soporte de $F_{\theta }$ depende del parámetro.

Vamos a ver un ejemplo de esta última situación. Sean $X_1,\dots ,X_n$ v.a.i.i.d. de una distribución uniforme en el intervalo $(0,\theta )$. El estimador de máxima verosimilitud de $\theta$ es $\hat{\theta }=X_{(n)}$, el máximo de las observaciones $X_1,\dots ,X_n$. Consideramos la distribución asintótica de $T_n=n(\theta -X_{(n)})$. Si $x\ge 0$ y $n$ sufientemente grande de manera que $x/n<\theta$, $${\text{P}}_F(T_n\le x)=1-{\text{P}}_F(X_{(n)}\le \theta -x/n)=1-{\left(\frac{\theta -x/n}{\theta }\right)}^n\to 1-e^{-x/\theta },$$ si $n\to \mathrm{\infty }$. Supongamos que $\theta =1$ y, por lo tanto, $T_n\stackrel{\text{d}}{\to }\text{exp}(1)$. Consideremos ahora la versión bootstrap $T_n^{\ast }=n(X_{(n)}-X_{(n)}^{\ast })$. Para $x\ge 0$ se verifica
$${\text{P}}_{F_n}(T_n^{\ast }\le x)\ge {\text{P}}_{F_n}(T_n^{\ast }=0)={\text{P}}_{F_n}(X_{(n)}^{\ast }=X_{(n)})=1-{\left(\frac{n-1}{n}\right)}^n\to 1-e^{-1}.$$ Podemos tomar, por ejemplo, $x=0.001$ para comprobar que ${\text{P}}_F(T_n\le x)$ y ${\text{P}}_{F_n}(T_n^{\ast }\le x)$ no pueden tener el mismo límite.

6.3 Intervalos de confianza bootstrap

Una vez que hemos estimado la distribución en el muestreo de $\sqrt{n}(\hat{\theta }-\theta )$, podemos usar la estimación para deducir intervalos de confianza para $\theta$. Existe una literatura muy amplia sobre el cálculo de intervalos de confianza mediante bootstrap. Aquí vamos a revisar brevemente algunos de los métodos más conocidos.

6.3.1 Método basado en la aproximación normal

Este es el método más sencillo de usar bootstrap para calcular un intervalo de confianza. Si la distribución de $\hat{\theta }$ es aproximadamente normal con media $\theta$ y desviación típica desconocida, basta usar el error típico bootstrap del estimador en la fórmula habitual del intervalo: $${\text{IC}}_{1-\alpha }(\theta )=[\hat{\theta }\mp z_{\alpha /2}{\text{et}}_{boot}(\hat{\theta })],$$ donde ${\text{et}}_{boot}(\hat{\theta })$ es un estimador bootstrap de la desviación típica (error típico bootstrap) de $\hat{\theta }$.

Ejemplo: distribución de Cauchy (continuación)

En el caso de la muestra de una distribución de Cauchy, cuyo error típico calculamos con bootstrap en un ejemplo anterior, este intervalo quedaría de la siguiente manera:

alfa <- 0.05
z <- qnorm(1 - alfa/2)
c(mediana_original - z*sd_mediana, mediana_original + z*sd_mediana)
#> [1] 0.4382376 1.7536529

Dado que sabemos que el valor del parámetro es $\theta =1$, en este caso el intervalo contiene al parámetro. Si repetimos todo el procedimiento muchas veces para muestras independientes, aproximadamente el 95% de ellas acabaremos obteniendo un intervalo que contiene a $\theta$.

6.3.2 Remuestreo de una cantidad pivotal

Si la distribución de $\sqrt{n}(\hat{\theta }-\theta )$, que estamos denotando por $H_n(x)$, fuese totalmente conocida, entonces se podría obtener un intervalo de confianza para $\theta$ de nivel exacto $1-\alpha$ despejando $\theta$ en la ecuación siguiente: $$1-\alpha ={\text{P}}_F\{H_n^{-1}(\alpha /2)\le \sqrt{n}(\hat{\theta }-\theta )\le H_n^{-1}(1-\alpha /2)\}.$$ El intervalo de confianza correspondiente es $$(\hat{\theta }-n^{-1/2}{H}_n^{-1}(1-\alpha /2),\hat{\theta }-n^{-1/2}{H}_n^{-1}(\alpha /2)).$$ Dado que $H_n$ no es conocida, resulta natural reemplazarla por el estimador bootstrap ${\hat{H}}_n$. En la práctica, esto requiere aplicar el siguiente procedimiento

• Ordenar todos los valores simulados $\sqrt{n}({\hat{\theta }}^{\ast b}-\hat{\theta })$
• Seleccionar los percentiles que dejan una proporción de valores $\alpha /2$ a su izquierda y a su derecha
• Utilizar estos valores en la fórmula anterior en el lugar de $H_n^{-1}(\alpha /2)$ y $H_n^{-1}(1-\alpha /2)$, respectivamente.

En algunos libros a este método se le llama método bootstrap híbrido.

Ejemplo: distribución de Cauchy (continuación)

En el siguiente ejemplo se aplica el método bootstrap híbrido para calcular un intervalo de confianza (nivel de confianza nominal $1-\alpha =0.95$) para la mediana de una distribución de Cauchy. Se replica el procedimiento $m=100$ muestras de tamaño $n=30$ y se determina el número de ellas en las que el intervalo contiene al verdadero valor del parámetro:

set.seed(100)

# Parámetros
R <- 1000
n <- 30
theta <- 1
m <- 100
alfa <- 0.05

# Cálculo de los intervalos
acierto <- NULL
intervalo <- NULL
for (i in 1:m){
  muestra_original <- rt(n, 1) + theta   
  mediana_original <- median(muestra_original)
  
  muestras_bootstrap <- sample(muestra_original, n*R, rep = TRUE)
  muestras_bootstrap <- matrix(muestras_bootstrap, nrow = n)
  medianas_bootstrap <- apply(muestras_bootstrap, 2, median)
  T_bootstrap <- sqrt(n) * (medianas_bootstrap - mediana_original)
  ic_min <- mediana_original -  quantile(T_bootstrap, 1-alfa/2)/sqrt(n)
  ic_max  <- mediana_original -  quantile(T_bootstrap, alfa/2)/sqrt(n)
  intervalo <- rbind(intervalo, c(ic_min, ic_max))
  acierto <- c(acierto, ic_min < theta & ic_max > theta)
}

# Gráfico
df <- data.frame(ic_min <- intervalo[,1],
                 ic_max <- intervalo[, 2],
                 ind = 1:m,
                 acierto = acierto)
ggplot(df) +
  geom_linerange(aes(xmin = ic_min, xmax = ic_max, y = ind, col = acierto)) +
  scale_color_hue(labels = c("NO", "SÍ")) +
  geom_vline(aes(xintercept = theta), linetype = 2) +
  theme_bw() +
  labs(y = 'Muestras', x = 'Intervalos (nivel 0.95)',
       title = 'IC (método bootstrap híbrido)')

En el gráfico anterior, se han representado gráficamente los intervalos obtenidos. Se usa el color rojizo en el caso de que el intervalo no contiene al parámetro. Esto ha ocurrido en 8 de los 100 intervalos (el nivel teórico es del 5 %).

Existen otros métodos más refinados para obtener intervalos de confianza mediante bootstrap. Una alternativa es basar el intervalo bootstrap en el remuestreo de una versión estudentizada del tipo $(\hat{\theta }-\theta )/{\hat{\sigma }}_n$, donde ${\hat{\sigma }}_n$ es un estimador de la desviación típica de $\hat{\theta }$. Nótese que si se usa bootstrap para obtener este estimador serían necesarios dos niveles de remuestreo.

6.3.3 Método del percentil bootstrap

Este método se basa en usar los percentiles de los valores bootstrap generados del estimador para construir el intervalo. Más concretamente,

• Generamos ${\hat{\theta }}_1^{\ast },\dots ,{\hat{\theta }}_B^{\ast }$
• Sea ${\hat{\theta }}^{\ast }(\alpha )$ el percentil $\alpha$ de los valores bootstrap, es decir, el valor tal que $\mathrm{\#}\{b:{\hat{\theta }}_b^{\ast }\le {\hat{\theta }}^{\ast }(\alpha )\}/B=\alpha$
• El intervalo es $[{\hat{\theta }}^{\ast }(\alpha /2),{\hat{\theta }}^{\ast }(1-\alpha /2)]$

Veamos una situación en la que el uso de este método es conveniente: supongamos que existe una transformación monótona creciente $g$ que normaliza y estabiliza la varianza, es decir, tal que $$\hat{\varphi }:=g(\hat{\theta })\cong \text{N}(\varphi ,c^2),\varphi =g(\theta ),c\in \mathbb{R},$$ En este caso, un intervalo de confianza aproximado viene dado por $[g^{-1}(\hat{\varphi }-c{z}_{\alpha /2}),g^{-1}(\hat{\varphi }+c{z}_{\alpha /2})]$. El método del percentil bootstrap es una aproximación a este intervalo que no requiere conocer ni $g$ ni $c$.

Supongamos que la transformación $g$ también es efectiva para la versión bootstrap del estimador, esto es, $$\begin{array}{}\text{(6.2)}& {\hat{\varphi }}^{\ast }:=g({\hat{\theta }}^{\ast })\cong \text{N}(\hat{\varphi },c^2),\hat{\varphi }=g(\hat{\theta }),c\in \mathbb{R}.\end{array}$$ Si aplicamos el método del percentil bootstrap al parámetro transformado $\varphi$ resulta el intervalo $$[g({\hat{\theta }}^{\ast }(\alpha /2)),g({\hat{\theta }}^{\ast }(1-\alpha /2))]\approx [\hat{\varphi }-c{z}_{\alpha /2},\hat{\varphi }+c{z}_{\alpha /2}].$$ Hemos usado que la transformación $g$ es creciente y la ecuación (6.2). Como consecuencia, $$[{\hat{\theta }}^{\ast }(\alpha /2),{\hat{\theta }}^{\ast }(1-\alpha /2)]\approx [g^{-1}(\hat{\varphi }-c{z}_{\alpha /2}),g^{-1}(\hat{\varphi }+c{z}_{\alpha /2})].$$

Ejemplo: un intervalo de confianza para el coeficiente de correlación

Este es un ejemplo clásico, en el sentido de que lo eligió Efron para ilustrar su método. Nosotros lo adaptaremos a un conjunto de datos de notas en 2009 y 2010 de una prueba al final de primaria en más de 1000 colegios de la Comunidad de Madrid. El siguiente código carga todos los datos y selecciona 100 de ellos aleatoriamente:

set.seed(100)
n <- 100

colegios <- read_table("http://verso.mat.uam.es/~joser.berrendero/datos/notas.txt", 
                     locale = locale(decimal_mark = ",")) %>% 
  mutate(tipo = factor(tipo)) %>% 
  slice_sample(n = n)   # selecciona n colegios aleatoriamente

El objetivo es calcular un intervalo de confianza del coeficiente de correlación entre las notas de 2009 y 2010 a partir de los 100 colegios seleccionados. La distribución del coeficiente de correlación muestral es bastante asimétrica, por lo que tradicionalmente la inferencia se basa en la transformación $z$ de Fisher, que normaliza y estabiliza la varianza de la distribución del coeficiente: $$\hat{\varphi }=g(\hat{\rho })=\frac{1}{2}\log \frac{1+\hat{\rho }}{1-\hat{\rho }},\varphi =g(\rho ),\hat{\varphi }\cong \text{N}(\varphi ,{\sigma }^2=\frac{1}{n-3}).$$

La correlación y la correlación transformada para los colegios de la muestra son:

datos_xy <- cbind(colegios$nota09, colegios$nota10)
correlacion <- cor(datos_xy)[1,2]
correlacion_fisherz <- 0.5 * log ((1+correlacion)/(1-correlacion)) 

round(c(correlacion, correlacion_fisherz), 2)
#> [1] 0.56 0.63

A continuación, calculamos un número grande de correlaciones bootstrap y representamos la distribución bootstrap de las correlaciones transformadas junto con la aproximación normal. Con ello comprobamos que la aproximación (6.2) se cumple en este caso, es decir, ${\hat{\varphi }}^{\ast }\cong \text{N}(\hat{\varphi },{\sigma }^2=1/(n-3))$. La línea roja corresponde a la correlación para la muestra original transformada.

R <- 1000     # número de remuestras

corr_bootstrap <- replicate(R, cor(datos_xy[sample(1:n, n, rep=TRUE),])[1,2])
corr_bootstrap_fisherz <- 0.5*log((1+corr_bootstrap)/(1-corr_bootstrap))

df <- data.frame(corr_bootstrap, corr_bootstrap_fisherz)

ggplot(df) +
  geom_histogram(aes(x=corr_bootstrap_fisherz, y=..density..),
                 fill='olivedrab4',
                 col='black',
                 bins = 20) +
  labs(x = 'Correlaciones bootstrap transformadas', y = NULL) + 
  geom_vline(xintercept = correlacion_fisherz, col = 'red') + 
  geom_function(fun = dnorm,
                args = list(mean=correlacion_fisherz, sd = 1/sqrt(n-3)),
                size = 1.2)

Quizá las colas pesan un poco más de lo que deberían pero vamos a considerar suficiente el grado de aproximación. Según hemos observado, en estas circunstancias el intervalo del percentil bootstrap debería ser similar al intervalo de la transformada $z$ de Fisher. El siguiente código calcula ambos:

alpha <- 0.05  # 1 - nivel de confianza

# IC basado en el percentil bootstrap
round(c(quantile(corr_bootstrap, alpha/2), quantile(corr_bootstrap, 1-alpha/2)), 2)
#>  2.5% 97.5% 
#>  0.37  0.70

# IC basado en la transformación z de Fisher
IC_phi <- c(correlacion_fisherz - qnorm(1-alpha/2)/sqrt(n-3),
            correlacion_fisherz + qnorm(1-alpha/2)/sqrt(n-3))
IC_rho <- (exp(2*IC_phi) - 1) / (exp(2*IC_phi) + 1) 
round(IC_rho, 2)
#> [1] 0.40 0.68

Como vemos, ambos intervalos son similares pero el bootstrap presenta la gran ventaja de que no es necesario conocer la transformada ni el valor de la varianza asintótica. El coeficiente para todos los colegios es igual a 0.53 con lo que los intervalos calculados con la muestra de 100 contienen al “verdadero valor del parámetro”.

6.4 Ejercicios

Ejercicio 6.1 Se extrae una remuestra bootstrap de una muestra de $n$ observaciones $X_1,\dots ,X_n$.

1. Calcula la probabilidad de que una observación prefijada, $X_j$, no aparezca en la muestra bootstrap.
2. Calcula el límite de esta probabilidad si $n\to \mathrm{\infty }$.

Ejercicio 6.2 Sea $X_1,\dots ,X_n$ una muestra de $n$ observaciones i.i.d. de una distribución $F$ con esperanza $\mu$ y varianza ${\sigma }^2$, y sea $X_1^{\ast },\dots ,X_n^{\ast }$ una muestra de $n$ observaciones i.i.d. de la distribución empírica de la muestra original $F_n$. Calcula las siguientes cantidades:

1. ${\text{E}}_{F_n}({\overline{X}}_n^{\ast }):=\text{E}({\overline{X}}_n^{\ast }|X_1,\dots ,X_n)$
2. ${\text{E}}_F({\overline{X}}_n^{\ast })$
3. ${\text{Var}}_{F_n}({\overline{X}}_n^{\ast }):=\text{Var}({\overline{X}}_n^{\ast }|X_1,\dots ,X_n)$
4. ${\text{Var}}_F({\overline{X}}_n^{\ast })$

Ejercicio 6.3 Dada una muestra de $n$ datos diferentes, calcula en función de $n$ el número de remuestras bootstrap distintas que es posible obtener.

Aplica la expresión obtenida al caso $n=15$. ¿Qué implicación práctica tiene el resultado?

Ejercicio 6.4 Sea $X_1,\dots ,X_n$ una muestra de v.a.i.i.d. con distribución uniforme en el intervalo $(0,\theta )$. Sea $\hat{\theta }=max\{X_1,\dots ,X_n\}$ el estimador de máxima verosimilitud de $\theta$. .Sea $X_1^{\ast },\dots ,X_n^{\ast }$ una muestra de v.a.i.i.d. con distribución uniforme en $(0,\hat{\theta })$ y sea ${\hat{\theta }}^{\ast }=max\{X_1^{\ast },\dots ,X_n^{\ast }\}$. Demuestra que $n(\hat{\theta }-{\hat{\theta }}^{\ast })\stackrel{\text{d}}{\to }\text{exp}(1/\theta )$ c.s.

Este ejercicio muestra la validez asintótica del llamado bootstrap paramétrico para el estimador de máxima verosimilitud, en un caso en el que el bootstrap no paramétrico no funciona.

Ejercicio 6.5 Sean $X_1\dots ,X_n$ v.a.i.i.d. de una distribución $F$ y sean $X_{1:n}\le \cdots \le X_{n:n}$ los correspondientes estadísticos de orden. Sea $X_1^{\ast }\dots ,X_n^{\ast }$ una remuestra bootstrap y sean $X_{1:n}^{\ast }\le \cdots \le X_{n:n}^{\ast }$ los correspondientes estadísticos de orden.

1. Para $i,k\in \{1,\dots ,n\}$, demuestra que $${\text{P}}_{F_n}(X_{k:n}^{\ast }>X_{i:n})=\sum _{j=0}^{k-1}(\genfrac{}{}{0ex}{}{n}{j}){\left(\frac{i}{n}\right)}^j{(1-\frac{i}{n})}^{n-j}.$$
2. Como consecuencia, demuestra que $${\text{P}}_{F_n}(X_{k:n}^{\ast }=X_{i:n})=k(\genfrac{}{}{0ex}{}{n}{k}){\int }_{(i-1)/n}^{i/n}{t}^{k-1}(1-t{)}^{n-k}dt.$$ Este ejercicio demuestra que para estimar la varianza y la distribución de los estadísticos de orden mediante bootstrap no es necesario el remuestreo. Esto se aplica en particular a la mediana.

(Indicación: consulta las propiedades de la función beta incompleta y su relación con la función de distribución de una v.a. binomial.)

Ejercicio 6.6 Consideremos la siguiente muestra de tamaño $n=10$:

muestra <- c(1, 2, 3.5, 4, 7, 7.3, 8.6, 12.4, 13.8, 18.1)

Sea $\hat{\theta }$ la media recortada al 40% que se obtiene al eliminar los dos mayores y los dos menores datos y calcular el promedio de los 6 datos restantes. Sea ${\hat{\sigma }}_R$ el estimador bootstrap de la desviación típica de $\hat{\theta }$ basado en $R$ remuestras.

Calcula ${\hat{\sigma }}_R$ para $R=10$ y para $R=1000$ usando 10 conjuntos independientes de $R$ remuestras. Comenta los resultados.

Ejercicio 6.7 Para la muestra del ejercicio 6.6 se verifica $S^2\approx 30.84$:

muestra <- c(1, 2, 3.5, 4, 7, 7.3, 8.6, 12.4, 13.8, 18.1)
var(muestra)
#> [1] 30.84233

1. Usa bootstrap para determinar el error típico de este estimador de ${\sigma }^2$.
2. Compara el resultado con el error típico que darías si, por ejemplo, supieras que los datos proceden de una distribución normal.
3. Calcula un intervalo de confianza para ${\sigma }^2$ usando el método bootstrap híbrido. Fija $1-\alpha =0.95$.

Ejercicio 6.8 Sea $F$ una distribución con media $\mu$, varianza ${\sigma }^2$ y coeficiente de asimetría $$\gamma ={\text{E}}_F[(X-\mu {)}^3]/{\sigma }^3.$$ Genera $R=1000$ muestras de observaciones i.i.d. $X_1,\dots ,X_n$ con $X_i\equiv \text{N}(0,1)$ para $n=100$. Para cada una de ellas, calcula tres intervalos de confianza bootstrap de nivel 95% para $\gamma$ usando el método híbrido, el método normal y el método percentil. Determina el porcentaje de intervalos que contienen al parámetro en cada caso. Repite el ejercicio con muestras procedentes de una distribución exponencial de parámetro $\lambda =1$.

6.5 Referencias

James et al. (2013) presenta una perspectiva aplicada del bootstrap. Para profundizar en las relaciones del bootstrap con otros métodos se puede consultar Efron y Hastie (2016). DasGupta (2008) contiene un buen resumen de las propiedades de consistencia. En general, la comparación teórica de las propiedades de los intervalos es técnicamente complicada y se basa en desarrollos de Edgeworth. La referencia clásica para este tema es Hall (1992). Una monografía clásica sobre bootstrap es Efron y Tibshirani (1994).

• DasGupta, A. (2008). Asymptotic theory of statistics and probability. Springer.
• Efron, B. y Hastie, T. (2016). Computer Age Statistical Inference. Cambridge University Press.
• Efron, B., y Tibshirani, R.J. (1994). An introduction to the bootstrap. CRC press.
• Hall, P. (1992). The bootstrap and Edgeworth expansion. Springer.
• James, G., Witten, D., Hastie, T. y Tibshirani, R. (2013). An introduction to statistical learning. Springer.