Parte II · Monte Carlo e MCMC — computar com aleatoriedade

draft

Se uma cadeia de Markov é uma máquina que caminha, Monte Carlo é a descoberta de que uma caminhada aleatória bem escolhida pode computar. Este capítulo segue a ideia do jogo de cartas de um matemático em convalescença até o motor que amostra quase todo modelo probabilístico moderno. A teoria da própria caminhada está no capítulo 01, Cadeias de Markov; aqui a caminhada é posta a trabalhar.

1. A ideia: amostrar em vez de resolver

Algumas quantidades são definidas por uma soma ou uma integral grande demais para avaliar diretamente — uma integral de alta dimensão, a probabilidade de uma configuração rara, o valor esperado de uma função sobre um espaço de estados enorme. A intuição de Monte Carlo é parar de tentar resolver a integral e, em vez disso, estimá-la por amostragem aleatória: sorteie muitas amostras, faça a média da coisa com que você se importa, e deixe a lei dos grandes números (ver cap. 01, §1) levar a média ao valor verdadeiro.

Para estimar uma integral \(I = \int f(x)\,dx\) sobre um domínio, sorteie \(N\) amostras independentes e faça a média de \(f\) (ponderada pela medida do domínio ou pela densidade de amostragem). A média amostral é um estimador não viesado de \(I\). A versão de brinquedo canônica: para estimar \(\pi\), jogue pontos uniformemente num quadrado unitário; a fração que cai dentro do quarto de círculo inscrito (\(x^2+y^2\le 1\)) se aproxima de \(\pi/4\).

Por que vence: o erro é cego à dimensão. O erro padrão de uma estimativa de Monte Carlo escala como \(\sigma/\sqrt{N}\) — erro \(\propto 1/\sqrt{N}\), então quadruplicar as amostras reduz o erro pela metade. Crucialmente, essa taxa não depende da dimensão do espaço sendo integrado. A quadratura baseada em grade sofre a maldição da dimensionalidade (erro \(\sim N^{-k/d}\) em \(d\) dimensões, sem esperança para \(d\) grande); Monte Carlo paga o mesmo \(1/\sqrt{N}\) quer o espaço tenha três dimensões ou três milhões. Esse único fato é por que ele domina a física estatística, as finanças e a inferência bayesiana — toda integração de alta dimensão disfarçada.

2. Origem: paciência, nêutrons e a bomba

O método nasceu dentro do Projeto Manhattan. Em 1946, se recuperando de uma doença (encefalite) e jogando paciência Canfield, Stanisław Ulam fez uma pergunta simples: qual é a probabilidade de um jogo montado dar certo? A combinatória era sem esperança — mas ele percebeu que podia simplesmente dar as cartas muitas vezes e contar as vitórias. Estimar a probabilidade por ensaios aleatórios repetidos.

Ulam descreveu a ideia a John von Neumann, e eles a apontaram para o problema real à frente deles: difusão e transporte de nêutrons em material físsil — um problema integro-diferencial de alta dimensão que resistia a todo ataque determinístico. Em março de 1947 von Neumann escreveu a Robert Richtmyer em Los Alamos com um esquema detalhado para rastrear histórias simuladas de nêutrons (posição, velocidade, colisão, absorção, fissão) por amostragem aleatória, explicitamente projetado para o ENIAC. As primeiras rodadas automatizadas de núcleo de fissão vieram em 1947–1948.

O nome é um codinome: Nicholas Metropolis propôs "Monte Carlo" em referência ao cassino de Mônaco — uma alusão a um tio de Ulam que "simplesmente tinha que ir a Monte Carlo" para jogar. A primeira descrição pública é:

Metropolis, N. & Ulam, S. (1949). The Monte Carlo Method. J. Amer. Statist. Assoc. 44(247):335–341.

Trabalhos posteriores em Los Alamos rodaram no MANIAC I (operacional em março de 1952), a máquina por trás do algoritmo de 1953 abaixo.

3. O problema que Monte Carlo puro não resolve — e a resposta do MCMC

Monte Carlo puro precisa sortear amostras independentes da distribuição-alvo. Mas as distribuições interessantes em física e estatística são justamente aquelas das quais você não pode amostrar diretamente: são de alta dimensão, e conhecidas apenas a menos de uma constante de normalização (você pode computar \(\propto \pi(x)\) mas não \(\pi(x)\) em si, porque o normalizador é sua própria integral intratável).

Markov Chain Monte Carlo (MCMC) vira o problema do avesso. Em vez de amostrar \(\pi\) diretamente, construa uma cadeia de Markov cuja distribuição estacionária seja \(\pi\) (distribuições estacionárias são definidas no cap. 01, §4). Rode a cadeia; uma vez que ela tenha misturado, os estados que visita são amostras de \(\pi\), e esperanças sob \(\pi\) se tornam médias temporais ao longo da trajetória (o teorema ergódico para cadeias de Markov). As amostras são correlacionadas em vez de independentes, mas vêm da distribuição certa — e você nunca precisou do normalizador.

Balanço detalhado é a maneira padrão de garantir que \(\pi\) é estacionária. Uma cadeia é reversível com respeito a \(\pi\) se, para todos os estados \(x, y\),

\[\pi(x)\,P(x\to y) = \pi(y)\,P(y\to x).\]

O balanço detalhado implica que \(\pi\) é invariante, então é o alvo de projeto de quase todo método clássico de MCMC.

