13.1: Formulación básica

La idea básica detrás del método MCMC es simple. Supongamos que tenemos un conjunto de estados etiquetados por un índice entero\(n \in \{0,1,2,\dots\}\), donde cada estado está asociado con una probabilidad\(\pi_{n}\). Por ejemplo, en mecánica estadística, para un sistema mantenido a temperatura constante\(T\), cada estado ocurre con probabilidad

\[\pi_n = \frac{\exp\left(-\frac{E_n}{k T}\right)}{Z},\]

donde\(E_{n}\) está la energía del estado,\(k\) es la constante de Boltzmann, y\(Z\) es la función de partición

\[Z = \sum_n \exp\left(-\frac{E_n}{k T}\right).\]

A partir de\(\{\pi_n\}\), nos gustaría calcular diversos valores de expectativa, que describen las propiedades termodinámicas del sistema. Por ejemplo, podríamos estar interesados en la energía promedio, que se define como

\[\langle E\rangle = \sum_n E_n \pi_n.\]

La forma más sencilla de encontrar\(\langle E\rangle\) es calcular explícitamente la suma anterior. Pero si el número de estados es muy grande, esto consume mucho tiempo prohibitivamente (a menos que haya una solución analítica tratable, y con frecuencia no la haya). Por ejemplo, si nos interesa describir\(1000\) distintos átomos cada uno teniendo\(2\) posibles niveles de energía, el número total de estados es\(2^{1000} \approx 10^{301}\). Tratar de calcular una suma sobre estos términos alucinantes tardaría más que la edad del universo.

El método MCMC evita este problema muestreando selectivamente los estados. Para lograrlo, diseñamos un proceso de Markov cuya distribución estacionaria es idéntica a las probabilidades dadas\(\pi_{n}\). Discutiremos cómo diseñar el proceso de Markov en la siguiente sección. Una vez que tenemos un proceso apropiado de Markov, podemos usarlo para generar una cadena larga de Markov, y usar esa cadena para calcular promedios móviles de nuestras cantidades deseadas, como\(\langle E\rangle\). Si la cadena de Markov es suficientemente larga, el promedio calculado de esta manera convergerá al verdadero valor de expectativa\(\langle E\rangle\).

El dato clave que hace que todo este trabajo sea que la longitud requerida de la cadena de Markov suele ser mucho menor que el número total de estados posibles. Para el problema antes mencionado de\(1000\) distintos átomos de dos niveles, hay\(10^{301}\) estados, pero una cadena de Markov de tan pocos como\(10^{6}\) pasos puede llegar dentro de varios por ciento del verdadero valor de\(\langle E \rangle\). (La precisión real variará de un sistema a otro). La razón de esto es que la gran mayoría de los estados son extremadamente improbables, y sus contribuciones a la suma que conduce\(\langle E \rangle\) son muy pequeñas. Una cadena de Markov puede obtener una buena estimación para\(\langle E \rangle\) muestreando los estados que tienen las probabilidades más altas, sin gastar mucho tiempo en estados de baja probabilidad.

13.1.1 El algoritmo de Metrópolis

El método MCMC requiere que diseñemos un proceso de Markov para que coincida con una distribución estacionaria dada\(\{\pi_n\}\). Este es un problema abierto, y generalmente hay muchas buenas maneras de lograr este objetivo. El método más común, llamado algoritmo Metrópolis, se basa en el principio de equilibrio detallado, que discutimos en el artículo sobre las cadenas de Markov. Para recapitular, el principio de equilibrio detallado establece que bajo circunstancias genéricas (que frecuentemente se cumplen en la física), las probabilidades de transición de un proceso de Markov están relacionadas con la distribución estacionaria por

\[P(n|m)\, \pi_m = P(m|n)\, \pi_n \quad \mathrm{for}\;\mathrm{all}\;m,n.\]

El algoritmo Metrópolis especifica el siguiente proceso de Markov:

1. Supongamos que en el paso\(k\), el sistema está en estado\(n\). Elegir aleatoriamente un estado candidato\(m\),, haciendo un paso aleatorio imparcial a través del espacio de los estados posibles. (Justo cómo se hace esta elección depende del sistema, y lo discutiremos a continuación).
2. Compara las probabilidades\(\pi_{n}\) y\(\pi_{m}\):
   • Si\(\pi_m \ge \pi_n\), acepta al candidato.
   • Si\(\pi_m < \pi_n\), acepta al candidato con probabilidad\(\pi_m/\pi_n\). De lo contrario, rechazar al candidato.
3. Si el candidato es aceptado, el estado en escalón\(k+1\) lo es\(m\). De lo contrario, el estado en escalón\(k+1\) permanece\(n\).
4. Repetir.

