Cuando aplicamos modelos probabilísticos a los datos normalmente es necesario integrar distribuciones de probabilidad complejas. El ejemplo más sencillo lo tenemos en el cálculo del valor esperado. En ocasiones estas integrales no se pueden calcular analíticamente y se usan métodos numéricos de muestreo para aproximar la integral.
El conocido método de Monte Carlo genera muestras aleatorias, independientes e idénticamente distribuidas para así aproximar los cálculos. Sin embargo, para distribuciones complejas no es posible generar estas muestras directamente y aquí es donde surge el método de Monte Carlo basado en cadenas de Markov.
Este método genera una cadena de Markov cuya distribución de equilibrio es la distribución objetivo. Si π(x) es la distribución de probabilidad deseada, el método genera una cadena de Markov con una probabilidad de transición p(x, y) de forma que π(x) sea la distribución estacionaria de la cadena. Al ser una cadena de Markov las muestras generadas ya no son independientes y están correlacionadas.
Uno de los algoritmos más utilizados para generar estas muestras es el Metropolis-Hastings que realiza los siguientes pasos:
- Partimos de una función f(x) que es proporcional a la distribución deseada P(X).
- Elige un primer valor arbitrario y(0) para la primera muestra.
- Repite los siguientes pasos para cada iteración:
- Genera el siguiente valor x usando la distribución g(x|y). Esta distribución podría ser una gaussiana o se podría realizar un incremento uniformemente distribuido.
- Calcula el ratio de aceptación r=f(x)/f(y) y si es mayor que un valor umbral (p.e. extraído de una distribución uniforme) acepta el valor propuesto x y si no repite el valor.
Debido a que el algoritmo acepta los movimientos a puntos más probables según la distribución P(x), la muestra se concentrará en regiones con alta densidad de probabilidad P(x) que es lo que deseamos.
A continuación vamos a ver una simulación del método usando Python. La distribución objetivo es la suma de dos distribuciones bidimensionales normales centradas en (-5, -5) y (5, 5) y con desviaciones típicas (2, 2) y (0.5, 2) respectivamente.
Cogiendo como punto inicial el (0, 0), un movimiento máximo en cada iteración de 3.5 y realizando 5000 iteraciones vemos como se aproxima bien la distribución:
En cambio, si realizamos la misma simulación pero empezando en el punto (5, 5) y bajando el movimiento máximo a 1 vemos como la secuencia no consigue salir de la región de la distribución gaussiana centrada en (5, 5).
Por lo tanto es muy importante controlar el punto de partida y el movimiento máximo para que la secuencia recorra todas las regiones con probabilidad relevante.
Deja una respuesta