4. Os algoritmos

Metropolis (1953). Nascido, de novo, na física estatística — uma simulação de 224 partículas de disco rígido para computar uma equação de estado, rodada no MANIAC I.

Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E. (1953). Equation of State Calculations by Fast Computing Machines. J. Chem. Phys. 21(6):1087–1092.

A regra: a partir do estado \(x\), proponha um movimento próximo para \(x'\). Com \(\pi \propto e^{-E/kT}\) (a distribuição de Boltzmann), seja \(\Delta E = E(x') - E(x)\) e aceite \(x'\) com probabilidade

\[\min\!\left(1,\; e^{-\Delta E/kT}\right) \;=\;\min\!\left(1,\; \frac{\pi(x')}{\pi(x)}\right).\]

Movimentos ladeira abaixo são sempre tomados; movimentos ladeira acima são tomados às vezes — apenas com frequência suficiente para satisfazer o balanço detalhado com respeito a \(\pi\). (Nota histórica: Arianna W. Rosenbluth escreveu a implementação no MANIAC e Marshall Rosenbluth desenvolveu boa parte do método; a atribuição a Metropolis é uma simplificação bem documentada.)

Metropolis–Hastings (1970). W. K. Hastings generalizou a regra para distribuições de proposta arbitrárias, possivelmente assimétricas \(q(x\to x')\):

Hastings, W. K. (1970). *Monte Carlo sampling methods using Markov chains and their applications.* Biometrika 57(1):97–109.

\[\alpha(x, x') = \min\!\left(1,\; \frac{\pi(x')\,q(x'\to x)}{\pi(x)\,q(x\to x')}\right).\]

A razão extra (a correção de Hastings) restaura o balanço detalhado quando a proposta não é simétrica; Metropolis é o caso especial de \(q\) simétrico. Este é o passo que levou o método para fora da física e para dentro da estatística de corrente principal.

Amostragem de Gibbs (1984). Geman & Geman, trabalhando em restauração bayesiana de imagens sobre campos aleatórios de Markov, introduziram a atualização de uma variável por vez a partir de sua condicional completa \(\pi(x_i \mid x_{-i})\) — um caso especial de Metropolis–Hastings cuja probabilidade de aceitação é sempre 1.

Geman, S. & Geman, D. (1984). *Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images.* IEEE Trans. PAMI 6(6):721–741.

Popularizado para trabalho bayesiano geral por Gelfand & Smith (1990), disparou a explosão do MCMC pela estatística nos anos 1990.

5. Onde o MCMC aparece hoje

  • Física estatística — o domínio fundador: modelos de rede, mais famosamente o

    modelo de Ising (spins numa rede, amostrados da distribuição de Boltzmann para achar magnetização e transições de fase), QCD de rede, equações de estado de fluidos.

  • Inferência bayesiana — o cavalo de batalha para sortear de uma posteriori

    \(p(\theta\mid \text{dados}) \propto p(\text{dados}\mid\theta)\,p(\theta)\) quando a evidência é intratável. Este é exatamente o caso "conhecido a menos de uma constante" para o qual o MCMC foi construído, e a razão pela qual a estatística bayesiana aplicada foi transformada depois de 1990.

  • Biologia computacional — filogenética bayesiana (MrBayes, BEAST), genética de populações,

    alinhamento de sequências.

  • Finanças — precificação Monte Carlo de derivativos dependentes de trajetória, value-at-risk,

    estimação de volatilidade estocástica.

  • ML probabilístico moderno — o motor geral para amostrar distribuições de alta dimensão não

    normalizadas: modelos de tópicos (LDA), modelos baseados em energia e generativos, e os descendentes cientes de gradiente que o escalam — Hamiltonian Monte Carlo, o No-U-Turn Sampler (NUTS, por trás do Stan/PyMC), e dinâmica de Langevin / SGLD.

Por que importa pra Stack

As camadas de dados e de IA da Stack se apoiam exatamente no objeto que o MCMC amostra: uma distribuição conhecida apenas a menos de uma constante, sobre um espaço grande demais para enumerar. Qualquer estimativa bayesiana, qualquer ranqueamento probabilístico, qualquer amostrador generativo em um componente Koder é, por baixo, uma cadeia de Markov engenheirada para ter a resposta como sua distribuição estacionária — a mesma construção que o PageRank usa para ranquear a web (capítulo 03) e a mesma que os modelos de linguagem especializam para amostrar texto (Compêndio de IA, cap. 35).

Fontes

  • Metropolis, N. & Ulam, S. (1949). The Monte Carlo Method. J. Amer. Statist. Assoc. 44(247):335–341.
  • Metropolis, Rosenbluth, Rosenbluth, Teller & Teller (1953). Equation of State Calculations by Fast Computing Machines. J. Chem. Phys. 21(6):1087–1092.
  • Hastings, W. K. (1970). Monte Carlo sampling methods using Markov chains and their applications. Biometrika 57(1):97–109.
  • Geman, S. & Geman, D. (1984). Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images. IEEE Trans. PAMI 6(6):721–741.
  • Eckhardt, R. (1987). Stan Ulam, John von Neumann, and the Monte Carlo Method. Los Alamos Science, Special Issue 15:131–141.
  • Gelfand, A. E. & Smith, A. F. M. (1990). Sampling-Based Approaches to Calculating Marginal Densities. J. Amer. Statist. Assoc. 85(410):398–409.