Uso intuitivo de los métodos de Monte Carlo con cadenas de Markov.

¿Es fácil? Lo intenté


Alexey Kuzmin, director de desarrollo de datos y trabajo en DomKlik, profesor de Data Science en Netology, tradujo un artículo de Rahul Agarwal sobre cómo los métodos de Monte Carlo trabajan con las cadenas de Markov para resolver problemas con un gran espacio de estado.

Todos los asociados con Data Science han oído hablar de los métodos de Monte Carlo con cadenas de Markov (MCMC). A veces se aborda el tema al estudiar estadísticas bayesianas, a veces al trabajar con herramientas como Prophet.

Pero MCMC es difícil de entender. Cada vez que leo sobre estos métodos, noto que la esencia de MCMC está oculta en las capas profundas de ruido matemático, y es difícil distinguir detrás de este ruido. Tuve que pasar muchas horas entendiendo este concepto.

En este artículo, está disponible un intento de explicar los métodos de Monte Carlo con cadenas de Markov, para que quede claro para qué se usan. Me centraré en algunas formas más de usar estos métodos en mi próxima publicación.

Entonces comencemos. MCMC consta de dos términos: cadenas Monte Carlo y Markov. Hablemos de cada uno de ellos.

Montecarlo




En los términos más simples , los métodos de Monte Carlo se pueden definir como simulaciones simples.

Monte Carlo Methods obtuvo su nombre del Casino Monte Carlo en Mónaco. En muchos juegos de cartas necesitas saber la probabilidad de ganar el crupier. Algunas veces el cálculo de esta probabilidad puede ser matemáticamente complicado o intratable. Pero siempre podemos ejecutar una simulación por computadora para jugar todo el juego muchas veces y considerar la probabilidad como el número de victorias dividido por el número de juegos jugados.
Esto es todo lo que necesita saber sobre los métodos de Monte Carlo. Sí, es solo una técnica de modelado simple con un nombre elegante.

Cadenas de Markov




Dado que el término MCMC consta de dos partes, aún debe comprender qué son las cadenas de Markov . Pero antes de pasar a las cadenas de Markov, hablemos un poco sobre las propiedades de Markov.

Supongamos que hay un sistema de estados M-posibles, y usted se mueve de un estado a otro. Que nada te confunda todavía. Un ejemplo específico de dicho sistema es el clima, que cambia de caliente a frío a moderado. Otro ejemplo es el mercado de valores, que está saltando de un estado bajista a un estado alcista y estancado.

La propiedad de Markov sugiere que para un proceso dado que se encuentra en el estado X n en un momento particular en el tiempo, la probabilidad X n + 1 = k (donde k es cualquiera de los estados M a los que puede ir el proceso) depende solo de ¿Cuál es esta condición en este momento? Y no sobre cómo llegó a su estado actual.
En términos matemáticos, podemos escribir esto en la forma de la siguiente fórmula:

Para mayor claridad, no le importa la secuencia de condiciones que tomó el mercado para volverse optimista. La probabilidad de que el próximo estado sea "bajista" está determinada solo por el hecho de que el mercado se encuentra actualmente en un estado "alcista". También tiene sentido en la práctica.

Un proceso con una propiedad de Markov se llama proceso de Markov. ¿Por qué es importante la cadena de Markov? Debido a su distribución estacionaria.

¿Qué es la distribución estacionaria?


Trataré de explicar la distribución estacionaria al calcularla para el siguiente ejemplo. Supongamos que tiene un proceso de Markov para el mercado de valores, como se muestra a continuación.

Tiene una matriz de probabilidad de transición que determina la probabilidad de una transición del estado X i a X j .

Matriz de probabilidad de transición, Q

En la matriz de probabilidad de transición dada Q, la probabilidad de que el próximo estado sea "toro", dado el estado actual de "toro" = 0.9; la probabilidad de que el próximo estado sea "bajista" si el estado actual es "toro" = 0.075. Y así sucesivamente.

Bueno, comencemos con un estado particular. Nuestro estado será establecido por el vector [toro, oso, estancamiento]. Si comenzamos con un estado "bajista", el vector será así: [0,1,0]. Podemos calcular la distribución de probabilidad para el siguiente estado multiplicando el vector de estado actual por la matriz de probabilidad de transición.

