9.2 El Bootstrap uniforme

Suponemos que X=(X1,,Xn) es una una muestra aleatoria simple (m.a.s.) de una población con distribución F y que estamos interesados en hacer inferencia sobre θ=θ(F) empleando un estimador ˆθ=T(X). Para ello nos gustaría conocer la distribución en el muestreo de un estadístico R(X,F), función del estimador (y por tanto de la muestra) y de la distribución poblacional. Por ejemplo el estimador studentizado: R=R(X,F)=ˆθθ^Var(ˆθ).

A veces podemos calcular directamente la distribución de R(X,F), aunque suele depender de cantidades poblacionales, no conocidas en la práctica. Otras veces sólo podemos llegar a aproximar la distribución de R(X,F) cuando n. Por ejemplo, bajo normalidad Xii.i.d.N(μ,σ2), si el parámetro de interés es la media θ(F)=μ=x dF(x)=xf(x) dx y consideramos como estimador la media muestral ˆθ=T(X)=θ(Fn)=x dFn(x)=ˉX. Como normalmente en la práctica la varianza es desconocida, podríamos considerar el estadístico: R=R(X,F)=nˉXμSn1tn1, donde S2n1 es la cuasivarianza muestral: S2n1=1n1nj=1(XjˉX)2. Si F no es normal entonces la distribución de R ya no es una tn1, aunque, bajo condiciones muy generales, RdN(0,1).

En el universo bootstrap se reemplaza la distribución poblacional (desconocida) F por una estimación, ˆF, de la misma. A partir de la aproximación ˆF podríamos generar, condicionalmente a la muestra observada, remuestras X=(X1,,Xn) con distribución XiˆF, que demoninaremos remuestras bootstrap. Por lo que podemos hablar de la distribución en el remuestreo de R=R(X,ˆF), denominada distribución bootstrap.

Una de las consideraciones más importantes al diseñar un buen método de remuestreo bootstrap es imitar por completo el procedimiento de muestreo en la población original (incluyendo el estadístico y las características de su distribución muestral).

Como ya se comentó anteriormente, en el bootstrap uniforme se emplea como aproximación la distribución empírica (ver Sección A.1.2): Fn(x)=1nni=11{Xix}. Es decir, ˆF=Fn, y por tanto R=R(X,Fn). En raras ocasiones (e.g. Sección 1.3 de Cao y Fernández-Casal, 2020) es posible calcular exactamente la distribución bootstrap de R. Normalmente se aproxima esa distribución mediante Monte Carlo, generando una gran cantidad, B, de réplicas de R. En el caso del bootstrap uniforme, el algoritmo es:

  1. Para cada i=1,,n generar Xi a partir de Fn, es decir P(Xi=Xj)=1n, j=1,,n

  2. Obtener X=(X1,,Xn)

  3. Calcular R=R(X,Fn)

  4. Repetir B veces los pasos 1-3 para obtener las réplicas bootstrap R(1),,R(B)

  5. Utilizar esas réplicas bootstrap para aproximar la distribución en el muestreo de R.

Para la elección del número de réplicas Monte Carlo B se aplicarían las mismas recomendaciones de la Sección 3.2 para el caso general de una aproximación por simulación.

Como ya se mostró anteriormente, el paso 1 se puede llevar a cabo simulando una distribución uniforme discreta mediante el método de la transformación cuantil (Algoritmo 5.1):

  1. Para cada i=1,,n arrojar UiU(0,1) y hacer Xi=XnUi+1

Aunque en R es recomendable33 emplear la función sample para generar muestras aleatorias con reemplazamiento del conjunto de datos original:

muestra_boot <- sample(muestra, replace = TRUE)

Ejemplo 9.2 (Inferencia sobre la media con varianza desconocida)

Como ejemplo consideramos la muestra de tiempos de vida de microorganismos lifetimes del paquete simres:

library(simres)
muestra <- lifetimes
summary(muestra)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.143   0.265   0.611   0.805   1.120   2.080
hist(muestra)
rug(muestra)
Distribución del tiempo de vida de una muestra de microorganismos.

Figura 9.6: Distribución del tiempo de vida de una muestra de microorganismos.

Supongamos que queremos obtener una estimación por intervalo de confianza de su vida media a partir de los 15 valores observados mediante bootstrap uniforme considerando el estadístico R=R(X,F)=nˉXμSn1, (en la Sección 11.2 se tratará con más detalle la construcción de intervalos de confianza).

En el bootstrap uniforme se emplea ˆF=Fn, con lo cual el análogo bootstrap del estadístico R será R=R(X,Fn)=nˉXˉXSn1, siendo ˉX=1nni=1Xi,S2n1=1n1ni=1(XiˉX)2.

El algoritmo bootstrap (aproximado por Monte Carlo) procedería así:

  1. Para cada i=1,,n arrojar UiU(0,1) y hacer Xi=XnUi+1

  2. Obtener ˉX y S2n1

  3. Calcular R=nˉXˉXSn1

  4. Repetir B veces los pasos 1-3 para obtener las réplicas bootstrap R(1),,R(B)

  5. Aproximar la distribución en el muestreo de R mediante la distribución empírica de R(1),,R(B)

Por ejemplo, podríamos emplear el siguiente código:

n <- length(muestra)
alfa <- 0.05
# Estimaciones muestrales
x_barra <- mean(muestra)
cuasi_dt <- sd(muestra)
# Remuestreo
set.seed(1)
B <- 1000
estadistico_boot <- numeric(B)
for (k in 1:B) {
  remuestra <- sample(muestra, n, replace = TRUE)
  x_barra_boot <- mean(remuestra)
  cuasi_dt_boot <- sd(remuestra)
  estadistico_boot[k] <- sqrt(n) * (x_barra_boot - x_barra)/cuasi_dt_boot
}

