1.2 El Bootstrap uniforme
Como ya se comentó anteriormente el bootstrap uniforme es aquel en el que se reemplaza la distribución poblacional (desconocida) por la distribución empírica: \[F_n\left( x \right) =\frac{1}{n}\sum_{i=1}^{n}\mathbf{1}\left\{ X_i\leq x\right\}.\]
Es decir \(\hat{F}=F_n\) y, por lo tanto, \(R^{\ast}=R\left( \mathbf{X}^{\ast},F_n \right)\).
Conviene recordar algunas propiedades de la distribución empírica: \[\begin{aligned} nF_n\left( x \right) &= \sum_{i=1}^{n}\mathbf{1}\left\{ X_i\leq x\right\} \sim \mathcal{B}\left( n,F\left( x \right) \right), \\ E\left( nF_n\left( x \right) \right) &= nF\left( x \right) \implies E\left( F_n\left( x \right) \right) =F\left( x \right), \\ Var\left( nF_n\left( x \right) \right) &= nF\left( x \right) \left( 1-F\left( x \right) \right) \\ &\implies Var\left( F_n\left( x \right) \right) =\frac{F\left( x \right) \left( 1-F\left( x \right) \right)}{n} \end{aligned}\]
Así pues, en este caso el algoritmo bootstrap uniforme (también llamado bootstrap naïve) es el siguiente:
Para cada \(i=1,\ldots ,n\) arrojar \(X_i^{\ast}\) a partir de \(F_n\), es decir \(P^{\ast}\left( X_i^{\ast}=X_j \right) =\frac{1}{n}\), \(j=1,\ldots ,n\)
Obtener \(\mathbf{X}^{\ast}=\left( X_1^{\ast},\ldots ,X_n^{\ast} \right)\)
Calcular \(R^{\ast}=R\left( \mathbf{X}^{\ast},F_n \right)\)
Como veremos más adelante, a veces (muy poco frecuentemente) es posible calcular exactamente la distribución bootstrap de \(R^{\ast}\). Cuando eso no es posible, esa distribución es fácilmente aproximable por Monte Carlo, arrojando una gran cantidad, \(B\), de réplicas de \(R^{\ast}\). En ese caso, el algoritmo se convierte en:
Para cada \(i=1,\ldots ,n\) arrojar \(X_i^{\ast}\) a partir de \(F_n\)
Obtener \(\mathbf{X}^{\ast}=\left( X_1^{\ast},\ldots ,X_n^{\ast} \right)\)
Calcular \(R^{\ast}=R\left( \mathbf{X}^{\ast},F_n \right)\)
Repetir \(B\) veces los pasos 1-3 para obtener las réplicas bootstrap \(R^{\ast (1)}, \ldots, R^{\ast (B)}\)
Utilizar esas réplicas bootstrap para aproximar la distribución en el muestreo de \(R\)
1.2.1 Ejemplos
En el Ejemplo 1.1 anteriormente visto de inferencia para la media con varianza conocida, el algoritmo bootstrap (basado en Monte Carlo) para aproximar la distribución en el muestreo de \(R\) empleado fue:
Para cada \(i=1,\ldots ,n\) arrojar \(U_i\sim \mathcal{U}\left( 0,1 \right)\) y hacer \(X_i^{\ast}=X_{\left\lfloor nU_i\right\rfloor +1}\)
Obtener \(\bar{X}^{\ast}=\frac{1}{n}\sum X_i^{\ast}\)
Calcular \(R^{\ast}=\sqrt{n}\frac{\bar{X}^{\ast}-\bar{X}}{ \sigma }\)
Repetir \(B\) veces los pasos 1-3 para obtener las réplicas bootstrap \(R^{\ast (1)}, \ldots, R^{\ast (B)}\)
Aproximar la distribución en el muestreo de \(R\) mediante la empírica de \(R^{\ast (1)}, \ldots, R^{\ast (B)}\)
Como curiosidad podemos calcular la esperanza y la varianza de \(R\) y la esperanza y varianza bootstrap de \(R^{\ast}\). Para \(R\) tenemos: \[\begin{aligned} E\left( R \right) &=\sqrt{n}\frac{E\left( \bar{X} \right) -\mu }{\sigma } =0, \\ Var\left( R \right) &=n\frac{Var\left( \bar{X} \right)}{\sigma^2}=n \frac{\frac{1}{n}\sigma^2}{\sigma^2}=1. \end{aligned}\]
Para calcular esos mismos momentos de \(R^{\ast}\), resultará útil obtener previamente la esperanza y varianza bootstrap de \(\bar{X}^{\ast}\): \[\begin{aligned} E^{\ast}\left( \bar{X}^{\ast} \right) &= \frac{1}{n} \sum_{i=1}^{n}E^{\ast}\left( X_i^{\ast} \right) =\frac{1}{n} \sum_{i=1}^{n}E^{\ast}\left( X_1^{\ast} \right) =E^{\ast}\left( X_1^{\ast} \right) =\bar{X}, \\ Var^{\ast}\left( \bar{X}^{\ast} \right) &= \frac{1}{n^2} \sum_{i=1}^{n}Var^{\ast}\left( X_i^{\ast} \right) =\frac{1}{n^2} \sum_{i=1}^{n}Var^{\ast}\left( X_1^{\ast} \right) =\frac{1}{n}Var^{\ast }\left( X_1^{\ast} \right) =\frac{S_n^2}{n}, \end{aligned}\] ya que \[\begin{aligned} E^{\ast}\left( X_1^{\ast} \right) &= \sum_{j=1}^{n}X_jP^{\ast}\left( X_1^{\ast}=X_j \right) =\sum_{j=1}^{n}\frac{1}{n}X_j=\bar{X}, \\ Var^{\ast}\left( X_1^{\ast} \right) &= E^{\ast}\left( X_1^{\ast 2} \right) -\left[ E^{\ast}\left( X_1^{\ast} \right) \right] ^2=\sum_{j=1}^{n}X_j^2P^{\ast}\left( X_1^{\ast}=X_j \right) -\bar{X} ^2 \\ &= \frac{1}{n}\sum_{j=1}^{n}X_j^2-\bar{X}^2=\frac{1}{n} \sum_{j=1}^{n}\left( X_j-\bar{X} \right)^2=S_n^2 \end{aligned}\]
Así pues, la esperanza y la varianza bootstrap de \(R^{\ast}\) resultan: \[\begin{aligned} E^{\ast}\left( R^{\ast} \right) &= \sqrt{n}\frac{E^{\ast}\left( \bar{X}^{\ast} \right) -\bar{X}}{\sigma }=0, \\ Var^{\ast}\left( R^{\ast} \right) &= n\frac{Var^{\ast}\left( \bar{X}^{\ast} \right)}{\sigma^2}=n\frac{\frac{1}{n}S_n^2}{\sigma^2}= \frac{S_n^2}{\sigma^2}. \end{aligned}\]
Es curioso observar que la esperanza de \(R\) y la esperanza bootstrap de \(R^{\ast}\) coinciden (son ambas cero), pero no ocurre lo mismo con sus varianzas: la de \(R\) es \(1\) y la varianza bootstrap de \(R^{\ast}\) es \(S_n^2/\sigma^2\), que, aunque tiende a \(1\) (en probabilidad o de forma casi segura, bajo las condiciones adecuadas) cuando \(n\rightarrow \infty\), no es igual a \(1\). Eso nos lleva a intuir que el método de remuestreo bootstrap propuesto quizá podría modificarse ligeramente para que imitase exactamente al caso no bootstrap también en la varianza. Puede comprobarse que eso se consigue remuestreando \(X^{\ast}\) de la distribución empírica de la muestra modificada: \(\left( \tilde{X}_1,\ldots ,\tilde{X}_n \right)\), siendo \[\tilde{X}_i=\bar{X}+\frac{\sigma }{S_n}\left( X_i-\bar{X} \right) \text{, }i=1,\ldots ,n.\]
Efectivamente, bajo ese nuevo remuestreo, se tiene \[\begin{aligned} E^{\ast}\left( \bar{X}^{\ast} \right) &= E^{\ast}\left( X_1^{\ast } \right) =\overline{\tilde{X}}=\frac{1}{n}\sum_{i=1}^{n}\left[ \bar{X}+ \frac{\sigma }{S_n}\left( X_i-\bar{X} \right) \right] \\ &= \bar{X}+\frac{1}{n}\frac{\sigma }{S_n}\sum_{i=1}^{n}\left( X_i- \bar{X} \right) =\bar{X}, \\ Var^{\ast}\left( \bar{X}^{\ast} \right) &= \frac{1}{n}Var^{\ast }\left( X_1^{\ast} \right) =\frac{\sigma^2}{n}, \end{aligned}\] ya que \[\begin{aligned} Var^{\ast}\left( X_1^{\ast} \right) &= E^{\ast}\left( X_1^{\ast 2} \right) -\left[ E^{\ast}\left( X_1^{\ast} \right) \right] ^2=\sum_{j=1}^{n}\tilde{X}_j^2P^{\ast}\left( X_1^{\ast}=\tilde{X} _j \right) -\overline{\tilde{X}}^2 \\ &= \frac{1}{n}\sum_{j=1}^{n}\tilde{X}_j^2-\overline{\tilde{X}}^2=\frac{ 1}{n}\sum_{j=1}^{n}\left( \tilde{X}_j-\overline{\tilde{X}} \right)^2= \frac{1}{n}\sum_{j=1}^{n}\left[ \frac{\sigma }{S_n}\left( X_j-\bar{X} \right) \right]^2 \\ &= \frac{\sigma^2}{S_n^2}\frac{1}{n}\sum_{j=1}^{n}\left( X_i- \bar{X} \right)^2=\frac{\sigma^2}{S_n^2}S_n^2=\sigma^2. \end{aligned}\] Como consecuencia \[\begin{aligned} E^{\ast}\left( R^{\ast} \right) &= \sqrt{n}\frac{E^{\ast}\left( \bar{X}^{\ast} \right) -\bar{X}}{\sigma }=0, \\ Var^{\ast}\left( R^{\ast} \right) &= n\frac{Var^{\ast}\left( \bar{X}^{\ast} \right)}{\sigma^2}=n\frac{\frac{\sigma^2}{n}}{\sigma^2} =1. \end{aligned}\]
Esto es muy coherente con lo que nos diría la intuición pues, si la varianza poblacional, \(\sigma^2\), es conocida (ese es el motivo de que podamos usarla directamente en la definición del estadístico \(R\)), el plan de remuestreo bootstrap también ha de conocer \(\sigma^2\), es decir ha de diseñarse de modo que la distribución bootstrap de \(X^{\ast}\) tenga también varianza bootstrap \(\sigma^2\). Eso ocurre con el remuestreo uniforme de la muestra transformada \(\left( \tilde{X}_1,\ldots ,\tilde{X}_n \right)\), pero no ocurre con el remuestreo naïve (a partir de la distribución empírica de la muestra original). Esto da pie a una de las consideraciones más importantes a la hora de diseñar un buen método de remuestreo bootstrap: ha de procurarse que el bootstrap imite todas las condiciones que cumple la población original.
El código para realizar remuestreo bootstrap uniforme sobre la empírica de la muestra perturbando es análogo:
# Remuestreo
B <- 1000
estadistico_boot <- numeric(B)
coeficiente <- sigma/sd(muestra)
muestra_perturbada <- x_barra + coeficiente * (muestra - x_barra)
for (k in 1:B) {
remuestra <- sample(muestra_perturbada, n, replace = TRUE)
x_barra_boot <- mean(remuestra)
estadistico_boot[k] <- sqrt(n) * (x_barra_boot - x_barra)/sigma
}
# Aproximación bootstrap de los ptos críticos
pto_crit <- quantile(estadistico_boot, c(alfa/2, 1 - alfa/2))
# Construcción del IC
ic_inf_boot <- x_barra - pto_crit[2] * sigma/sqrt(n)
ic_sup_boot <- x_barra - pto_crit[1] * sigma/sqrt(n)
IC_boot <- c(ic_inf_boot, ic_sup_boot)
names(IC_boot) <- paste0(100*c(alfa/2, 1-alfa/2), "%")
IC_boot## 2.5% 97.5%
## 0.5024398 1.0827052Ejemplo 1.3 (Inferencia sobre la mediana)
Consideramos la mediana poblacional como parámetro de interés: \[\theta = \theta \left( F \right) = F^{-1}\left( \frac{1}{2} \right) = \inf \left\{ x\in \mathbb{R} : F\left( x \right) \geq \frac{1}{2}\right\}.\] Dada una muestra \(\mathbf{X}=\left( X_1,\ldots ,X_n \right) \sim F\), \(\theta\) puede estimarse mediante la mediana muestral \[\begin{aligned} \hat{\theta} &= \theta \left( F_n \right) =F_n^{-1}\left( \frac{1}{2} \right) =\inf \left\{ x\in \mathbb{R} : F_n\left( x \right) \geq \frac{1}{2} \right\} \\ &= \left\{ \begin{array}{ll} X_{(m)} & \text{si } n=2m-1 \text{ es impar} \\ \frac{X_{(m)}+X_{\left( m+1 \right)}}{2} & \text{si } n=2m \text{ es par} \end{array} \right. \end{aligned}\] siendo \(X_{(1)},\ldots ,X_{(n)}\) los estadísticos ordenados.
El estadístico interesante para realizar inferencia en este contexto es \(R=\sqrt{n}\left( \hat{\theta}-\theta \right)\). Si la población de partida es continua, puede demostrarse que su distribución asintótica (i.e., cuando \(n \rightarrow \infty\)) viene dada por \[R=\sqrt{n}\left( \hat{\theta}-\theta \right) \overset{d}{\rightarrow } \mathcal{N}\left( 0,\frac{1}{f\left( \theta \right)^2} \right),\]donde \(f\) es la función de densidad de la población. Como consecuencia, la utilización de esta distribución límite, \(\mathcal{N}\left( 0, 1/f\left( \theta \right)^2 \right)\), para realizar inferencia sobre la mediana, además de comportar una aproximación de la distribución en el muestreo real, no puede utilizarse directamente porque la densidad (desconocida) aparece en la expresión de la varianza asintótica. Para ser utilizable en la práctica deberíamos estimar \(f\), lo cual es un problema añadido.
Esta es pues una situación muy natural en la que usar un método bootstrap para aproximar la distribución de \(R\). Consideremos como estimador de \(F\) la distribución empírica, \(F_n\), y procedamos según un bootstrap uniforme (supongamos \(n=2m-1\), impar, por simplicidad):
Para cada \(i=1,\ldots ,n\) arrojar \(U_i\sim \mathcal{U}\left( 0,1 \right)\) y hacer \(X_i^{\ast}=X_{\left\lfloor nU_i\right\rfloor +1}\)
Obtener \(X_{(1)}^{\ast},\ldots ,X_{(n)}^{\ast}\) los estadísticos ordenados de la remuestra bootstrap y quedarse con el que ocupa lugar central: \(\hat{\theta}^{\ast}=\theta \left( F_n^{\ast} \right) =X_{(m)}^{\ast}\)
Calcular \(R^{\ast}=\sqrt{n}\left( X_{(m)}^{\ast}-X_{\left(m \right)} \right)\)
Repetir \(B\) veces los pasos 1-3 para obtener las réplicas bootstrap \(R^{\ast (1)}, \ldots, R^{\ast (B)}\)
Aproximar la distribución en el muestreo de \(R\) mediante la empírica de \(R^{\ast (1)}, \ldots, R^{\ast (B)}\)
El código implementando este algoritmo sería muy similar al de los casos anteriores:
x_mediana<- median(muestra)
# Remuestreo
B <- 1000
estadistico_boot <- numeric(B)
coeficiente <- sigma/sd(muestra)
for (k in 1:B) {
remuestra <- sample(muestra, n, replace = TRUE)
x_mediana_boot <- median(remuestra)
estadistico_boot[k] <- sqrt(n) * (x_mediana_boot - x_mediana)
}
# Aproximación bootstrap de los ptos críticos
pto_crit <- quantile(estadistico_boot, c(alfa/2, 1 - alfa/2))
# Construcción del IC
ic_inf_boot <- x_mediana - pto_crit[2]/sqrt(n)
ic_sup_boot <- x_mediana - pto_crit[1]/sqrt(n)
IC_boot <- c(ic_inf_boot, ic_sup_boot)
names(IC_boot) <- paste0(100*c(alfa/2, 1-alfa/2), "%")
IC_boot## 2.5% 97.5%
## 0.132 0.962Sin embargo, como veremos más adelante, este caso de inferencia de la mediana es uno de los pocos casos en los que la distribución bootstrap se puede calcular de forma exacta, siendo dicha expresión utilizable en la práctica.