Tenga en cuenta que las probabilidades suman 1.

La fórmula puede encontrar la siguiente distribución de estados:


Y así sucesivamente. Al final, alcanzará un estado estacionario en el que el estado se estabiliza:


Para la matriz de probabilidad de transición Q descrita anteriormente, la distribución estacionaria s es

Puede obtener una distribución estacionaria con el siguiente código:

Q = np.matrix([[0.9,0.075,0.025],[0.15,0.8,0.05],[0.25,0.25,0.5]]) init_s = np.matrix([[0, 1 , 0]]) epsilon =1 while epsilon>10e-9:    next_s = np.dot(init_s,Q)    epsilon = np.sqrt(np.sum(np.square(next_s - init_s)))    init_s = next_s print(init_s) ------------------------------------------------------------------ matrix([[0.62499998, 0.31250002, 0.0625  ]]) 

También puede comenzar desde cualquier otro estado: logre la misma distribución estacionaria. Cambie el estado inicial en el código si desea asegurarse de esto.

Ahora podemos responder a la pregunta de por qué la distribución estacionaria es tan importante.

La distribución estacionaria es importante porque puede usarse para determinar la probabilidad de que un sistema esté en un cierto estado al azar.

Para nuestro ejemplo, podemos decir que en el 62.5% de los casos el mercado estará en un estado "alcista", el 31.25% en un estado "bajista" y el 6.25% en estancamiento.

Intuitivamente, puedes ver esto como un azar deambulando por la cadena.


Paseo al azar

Está en un cierto punto y elige el siguiente estado, observando la distribución de probabilidad del siguiente estado, teniendo en cuenta el estado actual. Podemos visitar algunos nodos con más frecuencia que otros, según las probabilidades de estos nodos.

Así es como Google resolvió el problema de búsqueda en los albores de Internet. El problema fue ordenar las páginas, dependiendo de su importancia. Google resolvió el problema utilizando el algoritmo Pagerank. El algoritmo Google Pagerank debería considerar el estado como una página, y la probabilidad de una página en una distribución estacionaria como su importancia relativa.

Ahora pasamos directamente a la consideración de los métodos MCMC.

¿Cuáles son los métodos de Monte Carlo con cadenas de Markov (MCMC)


Antes de responder qué es MCMC, déjame hacerte una pregunta. Sabemos acerca de la distribución beta. Conocemos su función de densidad de probabilidad. ¿Pero podemos tomar una muestra de esta distribución? ¿Puedes encontrar una manera de hacer esto?


Piensa ...

MCMC le permite elegir entre cualquier distribución de probabilidad. Esto es especialmente importante cuando necesita hacer una selección desde la distribución posterior.

La figura muestra el teorema de Bayes.

Por ejemplo, necesita hacer una muestra de una distribución posterior. Pero, ¿es fácil calcular el componente posterior junto con la constante de normalización (evidencia)? En la mayoría de los casos, puede encontrarlos en forma de un producto de probabilidad y probabilidad a priori. Pero calcular la constante de normalización (p (D)) no funciona. Por qué Echemos un vistazo más de cerca.

Suponga que H solo toma 3 valores:

p (D) = p (H = H1) .p (D | H = H1) + p (H = H2) .p (D | H = H2) + p (H = H3) .p (D | H = H3)

En este caso, p (D) es fácil de calcular. Pero, ¿qué pasa si el valor de H es continuo? ¿Sería posible calcular esto tan fácilmente, especialmente si H tomara valores infinitos? Para hacer esto, una integral compleja tendría que ser resuelta.

Queremos hacer una selección aleatoria de la distribución posterior, pero también queremos considerar p (D) como una constante.

Wikipedia escribe :

Los métodos de Monte Carlo con cadenas de Markov son una clase de algoritmos para el muestreo de una distribución de probabilidad, basada en la construcción de una cadena de Markov, que como distribución estacionaria tiene la forma deseada. El estado de la cadena después de una serie de pasos se usa como una selección de la distribución deseada. La calidad de muestreo mejora con el número creciente de pasos.

Veamos un ejemplo. Digamos que necesita una muestra de la distribución beta . Su densidad:


