Métodos de Montecarlo y Estadística Computacional

Con mi mayor tributo de agradecimiento y admiración a uno de mis mejores maestros, el doctor Jordi Ocaña del departamento de Estadística de la Facultad de Biología de la Universidad de Barcelona. El ha dedicado parte de su valioso tiempo, sustraído a sus investigaciones para formarnos en la disciplina. Solicito su benevolencia como excelente persona que es además de gran científico, para esas páginas, si llegara a leerlas, considerando que no es el conocimiento sino la enorme ilusión que en mi ha despertado esa bella disciplina lo que me puede servir de excusa para escribir sobre ella.

  1. Introducción. Bootstrap paramétrico y no paramétrico

  2. Intervalos de confianza bootstrap

  3. Otros métodos. Jacknife

  4. Remuestreo y contraste de hipótesis

  5. Test de permutaciones

  6. Métodos de Montecarlo  aplicados a modelos de Markov

 

  1. Introducción

La palabra bootstrap según comenta el propio Efron en su manual "An Introduction to the boostrap" redactado junto con el doctor Tibshirani, también de Stanford proviene del relato alemán Aventuras del Barón de Munchausen, donde se relata que estando el barón en el fondo de un lago, logro salir del mismo impulsándose tirando de los cordones de sus botas cuando todo parecía perdido. La primera sensación que tiene un estudiante en contacto con el bootstrap es precisamente de que se trata de un método igual de absurdo que el empleado por el original barón. Esta sensación primera se convierte luego en admiración  y finalmente en una especie de enamoramiento de un método cuyos resultados lo merecen. 

Fue Efron en el año 1979 quien creó la técnica que venia a mejorar el método Jacknife creado por Quenouville en el año  1949 seguido de  Tukey en 1958. Característica básica del método bootstrap es el pincipio "plug-in" que consiste en la substitución de la función subyacente de distribución desconocida F por un estimador  de la misma. Asimismo se emplea el remuestreo con substitución para obtener via computador gran número de remuestras sobre las que tener la base de estimación. 

 El bootstrap ofrece una aproximación para la inferencia y el contraste de hipótesis, mediante el uso de un remuestreo que nos permite estimar una distribución empírica del estadístico. En dicha técnica el poder computacional está sustituyendo al análisis teórico. El proceso se desarrolla mediante los pasos siguientes

  1.  Construir una distribución de probabilidad empírica (x), a partir de la muestra asignando una probabilidad de 1/n a cada punto, Xi, i= 1,...,n.

  2. Simular una muestra aleatoria simple de tamaño n con reposición. Esta es una remuestra bootstrap.

  3. Calcular el estadístico de interés, , a partir de esas remuestras, dando .

  4. Repetir los pasos 2 y 3 B veces, donde B es un número grande. La magnitud de B en la práctica depende de las pruebas que se van a aplicar a los datos. En general, B debería ser de entre 50 a 200 para estimar el error típico de  y, según Hall (1986), al menos de 1000 para estimar los intervalos de confianza alrededor de .

  5. Construir una distribución de probabilidad a partir de los B  asignando una probabilidad de 1/B a cada punto.

Un ejemplo de aplicación de esta  técnica que constituyó el primer trabajo de nuestro doctorado en bioestadística:

Bootstrap paramétrico y no paramétrico

Enunciado

Supongamos que la variable aleatoria X sigue una distribución uniforme entre 0 y un valor b positivo, desconocido, X~U(0, b). Supongamos muestras de tamaño n, asociadas a observaciones independientes e idénticamente distribuidas de X: X1, …,Xn. Un estimador razonable de b es e  (por ejemplo, es el estimador de máxima verosimilitud). Vamos a intentar estimar la distribución en el muestreo de este estimador mediante el método bootstrap.

1)      En primer lugar, convendría conocer cual es la verdadera distribución en el muestreo de , para así poder compararla con las estimaciones bootstrap que obtengamos. Podemos enfocar esta primera cuestión de dos maneras, en realidad complementarias:

