Métodos MCMC, un poco de historia.

Cadena de Markov de 3 estados. Cada estado tiene una probabilidad de pasar a alguno de los otros o quedarse en el mismo.

Hoy os traigo unos métodos a los que tengo mucho cariño ya que fue la base de mi proyecto final de carrera, los MCMC. Las siglas MCMC significan Markov Chain Monte Carlo, es decir, métodos Monte Carlo usando cadenas de Markov. Para entender estos métodos lo mejor será explicar por partes qué son las cadenas de Markov y qué son los métodos de Monte Carlo. Empecemos.

Cadenas de Markov

Las cadenas de Markov son un sistema numérico que sufre transiciones de un estado a otro dentro de un espacio de estados. Este proceso estocástico fue definido por Andréi Markov en 1907 [1] y se aplica en física estadística, termodinámica, meteorología (si hay nubes puede llover, nevar, quedarse como está o despejarse, desde soleado tenemos la posibilidad de que se nuble o que se quede como está, pero no va a empezar a nevar de golpe, etc.), epidemiología (usando el modelo de Galton-Watson podemos saber cómo se desarrollará una epidemia), economía, finanzas, genética, y un largo etcétera. La figura que tenéis a la derecha es una cadena de Markov simple de 3 estados definidos. Por ejemplo, partiendo del estado E2 podemos ir al estado E1 con una probabilidad de 0.4, ir al E3 con una probabilidad de 0.2 y quedarse en el estado E2 con una probabilidad de 0.4. El siguiente estado depende del estado anterior pero no de ninguno previo, a esta cualidad la conocemos como propiedad de Markov y se define matemáticamente de la siguiente forma: 

El origen del método

Stanisław Marcin Ulam desarrolló herramientas matemáticas para teoría de conjuntos, topología algebraica, etc.

Ya hemos definido las cadenas de Markov, ahora toca explicar los métodos de Monte Carlo. Los métodos se Monte Carlo son muy conocidos entre los físicos, ya que estos se utilizan para calcular, aproximar y simular expresiones o sistemas matemáticos complejos y difíciles de evaluar mediante un sistema no determinista o estadístico numérico. El nombre de Monte Carlo es en honor al casino de Monte Carlo en Mónaco ya que aquí, como en todos los casinos, se juega a juegos de azar.

El iniciador de este método fue Stanislaw Marcin Ulam, que tras haber pasado varios días jugando al solitario debido a una enfermedad, observó que era mucho más simple tener una idea del resultado general del juego haciendo pruebas múltiples y contando las proporciones de los resultados, que calcular todas las posibilidades. Esta idea la consiguió extrapolar a su trabajo en el Proyecto Manhattan en el que estaba teniendo problemas para calcular las ecuaciones integro-diferenciales que regían la difusión de neutrones. El Proyecto Manhattan se desarrolló durante la Segunda Guerra Mundial para el Gobierno de los Estados Unidos, con la colaboración de Reino Unido y de Canadá. El objetivo era desarrollar la primera bomba atómica antes que Alemania (Proyecto Uranio) y la URSS (Proyecto Borodino).

John von Neumann realizó contribuciones en física cuántica, análisis funcional, economía, análisis numérico, etc.

En las primeras bombas atómicas desarrolladas no se tuvo en cuenta la idea de Ulam. Años más tarde, una vez finalizada la guerra, hubo tiempo de aplicar este método recién salido del cascarón en las investigaciones de física nuclear y física cuántica. Una vez que Ulam convenció a su colega de trabajo John von Neumann del potencial de este nuevo método, ambos matemáticos trabajaron conjuntamente para el desarrollo de éste. El verdadero potencial del método lo obtuvieron Ulam, Nicholas Constantine Metropolis (matemático, físico y computador científico griego reclutado en el Proyecto Manhattan por Oppenheimer) y Enrico Fermi (físico italiano que realizó contribuciones en mecánica estadística, física cuántica , nuclear y de partículas) al obtener los valores estimados de la ecuación de Schrödinger – que describe la evolución temporal de una partícula masiva no relativista – para la captura de neutrones a nivel nuclear  [2].

Aunque la investigación con este método comenzó en los proyectos secretos que llevó a cabo Estados Unidos, la primera publicación donde se citó este método data de 1949 en un artículo de Metropolis y Ulam. La primera simulación de Monte Carlo fue llevada a cabo por un equipo encabezado por Metropolis. Esta simulación se realizó en la computadora ENIAC (considerada la primera computadora electrónica digital) en 1948 en la Universidad de Pennsylvania. Está valorado como uno de los 10 algoritmos más importantes desarrollados en el S.XX.

El método MCMC

Con los métodos y algoritmos MCMC se pueden resolver multitud de problemas como de optimización, cálculo de volúmenes complejos, minimización de costes, etc. En este post nos vamos a centrar en generar muestras aleatorias que correspondan a una determinada función de probabilidad.

Función de probabilidad de un dado de 6 caras suponiendo que cada una de las caras tiene la misma probabilidad de salir, es decir, que el dado no esté truncado.