donde C es la constante de normalización. En realidad, esta es una función de α y β, pero quiero mostrar que no es necesaria para una muestra de la distribución beta, por lo que la tomaremos como una constante.

El problema de distribución beta es realmente difícil, si no prácticamente insoluble. En realidad, es posible que deba trabajar con funciones de distribución más complejas y, a veces, no conocerá las constantes de normalización.

Los métodos MCMC facilitan la vida al proporcionar algoritmos que podrían crear una cadena de Markov que tiene una distribución beta como una distribución estacionaria, dado que podemos elegir entre una distribución uniforme (que es relativamente simple).

Si comenzamos con un estado aleatorio y pasamos al siguiente estado sobre la base de algún algoritmo varias veces, finalmente crearemos una cadena de Markov con una distribución beta como distribución estacionaria. Y los estados en los que nos encontramos durante mucho tiempo pueden usarse como muestra de la distribución beta.

Uno de estos algoritmos MCMC es el algoritmo Metropolis-Hastings.

Algoritmo Metrópolis-Hastings




Intuición


Entonces, ¿cuál es el propósito?

Intuitivamente, queremos caminar a lo largo de algún pedazo de superficie (nuestra cadena de Markov) de tal manera que la cantidad de tiempo que pasamos en cada ubicación sea proporcional a la altura de la superficie en ese lugar (la densidad de probabilidad deseada de la que queremos hacer una selección).

Entonces, por ejemplo, nos gustaría pasar el doble de tiempo en una colina de 100 metros de altura que en una colina vecina de 50 metros. Es bueno que podamos hacer esto incluso si no conocemos las alturas absolutas de los puntos en la superficie: todo lo que necesita saber son las alturas relativas. Por ejemplo, si la cima de la colina A es dos veces más alta que la cima de la colina B, entonces nos gustaría pasar el doble de tiempo en A que en B.

Existen esquemas más complejos para proponer nuevos lugares y reglas para su adopción, pero la idea principal es la siguiente:

  1. Elija una nueva ubicación "sugerida".
  2. Descubra qué tan alta o baja es esta ubicación en comparación con la actual.
  3. Permanecer en el lugar o mudarse a un nuevo lugar con una probabilidad proporcional a las alturas de los lugares.

El propósito del MCMC es seleccionar de alguna distribución de probabilidad sin tener que saber su altura exacta en ningún punto (no es necesario saber C).
Si el proceso de "deambular" está configurado correctamente, puede asegurarse de que se logre esta proporcionalidad (entre el tiempo empleado y la altura de distribución) .

Algoritmo


Ahora definamos y describamos la tarea en términos más formales. Sea s = (s1, s2, ..., sM) la distribución estacionaria deseada. Queremos crear una cadena de Markov con una distribución tan estacionaria. Comenzamos con una cadena de Markov arbitraria con estados M con la matriz de transición P, de modo que pij representa la probabilidad de transición del estado i a j.

Intuitivamente, sabemos cómo recorrer la cadena de Markov, pero la cadena de Markov no tiene la distribución estacionaria requerida. Esta cadena tiene alguna distribución estacionaria (que no necesitamos). Nuestro objetivo es cambiar la forma en que deambulamos por la cadena de Markov para que la cadena tenga la distribución estacionaria deseada.

Para hacer esto:

  1. Comience con un estado inicial aleatorio i.
  2. Seleccione aleatoriamente un nuevo estado asumido mirando las probabilidades de transición en la i-ésima fila de la matriz de transición P.
  3. Calcule una medida llamada probabilidad de decisión, que se define como: aij = min (sj.pji / si.pij, 1).
  4. Ahora lanza una moneda, que aterriza en la superficie del águila con probabilidad aij. Si un águila cae, acepte la oferta, es decir, vaya al siguiente estado, de lo contrario, rechace la oferta, es decir, permanezca en el estado actual.
  5. Repite muchas veces.

Después de una gran cantidad de pruebas, esta cadena convergerá y tendrá una distribución estacionaria. Entonces podemos usar estados de cadena como muestra de cualquier distribución.

Al hacer esto para muestrear la distribución beta, el único momento en que debe usar la densidad de probabilidad es buscar la probabilidad de tomar una decisión. Para hacer esto, divida sj por si (es decir, la constante de normalización C se cancela).

