Алгоритм Метрополиса-Гастингса: как работает MCMC

Алгоритм Метрополиса-Гастингса - базовый и до сих пор самый часто упоминаемый метод MCMC (Markov chain Monte Carlo). Его придумали Николас Метрополис с соавторами в 1953 году для расчётов в статистической физике, а в 1970 году Уилфред Кит Гастингс обобщил схему до несимметричных пропозиций - отсюда двойное название. Идея простая: построить цепь Маркова, чьё стационарное распределение совпадает с целевым $\pi(x)$, и собирать из неё реализации. Это нужно, когда плотность $\pi(x) \propto f(x)$ известна только с точностью до нормировочной константы - то есть мы умеем считать $f(x)$, но не интеграл $\int f$.

Зачем сэмплировать из $\pi(x) \propto f(x)$

Типовая постановка приходит из байесовской статистики. Апостериорная плотность параметра $\theta$ при наблюдаемых данных $D$ записывается через формулу Байеса:

$$\pi(\theta \mid D) = \frac{L(D \mid \theta)\,\pi_0(\theta)}{\int L(D \mid \theta)\,\pi_0(\theta)\,d\theta}.$$

Числитель $f(\theta) = L(D\mid\theta)\,\pi_0(\theta)$ считается явно - это произведение правдоподобия и априорного распределения. А знаменатель - интеграл по всему пространству параметров - может быть многомерным и аналитически не браться. Прямое сэмплирование (inverse-CDF, rejection sampling) либо требует нормировки, либо плохо масштабируется по размерности. MCMC обходит проблему: строит цепь, чьё инвариантное распределение и есть $\pi$, не вычисляя нормировку.

Если ты уже определился с конкретной целью и хочешь подобрать шаг и пропозицию - заполни поля ниже, и мы соберём acceptance ratio, подскажем sigma и проверим autocorrelation.

Общая схема: пропозиция и шаг принятия

Пусть текущее состояние цепи - $x$. На каждом шаге:

1. Сгенерировать кандидата $x' \sim q(x' \mid x)$ из произвольного пропозиционного распределения $q$, из которого мы умеем сэмплировать.
2. Посчитать acceptance ratio

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

