Estadísticas bayesianas, probablemente haya encontrado MCMC. Mientras el resto del mundo está obsesionado con las últimas novedades del LLM, Markov Chain Monte Carlo sigue siendo el silencioso caballo de batalla de las finanzas cuantitativas de alto nivel y la gestión de riesgos. Es la herramienta preferida cuando “adivinar” no es suficiente y es necesario mapear rigurosamente la incertidumbre.
A pesar del intimidante acrónimo, Markov Chain Monte Carlo es una combinación de dos conceptos sencillos:
Una Cadena de Markov es un proceso estocástico en el que el siguiente estado del sistema depende completamente de su estado actual y no de las secuencias de eventos que lo precedieron. Esta propiedad suele denominarse falta de memoria. Un método de Monte Carlo simplemente se refiere a cualquier algoritmo que se basa en un muestreo aleatorio repetido para obtener resultados numéricos.
En esta serie, presentaremos los algoritmos centrales utilizados en los marcos MCMC. Nos centramos principalmente en los utilizados para los métodos bayesianos.
Comenzamos con Metropolis-Hastings: el algoritmo fundamental que permitió los primeros avances en este campo. Pero antes de profundizar en la mecánica, analicemos el problema que los métodos MCMC ayudan a resolver.
El problema
Supongamos que queremos poder muestrear variables de una distribución de probabilidad cuya fórmula de densidad conocemos. En este ejemplo utilizamos la distribución normal estándar. Llamemos a una función que pueda tomar muestras de su norma.
Para que rnorm se considere útil, debe generar valores (x) a largo plazo que coincidan con las probabilidades de nuestra distribución objetivo. En otras palabras, si dejáramos que rnorm se ejecutara (100,000) veces y si recopiláramos estos valores y los graficaramos por la frecuencia con la que aparecieron (un histograma), la forma se parecería a la distribución normal estándar.
¿Cómo podemos lograr esto?
Comenzamos con la fórmula para la densidad no normalizada de la distribución normal:
[p(x) = e^{-frac{x^2}{2}}]
Esta función devuelve una densidad para un (x) dado en lugar de una probabilidad. Para obtener una probabilidad, necesitamos normalizar nuestra función de densidad mediante una constante, de modo que el área total bajo la curva se integre (suma) a (1).
Para encontrar esta constante necesitamos integrar la función de densidad en todos los valores posibles de (x).
[C = int^infty_{-infty}e^{-frac{x^2}{2}},dx]
No existe una solución de forma cerrada para la integral indefinida de (e^{-x^2}). Sin embargo, los matemáticos resolvieron la integral definida (de (-infty) a (infty)) moviéndose a coordenadas polares (porque aparentemente, convertir un problema (1D) en uno (2D) lo hace más fácil de resolver), y al darse cuenta de que el área total es (sqrt{2pi}).
Por lo tanto, para que el área bajo la curva sume (1), la constante debe ser la inversa:
[C = frac{1}{sqrt{2pi}}]
De aquí proviene la conocida constante de normalización (C) para la distribución normal.
Muy bien, tenemos la constante dada por los matemáticos que hace que nuestra distribución sea una distribución de probabilidad válida. Pero aún necesitamos poder tomar muestras de él.
Dado que nuestra escala es continua e infinita, la probabilidad de obtener exactamente un número específico (por ejemplo, (x = 1,2345…) con precisión infinita) es en realidad cero. Esto se debe a que un único punto no tiene ancho y, por lo tanto, no contiene ningún "área" debajo de la curva.
En cambio, debemos hablar en términos de rangos, es decir, ¿cuál es la probabilidad de obtener un valor (x) que se encuentre entre (a) y (b) ((a < x < b)), donde (a) y (b) son valores fijos?
En otras palabras, necesitamos encontrar el área bajo la curva entre (a) y (b). Y como habrás adivinado correctamente, para calcular esta área normalmente necesitamos integrar nuestra fórmula de densidad normalizada: la función integrada resultante se conoce como función de distribución acumulativa ((CDF)).
Como no podemos resolver la integral, no podemos derivar el (CDF) y por eso estamos estancados nuevamente. Los matemáticos inteligentes también resolvieron esto y pudieron usar la trigonometría (específicamente la transformada de Box-Muller) para convertir variables aleatorias uniformes en variables aleatorias normales.
Pero aquí está el truco: en el mundo real de la estadística bayesiana y el aprendizaje automático, tratamos con distribuciones multidimensionales complejas donde:
No podemos resolver la integral analíticamente. Por lo tanto, no podemos encontrar la constante de normalización (C). Finalmente, sin la integral, no podemos calcular la CDF, por lo que el muestreo estándar falla.
Estamos atrapados con una fórmula no normalizada y sin forma de calcular el área total. Aquí es donde entra en juego MCMC. Los métodos MCMC nos permiten tomar muestras de estas distribuciones sin necesidad de resolver esa integral imposible.
Introducción
Un proceso de Markov se define únicamente por sus probabilidades de transición (P(xrightarrow x')).
Por ejemplo, en un sistema con (4) estados:
[P(xrightarrow x') = begin{bmatrix} 0,5 y 0,3 y 0,05 y 0,15 \ 0,2 y 0,4 y 0,1 y 0,3 \ 0,4 y 0,4 y 0 y 0,2 \ 0,1 y 0,8 y 0,05 y 0,05 end{bmatrix}]
La probabilidad de pasar de cualquier estado (x) a (x') está dada por la entrada (i rightarrow j) en la matriz.
Mire la tercera fila, por ejemplo: ([0.4,0.4,0,0.2]).
Nos dice que si el sistema se encuentra actualmente en el Estado (3), tiene una probabilidad (40%) de pasar al Estado (1), una probabilidad (40%) de pasar al Estado (2), una probabilidad (0%) de permanecer en el Estado (3), y una probabilidad (20%) de pasar al Estado (4).
La matriz ha mapeado todos los caminos posibles con sus correspondientes probabilidades. Observe que cada fila (i) suma (1) de modo que nuestra matriz de transición representa probabilidades válidas.
Un proceso de Markov también requiere una distribución de estado inicial (pi_0) (¿comenzamos en el estado (1) con (100%) probabilidad o en cualquiera de los (4) estados con (25%) probabilidad cada uno?).
Por ejemplo, esto podría verse así:
[pi_0 = begin{bmatrix} 0,4 y 0,15 y 0,25 y 0,2 end{bmatrix}]
Esto simplemente significa que la probabilidad de comenzar desde el estado (1) es (0.4), el estado (2) es (0.15) y así sucesivamente.
Para encontrar la distribución de probabilidad de dónde estará el sistema después del primer paso (t_0 + 1), multiplicamos la distribución inicial por las probabilidades de transición:
[ pi_1 = pi_0P]
La multiplicación de matrices efectivamente nos da la probabilidad de todas las rutas que podemos tomar para llegar a un estado (j) sumando todas las probabilidades individuales de alcanzar (j) desde diferentes estados iniciales (i).
¿Por qué esto funciona?
Al utilizar la multiplicación de matrices, exploramos todos los caminos posibles hacia un destino y sumamos su probabilidad.
Observe que la operación también preserva maravillosamente el requisito de que la suma de las probabilidades de estar en un estado siempre será igual a (1).
Distribución estacionaria
Un proceso de Markov construido adecuadamente alcanza un estado de equilibrio cuando el número de pasos (t) se acerca al infinito:
[pi^* P = pi^*]
Este estado se conoce como equilibrio global.
(pi^*) se conoce como distribución estacionaria y representa un momento en el que la distribución de probabilidad después de una transición ((pi^*P)) es idéntica a la distribución de probabilidad antes de la transición ((pi^*)).
La existencia de tal estado resulta ser la base de todo método MCMC.
Al muestrear una distribución objetivo mediante un proceso estocástico, no nos preguntamos "¿Adónde vamos ahora?" sino más bien “¿Dónde acabaremos finalmente?”. Para responder a eso, necesitamos introducir previsibilidad a largo plazo en el sistema.
Esto garantiza que existe un estado teórico (t) en el que las probabilidades se "estabilizan" en lugar de moverse aleatoriamente por toda la eternidad. El punto en el que se “asientan” es el punto en el que esperamos poder comenzar a tomar muestras de nuestra distribución objetivo.
Por lo tanto, para poder estimar efectivamente una distribución de probabilidad usando un proceso de Markov necesitamos asegurarnos de que:
existe una distribución estacionaria. esta distribución estacionaria es única; de lo contrario, podríamos tener múltiples estados de equilibrio en un espacio muy alejado de nuestra distribución objetivo.
Las restricciones matemáticas impuestas por el algoritmo hacen que un proceso de Markov satisfaga estas condiciones, lo cual es fundamental para todos los métodos MCMC. La forma en que se logra esto puede variar.
Unicidad
En general, para garantizar la unicidad de la distribución estacionaria, necesitamos satisfacer tres condiciones. Expande la sección a continuación para verlos:
La Santísima Trinidad Irreducible: Un sistema es irreducible si en cualquier estado (x) existe una probabilidad distinta de cero de que se visite cualquier punto (x') en el espacio muestral. En pocas palabras, eventualmente puedes pasar de cualquier estado A a cualquier estado B. Aperiódico: El sistema no debe regresar a un estado particular en intervalos fijos. Una condición suficiente para la aperiodicidad es que exista al menos un estado donde la probabilidad de permanecer sea distinta de cero. Recurrente positivo: Un estado (x) es recurrente positivo si, a partir de ese estado, se garantiza que el sistema regresará a él y el número promedio de pasos que toma para regresar es finito. Esto lo aseguramos modelando un objetivo que tiene una integral finita y una distribución de probabilidad adecuada (el área bajo la curva debe sumar (1)).
Cualquier sistema que cumpla estas condiciones se conoce como sistema ergódico. Las tablas al final del artículo muestran cómo el algoritmo MH se ocupa de garantizar la ergodicidad y, por tanto, la unicidad.
Metrópolis-Hastings
El enfoque que adopta el algoritmo MH es comenzar con la definición de equilibrio detallado, una condición suficiente pero no necesaria para el equilibrio global. En pocas palabras, si nuestro algoritmo satisface el equilibrio detallado, garantizaremos que nuestra simulación tenga una distribución estacionaria.
Derivación
La definición de saldo detallado es:
[pi(x) P(x'|x) = pi(x') P(x|x') ,]
esto significa que el flujo de probabilidad de ir de (x) a (x') es el mismo que el flujo de probabilidad de ir de (x') a (x).
La idea es encontrar la distribución estacionaria construyendo iterativamente la matriz de transición, (P(x',x)) estableciendo (pi) como la distribución objetivo (P(x)) de la que queremos tomar la muestra.
Para implementar esto, descomponemos la probabilidad de transición (P(x'|x)) en dos pasos separados:
Propuesta ((g)): La probabilidad de proponer un movimiento a (x') dado que estamos en (x). Aceptación ((A)): La función de aceptación nos da la probabilidad de aceptar la propuesta.
De este modo,
[P(x'|x) = g(x'|x) a(x',x)]
La corrección de Hastings
Sustituyendo estos valores nuevamente en la ecuación anterior nos da:
[frac{pi(x)}{pi(x')} = frac{g(x|x') a(x,x')}{g(x'|x) a(x',x)}]
y finalmente reordenando obtenemos una expresión para nuestra aceptación como proporción:
[frac{a(x',x)}{a(x,x')} = frac{g(x|x')pi(x')}{g(x'|x)pi(x)}]
Esta relación representa la probabilidad de que aceptemos un movimiento hacia (x') en comparación con un regreso a (x).
El término (frac{g(x|x')pi(x')}{g(x'|x)pi(x)}) se conoce como corrección de Hastings.
Nota importante
Debido a que a menudo elegimos una distribución simétrica para la propuesta, la probabilidad de saltar desde (x rightarrow x') es la misma que saltar desde (x' rightarrow x). Por lo tanto, los términos de la propuesta se cancelan entre sí dejando solo la relación de las densidades objetivo.
Este caso especial donde la propuesta es simétrica y los términos (g) desaparecen se conoce históricamente como Algoritmo de Metropolis (1953). La versión más general que permite propuestas asimétricas (que requieren la relación (g), conocida como corrección de Hastings) es el algoritmo Metropolis-Hastings (1970).
El gran avance
Recuerde el problema original: no podemos calcular (pi(x)) porque no conocemos la constante de normalización (C) (la integral).
Sin embargo, observe de cerca la relación (frac{pi(x')}{pi(x)}). Si expandimos (pi(x)) a su densidad no normalizada (f(x)) y la constante (C):
[frac{pi(x')}{pi(x)} = frac{{f(x')} / C}{f(x) / C} = frac{f(x')}{f(x)}]
¡La constante (C) se cancela!
Este es el gran avance. Ahora podemos tomar muestras de una distribución compleja utilizando sólo la densidad no normalizada (que conocemos) y la distribución propuesta (que elegimos).
Todo lo que queda por hacer es encontrar una función de aceptación (A) que satisfaga el equilibrio detallado:
[frac{a(x',x)}{a(x,x')} = R ,]
donde (R) representa (frac{g(x|x')pi(x')}{g(x'|x)pi(x)}).
La aceptación de la metrópoli
La función de aceptación que utiliza el algoritmo es:
[a(x',x) = min(1,R)]
Esto asegura que la probabilidad de aceptación esté siempre entre (0) y (1).
Para ver por qué esta elección satisface el equilibrio detallado, debemos verificar que la ecuación también se cumple para el movimiento inverso. Necesitamos verificar que:
[frac{a(x',x)}{a(x,x')} = R,]
en dos casos:
Caso I: La jugada es ventajosa ((R ge 1))
Dado que (R ge 1), la inversa (frac{1}{R} le 1):
nuestra aceptación directa es (a(x',x) = min(1,R) = 1) nuestra aceptación inversa es (a(x,x') = min(1,frac{1}{R}) = frac{1}{R})
[frac{1}{a(x,x')} = R]
[frac{1}{frac{1}{R}} = R]
Caso II: La jugada no es ventajosa ((R < 1))
Dado que (R < 1), la inversa (frac{1}{R} > 1):
nuestra aceptación directa es (a(x',x) = min(1,R) = R) nuestra aceptación inversa es (a(x,x') = min(1,frac{1}{R}) = 1)
De este modo:
[frac{R}{a(x,x')} = R]
[frac{R}{1} = R,]
y la igualdad se cumple en ambos casos.
Implementación
Implementemos el algoritmo MH en Python en dos distribuciones de destino de ejemplo.
I. Estimación de una distribución gaussiana
Cuando trazamos las muestras en un gráfico frente a una distribución normal verdadera, esto es lo que obtenemos:
Ahora quizás te estés preguntando por qué nos molestamos en ejecutar un método MCMC para algo que podemos hacer usando np.random.normal(n_iterations). ¡Ese es un punto muy válido! De hecho, para un gaussiano unidimensional, la solución de transformada inversa (usando trignometría) es mucho más eficiente y es lo que realmente usa numpy.
¡Pero al menos sabemos que nuestro código funciona! Ahora, intentemos algo más interesante.
II. Estimación de la distribución del 'volcán'
Intentemos tomar muestras de una distribución mucho menos "estándar" que se construye en dos dimensiones, donde la tercera dimensión representa la densidad de la distribución.
Dado que el muestreo se realiza en el espacio (2D) (el algoritmo sólo conoce su ubicación xy, no la 'pendiente' del volcán), obtenemos un bonito anillo alrededor de la boca del volcán.
Resumen de condiciones matemáticas para MCMC
Ahora que hemos visto la implementación básica, aquí hay un resumen rápido de las condiciones matemáticas que requiere un método MCMC para funcionar realmente:
(Que existe un conjunto de probabilidades que, una vez alcanzadas, no cambiarán). Saldo detallado
El algoritmo está diseñado para satisfacer la ecuación de equilibrio detallada. Convergencia
(Garantizando que la cadena eventualmente converja a la distribución estacionaria). Ergodicidad
El sistema debe satisfacer las condiciones de la tabla 2 y ser ergódico. Unicidad de la distribución estacionaria
(Que existe sólo una solución para la ecuación de equilibrio detallada)Ergodicidad
Garantizado si el sistema es ergonómico.
Y así es como el algoritmo MH satisface los requisitos de ergodicidad:
(Capacidad de llegar a cualquier estado desde cualquier otro estado). Función de propuesta
A menudo se satisface utilizando una propuesta (como una gaussiana) que tiene una probabilidad distinta de cero en todas partes. Nota: Si no es posible saltar a algunas regiones, esta condición falla.2. Aperiódico
(El sistema no queda atrapado en un bucle). Paso de rechazo
El “lanzamiento de moneda” nos permite rechazar un movimiento y permanecer en el mismo estado, rompiendo cualquier periodicidad que se haya podido producir.3. Positivo Recurrente
(El tiempo de regreso esperado a cualquier estado es finito). Distribución de probabilidad adecuada
Garantizado por el hecho de que modelamos el objetivo como una distribución adecuada (es decir, integra/suma (1)).
Conclusión
En este artículo hemos visto cómo MCMC ayuda a resolver los dos desafíos principales del muestreo de una distribución dada solo su función de densidad (probabilidad) no normalizada:
El problema de normalización: para que una distribución sea una distribución de probabilidad válida, el área bajo la curva debe sumar (1). Para hacer esto necesitamos calcular el área total bajo la curva y luego dividir nuestros valores no normalizados por esa constante. Calcular el área implica integrar una función compleja y, en el caso de la distribución normal, por ejemplo, no existe una solución de forma cerrada. El problema de inversión: para generar una muestra necesitamos elegir una probabilidad aleatoria y preguntar qué valor de (x) corresponde a esta área. Para ello no sólo tenemos que resolver la integral sino también invertirla. Y como no podemos escribir la integral es imposible resolver su inversa.
Los métodos MCMC, comenzando con Metropolis-Hastings, nos permiten evitar estos problemas matemáticos imposibles mediante el uso de inteligentes paseos aleatorios y ratios de aceptación.
Para obtener una implementación más sólida del algoritmo Metropolis-Hastings y un ejemplo de muestreo utilizando una propuesta asimétrica (utilizando la corrección de Hastings), consulte el código de inicio aquí.
¿Qué sigue?
Hemos muestreado con éxito una distribución (2D) compleja sin siquiera calcular una integral. Sin embargo, si observa el código de Metropolis-Hastings, notará que nuestro paso de propuesta es esencialmente una suposición ciega y aleatoria (np.random.normal).
En dimensiones bajas, las conjeturas funcionan. Pero en los espacios de alta dimensión de los métodos bayesianos modernos, adivinar al azar es como tratar de obtener una mejor tasa de un usurero: casi todas las propuestas que hagas serán rechazadas.
En la Parte II, presentaremos el Hamiltoniano Monte Carlo (HMC), un algoritmo que nos permite explorar eficientemente el espacio de alta dimensión utilizando la geometría de la distribución para guiar nuestros pasos.