Метод Монте-Карло Метрополис: схема, баланс, сходимость

Метод Монте-Карло Метрополис - это способ оценивать средние по сложному вероятностному распределению, когда прямое интегрирование невозможно, а плотность известна только с точностью до нормировки. Его предложили Николас Метрополис, Арианна и Маршалл Розенблют и Эдвард и Августа Теллер в 1953 году для расчётов в статистической физике: нужно было считать термодинамические средние в системе многих частиц, где число конфигураций астрономически велико. Идея - не перебирать все состояния, а построить случайное блуждание (цепь Маркова), которое посещает конфигурации с частотой, пропорциональной их вероятности. Тогда обычное среднее вдоль траектории сходится к среднему по распределению.

Зачем нужен метод Монте-Карло Метрополис

Классический Монте-Карло оценивает интеграл $\langle A \rangle = \int A(x)\,\pi(x)\,dx$, генерируя точки $x$ из распределения $\pi$ и усредняя $A(x)$. Проблема в том, что в физике и байесовской статистике $\pi(x)$ обычно имеет вид $\pi(x) = \frac{1}{Z} e^{-\beta E(x)}$ - распределение Больцмана, где $Z = \sum_x e^{-\beta E(x)}$ - статистическая сумма (partition function). Считать $Z$ - значит просуммировать по всем состояниям, а их, например для модели Изинга на решётке $30\times 30$, уже $2^{900}$. Прямая выборка из такого $\pi$ невозможна.

