Bootstrap

Este procedimiento es utilizado para aproximar caracteristicas de las distribucion en el muestreo de un estadistico. Para ello se genera un numero de muestras, a partir de un remuestreo de la muestra original. Su ventaja principal es que no requiere hipotesis sobre el mecanismo generador de datos

Consideremos un conjuntos de datos de la distribucion uniforme

set.seed(1)
data <- runif(50)

Usando la funcion ecdf (empirical cumulative distribution function) se obtiene una distribucion bootstrap, a continuacion se grafica la distribucion teorica y bootstap.

 # Distribución bootstrap
curve(ecdf(data)(x), ylab = "FD", type = "s", lwd = 2)
# Distribución teórica
abline(a = 0, b = 1, lty = 2) 

Las características de la distribución empírica se pueden aproximar mediante simulación, realizando repetidamente muestreo aleatorio con reemplazamiento del conjunto de datos original (manteniendo el tamaño muestral):

xboot <- sample(data, replace=TRUE)
xboot

##  [1] 0.8983897 0.4112744 0.5530363 0.2672207 0.8983897 0.7237109 0.6470602
##  [8] 0.8983897 0.1255551 0.5995658 0.3841037 0.3721239 0.5297196 0.9919061
## [15] 0.2121425 0.3841037 0.1079436 0.2655087 0.4112744 0.8983897 0.6516738
## [22] 0.5530363 0.8983897 0.4112744 0.7237109 0.2059746 0.7176185 0.6684667
## [29] 0.6870228 0.4112744 0.2672207 0.4772301 0.2672207 0.6516738 0.7893562
## [36] 0.7774452 0.7237109 0.8696908 0.6870228 0.2121425 0.8696908 0.3823880
## [43] 0.5530363 0.4772301 0.4935413 0.5297196 0.9347052 0.4820801 0.1079436
## [50] 0.7176185

Ejemplo

Considere unos datos de tamaño 100 con una distribucion normal

dat = rnorm(n = 100,mean = 0,sd = 1)

Ahora se toman 1000 replicas de bootstrap

nboot = 1000
stat.dat = mean(dat)


set.seed(1)
stat.boot = c()
for (i in 1:nboot){
  dat.boot = sample(dat, replace=TRUE)
  stat.boot = c(stat.boot,mean(dat.boot))
}
mean.boot = mean(stat.boot)  
mean.boot

## [1] -0.07732813

hist(stat.boot, freq=FALSE, ylim = c(0,14))
abline(v=mean.boot, lwd=2)
# abline(v=stat.dat)

# Distribución poblacional
curve(dnorm(x, 0, 1/sqrt(1000)), lty=2, add=TRUE)
abline(v=0, lwd=2, lty=2)

Bootstrap natural/básico:

hist(stat.boot-stat.dat, freq=FALSE, ylim = c(0,14))
abline(v=mean.boot-stat.dat, lwd=2)

# Distribución poblacional
# Distribución teórica de stat.dat - stat.teor
curve(dnorm(x, 0, 1/sqrt(1000)), lty=2, add=TRUE)   
abline(v=0, lwd=2, lty=2)

Sesgo y error estándar bootstrap

# sesgo (teor=0)
mean.boot - stat.dat 

## [1] -0.001006526

# error estándar
sd(stat.boot) 

## [1] 0.08942403

# error estándar teórico
1/sqrt(1000) 

## [1] 0.03162278

Funcion optimizada para bootstrap

boot.strap <- function(dat, nboot=1000, statistic=mean)
{
  ndat <- length(dat)
  dat.boot <- sample(dat, ndat*nboot, replace=T)
  dat.boot <- matrix(dat.boot, ncol=nboot, nrow=ndat)
  stat.boot <- apply(dat.boot, 2, statistic)
}

fstatistic <- function(dat){
  #  mean(dat)
  mean(dat, trim=0.2)
  #  median(dat)
  #  max(dat)
}

set.seed(1)
stat.dat <- fstatistic(dat)
stat.boot <- boot.strap(dat, nboot, fstatistic)

res.boot <- c(stat.dat, mean(stat.boot)-stat.dat, sd(stat.boot))
names(res.boot) <- c("Estadístico", "Sesgo", "Err.Std")
res.boot

##  Estadístico        Sesgo      Err.Std 
## -0.132291855  0.001631157  0.106516854