Introducción a los métodos de Monte Carlo

Introducción

En los procedimientos de análisis de datos aparecen dos grandes clases de problemas numéricos: optimización e integración. No siempre es posible calcular de forma analítica los estimadores asociados a un modelo dado, así que a menudo recurrimos a soluciones numéricas. Una forma de evitar ese problema es usar simulación. La estimación Monte Carlo consiste en simular extracciones hipotéticas de una distribución de probabilidad para calcular cantidades relevantes de dicha distribución.

La idea básica de Monte Carlo consiste en expresar la integral como un valor esperado respecto a alguna distribución de probabilidad y aproximarla usando el estimador por el método de los momentos ($E[g(X)] \\approx \\overline{g(X)} = \\dfrac{1}{n}\\sum g(X_{i})$).

Si tenemos una función continua $g(\\theta)$ y queremos integrarla en el intervalo (a,b), podemos reescribir la integral como un valor esperado de una distribución uniforme $U \\sim U[a,b]$:

Usando el estimador por momentos, la aproximación de la integral es:

Donde los

son valores simulados de una distribución uniforme.

Ejemplo 1: aproximación de una integral exponencial

1) Dada una función $f(x)= e^{x}$, la integral en el intervalo [3,5] es:

La aproximación Monte Carlo de la integral es:

#       Declaración de la función deseada\nf = function(x){return(exp(x))}\n#       Declaración de la función de error absoluto\nerror = function(x,y){return(abs(x-y))}\n#       Valor real de la integral\nans = exp(5)-exp(3)\n\nset.seed(6971)\n#       número de iteraciones\nn = 10^2\n#       datos uniformes simulados\nx= runif(n,3,5)\n#       Aproximación Monte Carlo\nMCa= (5-3)*mean(f(x))\n#       Error de la aproximación\ne = error(ans,MCa)\n\nrest =  data.frame(n = n,MCapprox = MCa,error = e)\nset.seed(6971)\nfor(k in 3:6){\n  n = 10^k\n  x = runif(n,3,5)\n  mca = (5-3)*mean(f(x))\n  rest = rbind(rest,c(n,mca,error(ans,mca) ) )\n}\n\nkable(rest,digits = 5,align = 'c',caption = \"Resultados de la aproximación Monte Carlo de la integral\",\n      col.names =c(\"Número de simulaciones\",\"Aproximación Monte Carlo\",\"Error de la aproximación\"))

Aproximación Monte Carlo generalizada

En el caso general, la aproximación de la integral para una distribución f dada es:

Un algoritmo para construir $\\widehat{I}$ puede describirse con los siguientes pasos:

1) Genera de una distribución f

2) Calcula:

3) Obtén la media muestral: $$\\overline{I} = \\dfrac{1}{n}\\sum_{k=1}^{n}\\dfrac{g(\\theta_k)}{f(\\theta_k)}$$

En el siguiente bloque se presenta la función simple de aproximación Monte Carlo para mostrar cómo funciona el algoritmo, donde a y b son los parámetros de la densidad uniforme, n es el número de simulaciones deseadas y f es la función que queremos integrar.

# Función simple de Monte Carlo\nMCaf = function(n,a,b,f){\nx = runif(n,a,b)\nMCa = (b-a)*mean(f(x))\nreturn(MCa)\n}

Métodos de Monte Carlo en el análisis de datos bayesiano

La idea principal del análisis bayesiano de datos es ajustar un modelo (como una regresión o un modelo de series temporales) usando un enfoque de inferencia bayesiana. Suponemos que nuestros parámetros de interés tienen una distribución teórica; esta distribución (posterior) se actualiza usando la distribución de los datos observados (verosimilitud) y la información previa o externa sobre nuestros parámetros (distribución previa) mediante el teorema de Bayes.

$$p(\\theta / X) \\text{ } \\alpha\\text{ } p(X/\\theta)p(\\theta)$$ Donde:

$p(\\theta/X)$ es la distribución posterior del parámetro.

$p(X/\\theta)$ es la distribución muestral de los datos observados (verosimilitud).

$p(\\theta)$ es la distribución previa del parámetro.

El principal problema en el enfoque bayesiano es estimar la distribución posterior. Los métodos de Monte Carlo por cadenas de Markov (mcmc) generan una muestra de la posterior y aproximan valores esperados, probabilidades o cuantiles usando métodos de Monte Carlo.