Метод Метрополиса обходит это: вероятность принятия зависит только от отношения $\pi(x')/\pi(x)$, в котором нормировка $Z$ сокращается. Нам достаточно уметь считать энергию $E(x)$ - а не статистическую сумму.

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

Распределение Больцмана и принцип выборки по важности

Случайная выборка из равномерного распределения по конфигурациям почти бесполезна: подавляющее большинство состояний имеет ничтожный больцмановский вес $e^{-\beta E}$, и оценка $\langle A \rangle$ получается из единичных «удачных» точек с огромной дисперсией. Метод Монте-Карло Метрополис реализует выборку по важности (importance sampling): он генерирует конфигурации сразу с частотой, пропорциональной $\pi(x)$. Тогда среднее по выборке

$$\langle A \rangle \approx \frac{1}{N}\sum_{n=1}^{N} A(x_n)$$

не требует весов - каждый посещённый $x_n$ уже учтён с правильной частотой. Здесь $\beta = 1/(k_B T)$ - обратная температура: при больших $\beta$ (низкая $T$) распределение концентрируется около минимумов энергии, при малых $\beta$ становится почти равномерным.

Общая схема алгоритма Метрополиса

Пусть текущее состояние цепи - $x$. Один шаг (Monte Carlo step) состоит из трёх действий:

1. Предложить кандидата $x'$ симметричным пробным ходом - например, перевернуть случайный спин в модели Изинга или сдвинуть координату $x' = x + \varepsilon$, $\varepsilon \sim N(0, \sigma^2)$.
2. Вычислить изменение энергии $\Delta E = E(x') - E(x)$ и вероятность принятия

$$\alpha = \min\!\left(1,\; \frac{\pi(x')}{\pi(x)}\right) = \min\!\left(1,\; e^{-\beta \Delta E}\right).$$

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

Логика прозрачна: ход «вниз» по энергии ($\Delta E \le 0$) принимается всегда, ход «вверх» ($\Delta E > 0$) - с вероятностью $e^{-\beta\Delta E}$. Именно эти редкие подъёмы дают системе возможность выбираться из локальных ям и исследовать всё пространство. Это исходный метод Монте-Карло Метрополис; обобщение на несимметричные пробные ходы - алгоритм Метрополиса-Гастингса с дополнительным множителем $q(x\mid x')/q(x'\mid x)$.

Почему это работает: детальный баланс

Чтобы цепь Маркова имела стационарным распределением именно $\pi$, достаточно выполнения условия детального баланса (detailed balance):

$$\pi(x)\,P(x \to x') = \pi(x')\,P(x' \to x),$$

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

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

Равенство выполнено, значит $\pi$ - инвариантная мера. Если вдобавок цепь неприводима (из любого состояния достижимо любое) и апериодична, она эргодична: распределение $x_n$ при $n \to \infty$ стремится к $\pi$ независимо от старта. Детальный баланс - достаточное, но не необходимое условие; на нём держится корректность всего метода.

Модель Изинга: канонический пример

Историческое и до сих пор учебное применение - модель Изинга. Спины $s_i \in \{-1, +1\}$ на решётке, энергия

$$E(s) = -J\sum_{\langle i,j\rangle} s_i s_j - h\sum_i s_i,$$

где сумма по соседним парам, $J$ - обменное взаимодействие, $h$ - внешнее поле. Метод Монте-Карло Метрополис здесь работает так: выбираем случайный спин, считаем $\Delta E$ от его переворота (зависит только от ближайших соседей - это локально и дёшево), принимаем переворот с вероятностью $\min(1, e^{-\beta\Delta E})$. Усредняя намагниченность $M = \sum_i s_i$ и энергию вдоль траектории, получаем термодинамические средние. Вблизи критической температуры $\beta_c$ цепь резко замедляется (критическое замедление) - соседние конфигурации сильно скоррелированы, и тут переходят к кластерным алгоритмам (Вольфа, Свендсена-Ванга).

Burn-in, автокорреляция и оценка погрешности

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

• Автокорреляция. Соседние $x_n$ зависимы, поэтому оценивают интегральное время автокорреляции $\tau$. Эффективное число независимых выборок $\mathrm{ESS} = N/(2\tau + 1)$ много меньше $N$. Стандартная ошибка среднего - не $\sigma/\sqrt{N}$, а $\sigma/\sqrt{\mathrm{ESS}}$.
• Acceptance rate. Доля принятых ходов. Слишком высокая (близко к 1) - шаг $\sigma$ мал, цепь еле движется; слишком низкая - шаг велик, почти всё отвергается. Для гауссовых пробных ходов оптимум около $0.234$ в многомерном случае и $\approx 0.44$ в одномерном.

Чтобы поймать застревание в одной моде мультимодального распределения, запускают несколько независимых цепей и сравнивают их статистику ($\hat{R}$ Гельмана-Рубина $< 1.01$).

Связь с имитацией отжига

Если медленно увеличивать $\beta$ (понижать «температуру») по ходу выборки, метод Метрополиса превращается в имитацию отжига (simulated annealing) - алгоритм оптимизации. При высокой $T$ система свободно гуляет и не застревает в локальных минимумах энергии; при постепенном охлаждении она «оседает» в глобальном минимуме. Это прямое следствие того, что распределение Больцмана при $\beta \to \infty$ концентрируется на состояниях с минимальной $E$. Один и тот же шаг принятия $\min(1, e^{-\beta\Delta E})$ служит и для расчёта средних, и для поиска оптимума.

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

• Пытаются вычислить статистическую сумму $Z$. Весь смысл метода в том, что $Z$ сокращается в отношении $\pi(x')/\pi(x)$ - её считать не нужно и обычно невозможно.
• Не отбрасывают burn-in. Начальные нетермализованные конфигурации смещают оценку средних, особенно если старт далёк от равновесия.
• Считают шаги цепи независимыми. Используют $\sigma/\sqrt{N}$ для погрешности, игнорируя автокорреляцию, - реальная неопределённость занижается в $\sqrt{2\tau+1}$ раз.
• Гонятся за высоким acceptance rate. Приём $>0.9$ означает, что цепь почти стоит на месте; целиться надо в умеренный диапазон, а не в максимум.
• Запускают одну цепь у критической точки. Вблизи $\beta_c$ из-за критического замедления одиночная цепь не успевает разойтись - нужны кластерные методы или несколько стартов.

FAQ

Чем метод Монте-Карло Метрополис отличается от обычного Монте-Карло? Обычный Монте-Карло генерирует независимые точки из распределения напрямую и усредняет. Метод Метрополиса строит зависимую цепь Маркова - это нужно, когда напрямую сэмплировать $\pi$ нельзя, а известно лишь $\pi$ с точностью до нормировки.

Почему ход с ростом энергии вообще принимается? Иначе цепь скатилась бы в ближайший локальный минимум и осталась там. Подъёмы с вероятностью $e^{-\beta\Delta E}$ обеспечивают неприводимость и дают правильный больцмановский вес высокоэнергетическим состояниям при конечной температуре.

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

Коротко

Метод Монте-Карло Метрополис оценивает средние по распределению Больцмана $\pi(x)\propto e^{-\beta E(x)}$, строя цепь Маркова: предлагается симметричный пробный ход, кандидат принимается с вероятностью $\min(1, e^{-\beta\Delta E})$. Нормировка $Z$ при этом сокращается, поэтому достаточно считать только энергию. Корректность гарантирует детальный баланс, а $\pi$ оказывается стационарным распределением цепи. Канонический пример - расчёт термодинамических средних в модели Изинга; качество выборки контролируют burn-in, автокорреляция и acceptance rate, а понижение температуры превращает метод в имитацию отжига.