Con base en la descripción anterior, verifiquemos que la distribución estacionaria del proceso de Markov satisface la condición deseada de balance detallado. Considerar dos estados cualesquiera\(a,\: b\), y asumir sin pérdida de generalidad que\(\pi_a \le \pi_b\). Si partimos de\(a\), supongamos que elegimos un paso candidato\(a\rightarrow b\) con cierta probabilidad\(q\). Entonces, según las reglas de Metrópolis, la probabilidad de hacer realmente esta transición,\(a\rightarrow b\), es\(q\) veces la probabilidad de aceptación\(1\). Por otro lado, supongamos que partimos de\(b\) su lugar. Debido a que la elección del candidato es imparcial, elegiremos un paso candidato\(b\rightarrow a\) con la\(q\) misma probabilidad que en el caso anterior. De ahí que la probabilidad de transición para\(b\rightarrow a\) es\(q\) veces la probabilidad de aceptación de\(\pi_a/\pi_b\). Como resultado,

\[\begin{align}\left\{\begin{aligned}P(b|a) &= q \times 1 \\ P(a|b) &= q \times \frac{\pi_a}{\pi_b}.\end{aligned}\right. \quad \Rightarrow \quad P(a|b) \,\pi_b = P(b|a)\, \pi_a\end{align}\]

Dado que este razonamiento es arbitrario\(a,\: b\), el principio de equilibrio detallado implica que la distribución estacionaria de nuestro proceso de Markov sigue la distribución deseada\(\{\pi_n\}\).

Expresión en términos de energías

En física, el método MCMC se aplica comúnmente a los estados termodinámicos, para lo cual

\[\pi_n = \frac{\exp\left(-\frac{E_n}{kT}\right)}{Z}.\]

En tales casos, el algoritmo Metrópolis puede expresarse de manera equivalente en términos de las energías estatales:

1. Supongamos que en el paso\(k\), el sistema está en estado\(n\). Elegir aleatoriamente un estado candidato\(m\),, haciendo un paso aleatorio imparcial a través del espacio de los estados posibles.
2. Compara las energías\(E_{n}\) y\(E_{m}\):
   • Si\(E_m \le E_n\), acepta al candidato.
   • Si\(E_m > E_n\), acepta al candidato con probabilidad\(\exp\left[(E_n - E_m)/kT\right]\). De lo contrario, rechazar al candidato.
3. Si el candidato es aceptado, el estado en escalón\(k+1\) lo es\(m\). De lo contrario, el estado en escalón\(k+1\) permanece\(n\).
4. Repetir.

13.1.2 Paso a través del espacio estatal

Una forma de pensar sobre el algoritmo de Metrópolis es que toma un esquema para realizar un recorrido aleatorio imparcial por el espacio de los posibles estados (representados por nuestras elecciones de candidatos), y lo convierte en un esquema para realizar una caminata aleatoria sesgada. La caminata aleatoria sesgada corresponde a un proceso de Markov con la distribución estacionaria que nos interesa.

La forma en que se eligen los candidatos de Metrópolis (es decir, la parte de “caminata aleatoria imparcial”) varía de un sistema a otro, y una vez más hay múltiples esquemas válidos que podríamos emplear. Por ejemplo, supongamos que tenemos una colección de 6 átomos, donde cada átomo puede estar en el nivel etiquetado\(0\) o el nivel etiquetado\(1\). Cada estado del sistema general se describe mediante una lista de\(6\) símbolos,\(011001\) p. Entonces podemos hacer una caminata imparcial por el “espacio de estado” eligiendo aleatoriamente uno de los\(6\) átomos (con igual probabilidad), y volteándolo. Por ejemplo, podríamos optar por voltear el segundo átomo:

\[011001 \quad \stackrel{\mathrm{flip}\; \mathrm{second}\; \mathrm{atom}}{\longrightarrow}\quad 001001 \;\;\;\;\; (q = 1/6).\]

Si partimos del otro estado, el proceso inverso tiene la misma probabilidad:

\[001001 \quad \stackrel{\mathrm{flip}\; \mathrm{second}\; \mathrm{atom}}{\longrightarrow}\quad 011001 \;\;\;\;\; (q = 1/6).\]

De ahí que se diga que este esquema para caminar por el “espacio estatal” es imparcial. Tenga en cuenta que, para un esquema de caminata dado, no siempre es posible conectar cada dos estados por un solo paso; por ejemplo, en este caso no podemos pasar de\(000000\) a\(111111\) en un solo paso.

Hay más de un esquema de caminata posible; por ejemplo, un esquema diferente podría implicar elegir aleatoriamente dos átomos y voltearlos. Sea cual sea el esquema que escojamos, sin embargo, lo más importante es que la caminata debe ser imparcial: cada paso posible debe ocurrir con la misma probabilidad que el paso inverso. De lo contrario, la prueba anterior de que el algoritmo de Metrópolis satisface el balance detallado no funcionaría.