1. Бросить $u \sim \mathrm{Uniform}(0, 1)$ и принять кандидата ($x_{n+1} = x'$), если $u \le \alpha$. Иначе остаться: $x_{n+1} = x$.

Важный момент: $\pi$ входит в $\alpha$ только отношением $\pi(x')/\pi(x)$. Если $\pi \propto f$, то нормировочная константа сокращается - поэтому достаточно считать $f$. Это ключ ко всему MCMC: мы никогда не интегрируем плотность, мы только сравниваем её значения в двух точках.

Стационарность через detailed balance

Цепь имеет инвариантным распределение $\pi$, если выполнено условие детального баланса (detailed balance):

$$\pi(x)\, P(x \to x') = \pi(x')\, P(x' \to x)\quad \text{для всех } x, x',$$

где $P(x \to x') = q(x' \mid x)\,\alpha(x, x')$ - вероятность перехода. Проверим. Без ограничения общности пусть $\pi(x')\,q(x\mid x') \le \pi(x)\,q(x'\mid x)$ - тогда $\alpha(x, x') = \frac{\pi(x')\,q(x\mid x')}{\pi(x)\,q(x'\mid x)}$ и $\alpha(x', x) = 1$. Подставляем:

$\pi(x)\,q(x'\mid x)\,\alpha(x,x') = \pi(x')\,q(x\mid x') = \pi(x')\,q(x\mid x')\,\alpha(x',x).$

Условие detailed balance выполнено, значит $\pi$ - инвариантная мера цепи. При неприводимости и апериодичности (для $q > 0$ это почти всегда так) цепь эргодична и распределение $x_n$ при $n \to \infty$ стремится к $\pi$.

Частный случай: алгоритм Метрополиса

Если пропозиция симметрична - $q(x' \mid x) = q(x \mid x')$, как у нормального шага $x' = x + \varepsilon$, $\varepsilon \sim N(0, \sigma^2)$ - отношение пропозиций сокращается, и формула упрощается до

$$\alpha(x, x') = \min\!\left(1,\; \frac{\pi(x')}{\pi(x)}\right).$$

Это исходный алгоритм Метрополиса 1953 года. Логика интуитивная: если кандидат в более «плотной» области, принимаем всегда; если в менее плотной, принимаем с вероятностью отношения плотностей. Гастингс снял требование симметрии: добавил поправочный множитель $q(x\mid x')/q(x'\mid x)$, и стало можно использовать любые пропозиции, в том числе direction-biased (Langevin, independence sampler).

Типовые пропозиции

Random Walk Gaussian (нормальное блуждание). $q(x'\mid x) = N(x'; x, \sigma^2 I)$. Симметрично, простейший выбор. Подбор $\sigma$ критичен: слишком маленький - цепь почти всё принимает, но медленно гуляет (высокий autocorrelation); слишком большой - почти всё отвергает и стоит на месте. Оптимальный target acceptance rate для $d$-мерного гауссовского таргета - $\approx 0.234$ (Roberts, Gelman, Gilks 1997), на $d = 1$ ближе к $0.44$.

Langevin (MALA - Metropolis-adjusted Langevin algorithm). Пропозиция $x' = x + \tfrac{\sigma^2}{2}\nabla \log \pi(x) + \sigma \varepsilon$ использует градиент log-плотности и сдвигает кандидата в сторону мод. Поправка Метрополиса-Гастингса нужна, потому что дискретизация SDE Ланжевена смещена. Оптимальный target acceptance rate - $0.574$.

Hamiltonian Monte Carlo (HMC). Вводится вспомогательный импульс $p \sim N(0, M)$, и цепь движется по гамильтоновой траектории $H = -\log\pi(x) + \tfrac12 p^T M^{-1} p$ интегратором Верле. Шаг принятия компенсирует ошибку интегратора. На гладких многомерных плотностях HMC даёт ESS на одну итерацию на порядки выше, чем Random Walk, и стандартно используется в Stan, PyMC, NumPyro. Target acceptance - $0.65$–$0.80$.

Burn-in и проверки сходимости

Цепь стартует из произвольной точки, и первые сотни-тысячи итераций система ещё не «забыла» начальное состояние - это burn-in (или warm-up), и его выкидывают. Длина зависит от шага и от того, насколько $x_0$ далеко от типичных значений $\pi$. Эмпирически: 10–50% выборки.

Главные численные диагностики сходимости и качества:

• $\hat{R}$ Гельмана-Рубина. Запускают $m$ независимых цепей, сравнивают внутрицепную дисперсию $W$ с межцепной $B$ и считают $\hat{R} = \sqrt{(W + B/n)/W}$. Хорошее значение - $\hat{R} < 1.01$. Если $\hat{R}$ заметно больше единицы, цепи ещё не пришли к общему распределению.
• Autocorrelation $\rho_k$. Чем медленнее $\rho_k$ убывает с лагом $k$, тем сильнее зависимы соседние реализации. Интегральное время автокорреляции $\tau = 1 + 2\sum_{k\ge 1}\rho_k$ показывает, сколько шагов между «эффективно независимыми» сэмплами.
• Effective sample size. $\mathrm{ESS} = N/\tau$ - число «эффективных» независимых наблюдений из $N$ собранных. Если $\mathrm{ESS}$ мал по сравнению с $N$, цепь либо плохо настроена, либо смешивается медленно из-за самой геометрии $\pi$ (мультимодальность, узкие хребты).

Где применяется

• Байесовская статистика. Любая нетривиальная апостериорная плотность $\pi(\theta\mid D)$, где интеграл нормировки не берётся. Иерархические модели, GLM, спецификации с непрерывно-дискретными смесями.
• Статистическая физика. Исходная задача 1953 года - расчёт термодинамических средних в модели Изинга по распределению Больцмана $\pi(s) \propto e^{-\beta H(s)}$, отсюда же выросло и название.
• Машинное обучение. Выборка из энергетических моделей (Restricted Boltzmann Machines), регуляризация генеративных моделей, посэмплинг весов в Bayesian Neural Networks.
• Метрическая теория, теория чисел, биоинформатика. Везде, где нужно усреднить функцию по сложному распределению с явным $f(x)$.

Типовые задачи

1. Сгенерировать $10^4$ реализаций из $\pi(x) \propto e^{-x^4/4}$, оценить $E[X^2]$ - целевая плотность из физики, нормировка не выражается в элементарных функциях. Подходит Random Walk с $\sigma \approx 1.5$.
2. Сэмплировать апостериор $\pi(\theta\mid D)$ для логистической регрессии с нормальным prior. MALA с $\sigma = 0.05$ при $d = 20$ обычно даёт ESS $\ge 0.3 N$.
3. Симулировать конфигурации модели Изинга при $\beta = 0.4$ на решётке $32\times 32$ - Метрополис с одиночными flip-ами, acceptance $\approx 0.3$.

Частые ошибки

• Забывают, что $\pi$ нужна только с точностью до константы. Пытаются нормировать $\pi$ численно - ломают саму причину, по которой MCMC нужен.
• Не выбрасывают burn-in. Включают начальные итерации в оценки - смещение математического ожидания и недооценка дисперсии.
• Считают независимыми соседние реализации цепи. $\rho_k$ не нулевой - стандартная ошибка $\sigma/\sqrt{N}$ занижает реальную неопределённость; вместо $N$ нужно $\mathrm{ESS}$.
• Подгоняют $\sigma$ под максимум acceptance rate. Acceptance $> 0.9$ - это плохо: цепь почти не двигается. Целиться надо в $0.234$ для Random Walk, $0.574$ для MALA, $0.65$–$0.8$ для HMC.
• Запускают одну цепь и доверяют ей. Без $\hat{R}$ по нескольким независимым стартам нельзя обнаружить, что цепь застряла в одной моде мультимодального распределения.

FAQ

Чем отличается Метрополис от Метрополиса-Гастингса? Метрополис требует симметричной пропозиции $q(x'\mid x) = q(x\mid x')$, и тогда $\alpha = \min(1, \pi(x')/\pi(x))$. Гастингс добавил поправочный множитель $q(x\mid x')/q(x'\mid x)$ и снял ограничение симметрии. На практике под общим именем «Метрополиса-Гастингса» понимают полную версию.

Зачем acceptance ratio именно $\min(1, \cdot)$, а не $\cdot$ напрямую? Чтобы вероятность принятия не превышала единицы. Условие detailed balance выполняется и для не-усечённого отношения, но как вероятность принятия отношение $> 1$ интерпретировать нельзя - берут отсечение.

Можно ли использовать MCMC, если плотность $\pi$ имеет мультимодальную форму? Можно, но Random Walk будет надолго застревать в одной моде. Помогают tempering-схемы (parallel tempering, simulated annealing), HMC с длинными траекториями или специальные jumps - иначе ESS из моды в моду оказывается катастрофически мал.

Коротко

Алгоритм Метрополиса-Гастингса - это рецепт построения цепи Маркова, чьё стационарное распределение совпадает с заданным $\pi$. На каждом шаге генерируется кандидат из пропозиции $q$, считается acceptance ratio $\alpha = \min(1, \pi(x')q(x\mid x')/\pi(x)q(x'\mid x))$ и принимается решение принять или отвергнуть. Detailed balance обеспечивает инвариантность $\pi$. Метрополис - частный случай с симметричной $q$. Качество выборки контролируют autocorrelation, ESS и $\hat{R}$ Гельмана-Рубина. Алгоритм стоит в основе байесовской статистики и MCMC-инструментов вроде Stan и PyMC.