a)      Es fácil determinar teóricamente cuál es esta distribución, es decir, determinar la distribución de  para n variables aleatorias independientes, todas ellas con distribución U(0, b). Intentad determinarla. Representad gráficamente la función de densidad obtenida, para algunos valores de b concretos, por ejemplo 5 y 10, y para algunos tamaños muestrales n también concretos, por ejemplo 10 y 50.

b)      Tanto si hemos completado el paso anterior con éxito como si no, otra posibilidad es hacernos una idea de cuál es esta distribución mediante simulación, es decir, aplicando lo que hemos llamado “el principio de Montecarlo”. Intentad realizar dicha simulación, para los mismos valores del parámetro y del tamaño muestral anteriores, y representad la estimación de la función de densidad (p.e. función “density” de R). Comparadla con la teórica, si la habéis obtenido. Deberían ser muy parecidas.

2)      Generad una muestra aleatoria simple para cada una de las combinaciones del valor del parámetro y del tamaño muestral anteriores. Cada una de estas muestras serán nuestros “datos” –que en este caso sabremos perfectamente de donde proceden, lógicamente en un caso real no sería así.

3)      Estimad la distribución muestral mediante el método bootstrap, para cada una de las cuatro muestras por separado. Representad gráficamente la densidad obtenida y comparadla con la verdadera distribución muestral. En cada caso, el remuestreo bootstrap puede ser (probad los dos enfoques):

a)      No paramétrico.

b)      Paramétrico (utilizando la información adicional de que son uniformes).

Las distribuciones bootstrap obtenidas ¿Realmente estiman de forma adecuada la verdadera distribución en el muestreo de ?

 

  1. Descripción del trabajo llevado a cabo:

    Sabemos que la variable aleatoria X sigue una distribución uniforme de parámetros 0 y b

    X ~ U  (0, b)

    Donde  b  > 0

    Tomando muestras de tamaño n, de observaciones independientes e idénticamente distribuidas de X, estimaremos  el valor de b, representado por el estimador :

     = max (X1, X2, …, Xn)

    A cuyo fin emplearemos el método bootstrap

    1.      “Verdadera” distribución  de  en el muestreo.

    a)      Determinación teórica de dicha distribución:

    Por definición tenemos que:

     = max (X1, X2, …, Xn) para n variables aleatorias independientes, donde es    U ~  (0,b)

     

    Sea D la distribución de probabilidad de .

    D (b) = P ( ≤ b)  = P (max {X1, X2, …, Xn} ≤ b)

    Lo que implica que:

                                  = P (X1 ≤ b) · P (X2 ≤ b) · … · P (Xn ≤ b)

                                     = F(X1) · F(X2) · … · F(Xn)

    o sea

                                        = ( x/ b )n

    De donde podemos hallar la función de densidad derivando respecto de x

     

    siendo b :   0 ≤ x ≤ b

    Una vez deducida la “verdadera” función de densidad de la distribución del estadístico se ha introducido en R siguiendo el siguiente esquema:

    #############################################################

    ### DETERMINACIÓN TEÓRICA DE LA VERDADERA FUNCIÓN DE DENSIDAD

    ### DE LA DISTRIBUCIÓN DEL ESTADÍSTICO

    #############################################################

    # Densidad del estadístico beta.est

    # f(x)=(n/b^n)*x^(n-1)

    ###__________________________________________________________

    ###               

    ###                     CASO 1: b=5 , n=10

    ###__________________________________________________________

    beta <- 5; n <- 10

    rango.x <- seq(from=0, to=beta, by=0.1)

    max <- (n/beta^n)*rango.x^(n-1)

    plot(rango.x, max, type="l", col=1, ylim=c(0,1.6))

    #

    ###__________________________________________________________

    ###               

    ###                     CASO 2: b=5 , n=50

    ###__________________________________________________________

     

    beta <- 5; n <- 50

    rango.x <- seq(from=0, to=beta, by=0.1)

    max <- (n/beta^n)*rango.x^(n-1)

    plot(rango.x, max, type="l", col=1, ylim=c(0,1.6))

    ###__________________________________________________________

    ###               

    ###                     CASO 3: b=10 , n=10

    ###__________________________________________________________

    beta <- 10; n <- 10

    rango.x <- seq(from=0, to=beta, by=0.1)

    max <- (n/beta^n)*rango.x^(n-1)

    plot(rango.x, max, type="l", col=1, ylim=c(0,1.6))

    #

    ###__________________________________________________________

    ###                

    ###                     CASO 4: b=10 , n=50

    ###__________________________________________________________

    beta <- 10; n <- 50

    rango.x <- seq(from=0, to=beta, by=0.1)

    max <- (n/beta^n)*rango.x^(n-1)

    plot(rango.x, max, type="l", col=1, ylim=c(0,1.6))

    De donde se han obtenido las siguientes gráficas correspondientes a cada uno de los casos que se prepresentan a continuación:

     

    CASO 1: b=5 , n=10

    CASO 2: b=5 , n=50

    CASO 3: b=10 , n=10

    CASO 4: b=10 , n=50

     

     

     

     

    De la observación de las gráficas se deduce la mayor influencia del tamaño de muestra n que del parámetro b en la forma de la función de densidad cuya pendiente aumenta cuando lo hace n