Selección Beta




Ahora pasamos al problema del muestreo de la distribución beta.

Una distribución beta es una distribución continua en [0,1] y puede tener valores infinitos en [0,1]. Suponga que una cadena de Markov arbitraria P con estados infinitos en [0,1] tiene una matriz de transición P tal que pij = pji = todos los elementos en la matriz.

No necesitamos la Matriz P, como veremos más adelante, pero quiero que la descripción del problema sea lo más cercana posible al algoritmo que propusimos.

  • Comience con un estado inicial aleatorio que obtuve de una distribución uniforme en (0,1).
  • Seleccione aleatoriamente un nuevo estado asumido observando las probabilidades de transición en la i-ésima fila de la matriz de transición P. Supongamos que elegimos otro estado Unif (0,1) como el estado asumido j.
  • Calcule la medida, que se llama la probabilidad de tomar una decisión:


Lo que se simplifica a:

Desde pji = pij, y donde

  • Ahora tira una moneda. Con probabilidad aij caerá un águila. Si cae un águila, entonces debe aceptar la oferta, es decir, pasar al siguiente estado. De lo contrario, vale la pena rechazar la oferta, es decir, permanecer en el mismo estado.
  • Repita la prueba muchas veces.

Código:


Es hora de pasar de la teoría a la práctica. Escribiremos nuestra muestra beta en Python.

 impo rt rand om # Lets define our Beta Function to generate s for any particular state. We don't care for the normalizing constant here. def beta_s(w,a,b): return w**(a-1)*(1-w)**(b-1) # This Function returns True if the coin with probability P of heads comes heads when flipped. def random_coin(p): unif = random.uniform(0,1) if unif>=p: return False else: return True # This Function runs the MCMC chain for Beta Distribution. def beta_mcmc(N_hops,a,b): states = [] cur = random.uniform(0,1) for i in range(0,N_hops): states.append(cur) next = random.uniform(0,1) ap = min(beta_s(next,a,b)/beta_s(cur,a,b),1) # Calculate the acceptance probability if random_coin(ap): cur = next return states[-1000:] # Returns the last 100 states of the chain 

Compare los resultados con la distribución beta real.

 impo rt num py as np import pylab as pl import scipy.special as ss %matplotlib inline pl.rcParams['figure.figsize'] = (17.0, 4.0) # Actual Beta PDF. def beta(a, b, i): e1 = ss.gamma(a + b) e2 = ss.gamma(a) e3 = ss.gamma(b) e4 = i ** (a - 1) e5 = (1 - i) ** (b - 1) return (e1/(e2*e3)) * e4 * e5 # Create a function to plot Actual Beta PDF with the Beta Sampled from MCMC Chain. def plot_beta(a, b): Ly = [] Lx = [] i_list = np.mgrid[0:1:100j] for i in i_list: Lx.append(i) Ly.append(beta(a, b, i)) pl.plot(Lx, Ly, label="Real Distribution: a="+str(a)+", b="+str(b)) pl.hist(beta_mcmc(100000,a,b),normed=True,bins =25, histtype='step',label="Simulated_MCMC: a="+str(a)+", b="+str(b)) pl.legend() pl.show() plot_beta(0.1, 0.1) plot_beta(1, 1) plot_beta(2, 3) 




Como puede ver, los valores son muy similares a la distribución beta. Por lo tanto, la red MCMC ha alcanzado un estado estacionario

En el código anterior, creamos una muestra beta, pero el mismo concepto se aplica a cualquier otra distribución de la que queramos hacer una selección.

Conclusiones




Fue una gran publicación. Felicitaciones si lo lees hasta el final.

En esencia, los métodos MCMC pueden ser complejos, pero nos brindan una gran flexibilidad. Puede seleccionar desde cualquier función de distribución utilizando la selección a través de MCMC. Típicamente, estos métodos se usan para tomar muestras de distribuciones posteriores.

También puede usar MCMC para resolver problemas con un gran espacio de estado. Por ejemplo, en un problema de mochila o para descifrar. Intentaré brindarte ejemplos más interesantes en la próxima publicación. Estén atentos.

De los editores


Source: https://habr.com/ru/post/460497/


All Articles