Las características de interés de la distribución en el muestreo de R se aproximan por las correspondientes de la distribución bootstrap de R. En este caso nos interesa aproximar los puntos críticos xα/2 y x1α/2, tales que: P(xα/2<R<x1α/2)=1α. Para lo que podemos emplear los cuantiles muestrales34:

pto_crit <- quantile(estadistico_boot, c(alfa/2, 1 - alfa/2))
pto_crit
##    2.5%   97.5% 
## -3.0022  1.8773

A partir de los cuales obtenemos la correspondiente estimación por IC boostrap: ^ICboot1α(μ)=(¯Xx1α/2Sn1n, ¯Xxα/2Sn1n).

ic_inf_boot <- x_barra - pto_crit[2] * cuasi_dt/sqrt(n)
ic_sup_boot <- x_barra - pto_crit[1] * cuasi_dt/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.50301 1.28881

Este procedimiento para la construcción de intervalos de confianza se denomina método percentil-t y se tratará en la Sección 11.2.4.

Como ejemplo adicional podemos comparar la aproximación de la distribución bootstrap del estadístico con la aproximación tn1 basada en normalidad.

hist(estadistico_boot, freq = FALSE, ylim = c(0, 0.4))
abline(v = pto_crit, lty = 2)
curve(dt(x, n - 1), add = TRUE, col = "blue")
pto_crit_t <- qt(1 - alfa/2, n - 1)
abline(v = c(-pto_crit_t, pto_crit_t), col = "blue", lty = 2)
Aproximación de la distribución de la media muestral studentizada mediante bootstrap uniforme.

Figura 9.7: Aproximación de la distribución de la media muestral studentizada mediante bootstrap uniforme.

En este caso la distribución bootstrap del estadístico es más asimétrica, por lo que el intervalo de confianza no está centrado en la media, al contrario que el obtenido con la aproximación tradicional. Por ejemplo, podemos obtener la estimación basada en normalidad mediante la función t.test():

t.test(muestra)$conf.int
## [1] 0.45994 1.15073
## attr(,"conf.level")
## [1] 0.95

En el caso multidimensional, cuando trabajamos con un conjunto de datos con múltiples variables, podríamos emplear un procedimiento análogo, a partir de remuestras del vector de índices. Por ejemplo:

data(iris)
str(iris)
## 'data.frame':    150 obs. of  5 variables:
##  $ Sepal.Length: num  5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
##  $ Sepal.Width : num  3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
##  $ Petal.Length: num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
##  $ Petal.Width : num  0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
##  $ Species     : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1..
n <- nrow(iris)
# i_boot <- floor(n*runif(n)) + 1
# i_boot <- sample.int(n, replace = TRUE)
i_boot <- sample(n, replace = TRUE)
data_boot <- iris[i_boot, ]
str(data_boot)
## 'data.frame':    150 obs. of  5 variables:
##  $ Sepal.Length: num  5.1 5.6 6.2 4.8 5.5 6.2 5.5 5.6 5 6.5 ...
##  $ Sepal.Width : num  3.8 2.5 2.9 3.1 2.3 2.9 2.6 2.8 3.6 3 ...
##  $ Petal.Length: num  1.9 3.9 4.3 1.6 4 4.3 4.4 4.9 1.4 5.8 ...
##  $ Petal.Width : num  0.4 1.1 1.3 0.2 1.3 1.3 1.2 2 0.2 2.2 ...
##  $ Species     : Factor w/ 3 levels "setosa","versicolor",..: 1 2 2 1 2..

Esta forma de proceder es la que emplea por defecto el paquete boot que describiremos más adelante (Sección 9.3.1).

Ejercicio 9.1 (Bootstrap uniforme multidimensional)
Considerando el conjunto de datos Prestige del paquete carData, supongamos que queremos realizar inferencias sobre el coeficiente de correlación entre prestige (puntuación de ocupaciones obtenidas a partir de una encuesta) e income (media de ingresos en la ocupación). Para ello podemos considerar el coeficiente de correlación lineal de Pearson: ρ=Cov(X,Y)σ(X)σ(Y) Su estimador es el coeficiente de correlación muestral: r=ni=1(xi¯x)(yi¯y)ni=1(xi¯x)2ni=1(yi¯y)2, que podemos calcular en R empleando la función cor():

data(Prestige, package = "carData")
# with(Prestige, cor(income, prestige))
cor(Prestige$income, Prestige$prestige)
## [1] 0.71491

Para realizar inferencias sobre el coeficiente de correlación, como aproximación más simple, se puede considerar que la distribución muestral de r es aproximadamente normal de media ρ y varianza Var(r)1ρ2n2.

Aproximar mediante bootstrap uniforme (multididimensional) la distribución del estadístico R=rρ, empleando B=1000 réplicas, y compararla con la aproximación normal, considerando ^Var(r)=1r2n2.


  1. De esta forma se evitan posibles problemas numéricos al emplear el método de la transformación cuantil cuando n es extremadamente grande (e.g. https://stat.ethz.ch/pipermail/r-devel/2018-September/076817.html).↩︎

  2. Se podrían considerar distintos estimadores del cuantil xα (ver p.e. la ayuda de la función quantile()). Si empleamos directamente la distribución empírica, el cuantil se correspondería con la observación ordenada en la posición Bα (se suele hacer una interpolación lineal si este valor no es entero), lo que equivale a emplear la función quantile() de R con el parámetro type = 1. Esta función considera por defecto la posición 1+(B1)α (type = 7). En el libro de Davison y Hinkley (1997), y en el paquete boot, se emplea (B+1)α (equivalente a type = 6; lo que justifica que consideren habitualmente 99, 199 ó 999 réplicas bootstrap).↩︎