b)      Completado el paso anterior intentaremos determinar  la distribucióndel estadístico  mediante simulación, aplicando el principio de Montecarlo:

A tal fin realizaremos la simulación para los mismos valores del parámetro y tamaño muestral anteiores y representaremos a continuación la estimación de la función de densidad mediante la función “density” de R para compararla con la anterior obtenida.

Las condiciones con las que efectuamos el muestreo con las siguientes:

Variable X ~ U  (0, b)

Generamos una muestra de 10.000 valores del estadístico a estimar con el siguiente esquema habitual:

X1  = ( x1,  x2,  ...  x10 )            Þ     max  (X1)   =   1

X2  = ( x1,  x2,  ...  x10 )            Þ     max  (X1)   =   1

X10.000  = ( x1,  x2,  ...  x10 )        Þ    max  (X1)   =   1

###############################################################

### DETERMINACIÓN DE LA DISTRIBUCIÓN DEL ESTADÍSTICO b ESTIMADA

### APLICANDO EL PRINCIPIO DE MONTECARLO

###############################################################

Comenzaremos dividiendo la pantalla gráfica con el fin de obtener en una sola gráfica los cuatro casos. Para lo cual hacemos:

###_________________________________________________________

###               

###         División de la pantalla de gráficos

###__________________________________________________________

m <- matrix(1:4, 2,2)

layout (m, widths=c(1,1),

          heights= c(1,1))

layout.show(4)

###__________________________________________________________               

###

###                     CASO 1: b=5 , n=10

###__________________________________________________________

beta <- 5; n <- 10

# generación de 10000 muestras aleatorias con distribución uniforme

M <- matrix(data=runif(n*10000, 0, beta), nrow=10000, ncol=10, byrow=T)

#obtención de los valores máximos de cada muestra

betas.est <- apply(M, 1, max)

fdens.beta <- density(betas.est, from=0, to=beta)

plot(fdens.beta, type="l", col=1, ylim=c(0,1.6))

###__________________________________________________________

###               

###                     CASO 2: b=5 , n=50

###__________________________________________________________

beta <- 5; n <- 50

#generación de 10000 muestras aleatorias con distribución uniforme

M <- matrix(data=runif(n*10000, 0, beta), nrow=10000, ncol=10, byrow=T)

#obtención de los valores máximos de cada muestra

betas.est <- apply(M, 1, max)

fdens.beta <- density(betas.est, from=0, to=beta)

plot(fdens.beta, type="l", col=1, ylim=c(0,1.6))

###__________________________________________________________

###               

###                     CASO 3: b=10 , n=10

###__________________________________________________________

beta <- 10; n <- 10

#generación de 10000 muestras aleatorias con distribución uniforme

M <- matrix(data=runif(n*10000, 0, beta), nrow=10000, ncol=10, byrow=T)

#obtención de los valores máximos de cada muestra

betas.est <- apply(M, 1, max)

fdens.beta <- density(betas.est, from=0, to=beta)