Pero, ¿qué es una función de probabilidad? Por ejemplo, la función de probabilidad de un dado de seis caras daría valores de 1/6 distribuidos en 1, 2, 3, 4, 5 y 6. Si sumamos todas las probabilidades, el resultado es 1. En este caso la función solo toma valores discretos, pero también podemos tener funciones de probabilidad continuas, como por ejemplo la altura de la población. Si calculamos el área de la función de probabilidad de una variable aleatoria continua (integrando la función), lo que obtenemos es también el valor 1.

En azul podemos observar una función objetivo p(x) de la que queremos sacar muestras, en rojo podemos observar la función tentativa de la que sabemos sacar muestras o generar números aleatorios.

Volvamos a lo que nos atañe, imagina que queremos generar muestras aleatorias con una determinada distribución pero esta es muy rara, con muchos picos, valles, etc. y no disponemos de un generador de números aleatorios que se adecue a esa distribución. Existen generadores de gaussianas (los procesos naturales se adaptan a este tipo de distribuciones) y de distribuciones uniformes (el ejemplo del dado o de las ruletas de los casinos). Si la distribución a generar, llamémosla p(x), no la podemos dividir o no sabemos dividirla como suma de estas distribuciones que conocemos, bien cambiando medias, varianzas y pesos, estamos fastidiados. Aquí es donde entran los métodos MCMC: cogemos la distribución p(x) de la que queremos obtener muestras (números aleatorios), y la utilizamos para evaluar muestras que generamos desde un generador π(x) que conocemos bien, por ejemplo un generador de muestras con distribución gaussiana. 

El algoritmo

Os presento a continuación el algoritmo Metropolis-Hastings, que es el más aceptado y conocido dentro de los algoritmos MCMC. La explicación de su funcionamiento es sencilla y la idea es válida para el resto de los MCMC:

1. Generamos una muestra x’ de la densidad π(x_t|x_t-1), siendo x_t-1 la última muestra válida.
2. Calculamos el ratio de aceptación 
3.  Elegimos x_t=x’ con una probabilidad α(x_t-1,x’), o mantenemos x_t=x_t-1 con una probabilidad de 1-α(x_t-1,x’). Volvemos al paso número 1.

Hay más funciones de aceptación que la mostrada anteriormente, pero la óptima es la que os he dispuesto según Peskun[3]. La demostración de que esto funciona y todas las muestras x₀,x₁,x₂,…,xn que obtenemos se distribuyen como p(x) es larga y tediosa. Si queréis verla podéis descargaros mi PFC [4].

En la figura sólo está representada p(x). También se ha forzado a que π(xt|xt-1) = π(xt-1|xt), de esta forma el ratio de aceptación es mucho más sencillo de calcular.

El resultado

Un ejemplo real: en la figura de abajo podéis ver en azul la función de la que queremos generar muestras. En rojo se encuentra la representación gráfica de los valores obtenidos donde la superficie de cada barra es proporcional a la frecuencia de los valores representados.  

Función de densidad de probabilidad en azul. Histrograma de las muestras de MCMC en rojo.

Si a la salida del algoritmo hubiésemos puesto algo que dibujase los valores, obtendríamos algo como la imagen de abajo. Los valores del algoritmo se van a repetir en ocasiones o modificarse muy poco durante pequeños intervalos de tiempo, pero de vez en cuando se producen saltos bruscos que permiten que el método recorra toda la función de probabilidad y no se quede estancada.

En verde se encuentran todas las muestras que ha generado la cadena de Markov. Esta ha necesitado unas 5000 muestras para converger a la densidad de probabilidad tentativa.

Los valores obtenidos no son muestras aleatorias e independientes, así que lo que se hace para obtenerlas se elige de forma aleatoria los valores obtenidos de la ejecución del algoritmo. Por ejemplo, si quieres 100 muestras aleatorias con una determinada función de probabilidad, genera 1000 valores con el método MCMC y luego elige una muestra de cada diez de forma aleatoria.

Como os podéis imaginar, no son métodos muy rápidos que digamos. Algunas simulaciones tardan horas o días en finalizarse. Hay que tener en cuenta que, dependiendo del tipo de problema a resolver, se pueden lanzar varios algoritmos MCMC al mismo tiempo y en paralelo.

Estos métodos son fácilmente modificables y adaptables. Diferentes investigadores han mezclado estos métodos con algoritmos genéticos, de aprendizaje, usando cadenas con diferentes temperaturas (varianza), mezclando varias cadenas, usando una cadena y generando muestras de las que parten otras cadenas y eliminando muestras, y un largo etcétera. 

Estos métodos son muy buenos cuando necesitas un resultado aproximado, y sabes que obtener un resultado exacto no va a cambiar mucho el resultado final. 

Víctor Pascual del Olmo

Referencias:

[1] Basharin, Gely P.; Langville, Amy N.; Naumov, Valeriy A. (2004). «The Life and Work of A. A. Markov» (en inglés). Linear Algebra and its Applications 386: pp. 3-26.

[2] [23] Mackay, D.J.C.: Introduction to Monte Carlo Methods. Springer (1986)

[3] Peskun, P.H.: Optimun Monte Carlo sampling using Markov chains. Biometrika Vol. 60, Issue 3, (May 1973), pp. 607-612

[4] http://e-archivo.uc3m.es/handle/10016/13735

Víctor Pascual