En las dos secciones siguientes, mostramos dos ejemplos para aproximar probabilidades y cuantiles de una distribución teórica. El procedimiento presentado arriba es la metodología habitual en un enfoque bayesiano.

Ejemplo 2: aproximación de probabilidades de una distribución gamma

Supongamos que queremos calcular la probabilidad de que una variable aleatoria $\\theta$ esté entre cero y 5, $P(0 < \\theta < 5)$, donde $\\theta$ tiene distribución gamma con parámetros a = 2 y b = 1/3 ($\\theta \\sim Gamma(a = 2,b = 1/3)$). Entonces, la probabilidad es:

Donde $I_{[0,5]}(\\theta) = 1$ si $\\theta$ pertenece al intervalo [0,5]. La idea de la aproximación Monte Carlo es contar cuántas observaciones pertenecen al intervalo [0,5] y dividirlo entre el total de datos simulados.

set.seed(6972)\n#       número de iteraciones\nn = 10^2\n#       datos simulados de gamma\nx= rgamma(n,shape = 2,1/3)\n#       Aproximación Monte Carlo\nMCa= mean(x <= 5)\n#       Error de la aproximación\ne = error(pgamma(5,2,1/3),MCa)\n\nrest =  data.frame(n = n,MCapprox = MCa,error = e)\n\nfor(k in 3:6){\n  n = 10^k\n  x= rgamma(n,shape = 2,1/3)\n  mca= mean(x <= 5)\n rest = rbind(rest,c(n,mca,error(pgamma(5,2,1/3),mca) ) )\n}\n\nkable(rest,digits = 5,align = 'c',caption = \"Resultados de la aproximación Monte Carlo de la probabilidad\",\n      col.names =c(\"Número de simulaciones\",\"Aproximación Monte Carlo\",\"Error de la aproximación\"))

Ejemplo 3: aproximación de un cuantil de una distribución normal

Supongamos que queremos calcular el cuantil 0,95 de una variable aleatoria $\\theta$ con distribución normal con parámetros $\\mu = 20$ y $\\sigma = 3$ ($\\theta \\sim normal(\\mu = 20,\\sigma^{2} = 9)$). Entonces, el cuantil 0,95 es:

La idea es encontrar el mayor valor muestral cuya probabilidad sea igual o menor que 0,95. La aproximación Monte Carlo del cuantil se estima con la función quantile() aplicada a los datos simulados.

set.seed(6973)\n#       número de iteraciones\nn = 10^2\n#       datos simulados normales\n  x= rnorm(n,20,3)\n#       Aproximación Monte Carlo\nMCa= quantile(x,0.95)\n#       Error de la aproximación\ne = error(qnorm(0.95,20,3),MCa)\n\nrest =  data.frame(n = n,MCapprox = MCa,error = e)\n\nfor(k in 3:6){\n  n = 10^k\n  x= rnorm(n,20,3)\n  mca= quantile(x,0.95)\n rest = rbind(rest,c(n,mca,error(qnorm(0.95,20,3),mca) ) )\n}\n\nkable(rest,digits = 5,align = 'c',caption = \"Resultados de la aproximación Monte Carlo del cuantil\",\n      col.names =c(\"Número de simulaciones\",\"Aproximación Monte Carlo\",\"Error de la aproximación\"),row.names = FALSE)

Debate y conclusiones

Los métodos de aproximación de Monte Carlo ofrecen una herramienta alternativa para aproximar integrales y son fundamentales en el enfoque de inferencia bayesiana, especialmente cuando trabajamos con modelos sofisticados y complejos. Como se aprecia en los tres ejemplos, los métodos de Monte Carlo proporcionan una excelente aproximación, pero requieren un número muy elevado de simulaciones para lograr un error cercano a cero.

Referencias

1) Introducing Monte Carlo methods with R, Springer 2004, Christian P. Robert and George Casella.

2) Handbook of Markov Chain Monte Carlo, Chapman and Hall, Steve Brooks, Andrew Gelman, Galin L. Jones, and Xiao-Li Meng.

3) Introduction to mathematical Statistics, Pearson, Robert V. Hogg, Joseph W. Mckean, and Allen T. Craig.

4) Statistical Inference An Integrated Approach, Chapman and Hall, Helio S. Migon, Dani Gamerman, Francisco Louzada.