plot(fdens.beta, type="l", col=1, ylim=c(0,1.6))

###__________________________________________________________

###               

###                     CASO 4: b=10 , n=50

###__________________________________________________________

beta <- 10; n <- 50

#generación de 10000 muestras aleatorias con distribución uniforme

M <- matrix(data=runif(n*10000, 0, beta), nrow=10000, ncol=10, byrow=T)

#obtención de los valores máximos de cada muestra

betas.est <- apply(M, 1, max)

fdens.beta <- density(betas.est, from=0, to=beta)

plot(fdens.beta, type="l", col=1, ylim=c(0,1.6))

De cada uno de los cuatro casos obtenemos graficada la función de densidad, que se muestra a continuación donde se observa una cierta similitud en las mismas


1

3

2

4

Gráfica de los casos

1: b=5  , n=10     3: b=10 , n=10

2: b=5  , n=50     4: b=10 , n=50

Así mismo resulta una similitud entre la distribución “verdadera” del estadístico obtenida mediante gráfico de la función de densidad teórica para cada uno de los cuatro casos con la obtenida por simulación aplicando el principio de Montecarlo, como puede observarse comparándolas

2.      Generamos una muestra aleatoria simple para cada una de las combinaciones del valor de parámetro y tamaño muestral anterior y la tomamos como “datos”.

Este paso, para mayor comodidad lo realizamos dentro de cada uno de las primeras líneas de programa en cada uno de los cuatro casos casos de bootstrap paramétrico y no paramétrico que siguen a continuación

3.      Estimaremos ahora la distribución muestral mediante el método bootstrap para cada una de las cuatro muestras, empezando por bootstrap no paramétrico.

a)bootstrap no paramétrico

Generamos 10.000 remuestras con reemplazamiento a partir de cada muestra concreta para, a continuación, obtener los valores máximos de cada remuestra y la función de densidad buscada

################################################################

### DETERMINACIÓN DE LA DISTRIBUCIÓN DEL ESTADÍSTICO b ESTIMADA

### APLICANDO BOOTSTRAP NO PARAMÉTRICO

##############################################################

###______________________________________________________________

###               

###                     CASO 1: b=5 , n=10

###___________________________________________________________

# Generación de una muestra aleatoria (nuestra muestra)

x<- runif(n, 0, beta)

#generación de 10000 remuestras de x con reemplazamiento de elementos

muestras.boot <- matrix(sample(x, replace=T, size=n*10000), nrow=10000, ncol=10, byrow=T)

#Obtención de los valores máximos de cada remuestra

betas.est <- apply(muestras.boot, 1, max)

fdens.boot <- density(betas.est, from=0, to=max(betas.est))

plot(fdens.boot, type="l", col=1, ylim=c(0,1.6))

###__________________________________________________________

###               

###                     CASO 2: b=5 , n=50

###__________________________________________________________

# Generación de una muestra aleatoria (nuestra muestra)

x<- runif(n, 0, beta)

#generación de 10000 remuestras de x con reemplazamiento de elementos

muestras.boot <- matrix(sample(x, replace=T, size=n*10000), nrow=10000, ncol=10, byrow=T)

#Obtención de los valores máximos de cada remuestra

betas.est <- apply(muestras.boot, 1, max)

fdens.boot <- density(betas.est, from=0, to=max(betas.est))

plot(fdens.boot, type="l", col=1, ylim=c(0,1.6))

###__________________________________________________________

###               

###                     CASO 3: b=10 , n=10

###__________________________________________________________

# Bootstrap no paramétrico

# Generación de una muestra aleatoria (nuestra muestra)

x<- runif(n, 0, beta)

#generación de 10000 remuestras de x con reemplazamiento de elementos

muestras.boot <- matrix(sample(x, replace=T, size=n*10000), nrow=10000, ncol=10, byrow=T)

#Obtención de los valores máximos de cada remuestra

betas.est <- apply(muestras.boot, 1, max)

fdens.boot <- density(betas.est, from=0, to=max(betas.est))

plot(fdens.boot, type="l", col=1, ylim=c(0,1.6))

###__________________________________________________________

###      &