MCMC и сэмплирование

1953 год: Метрополис и коллеги запустили первый MCMC-алгоритм на компьютере MANIAC, решая задачу из ядерной физики с 56 частицами. Сегодня тот же принцип используется в Stan для моделей с миллионом параметров, в LIGO для оценки параметров гравитационных волн, и в каждом диффузионном генеративном AI. MCMC - рабочая лошадка современной байесовской статистики.

• **Байесовские модели (Stan, PyMC)** - иерархические модели, смешанные эффекты - сэмплируются через NUTS (HMC)
• **Физика частиц и астрофизика** - оценка параметров из огромных датасетов через MCMC; обработка неопределённостей
• **Генеративные модели (Langevin Diffusion)** - сэмплирование из апостериорного распределения нейронной сети = уравнение Ланжевена

Предварительные знания

• Discrete-Time Markov Chains
• Martingale Inequalities and Convergence

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

Как сэмплировать из сложного распределения π(x), которое мы знаем с точностью до константы нормировки? Алгоритм Метрополиса-Гастингса строит цепь Маркова, стационарным распределением которой является π. Идея проста: предлагать случайные шаги и принимать или отклонять их с вероятностью, зависящей от отношения плотностей.

**Алгоритм Метрополиса-Гастингса:** 1. Из текущего x предложить y ~ q(y|x) (предложенный шаг) 2. Принять y с вероятностью α = min(1, π(y)q(x|y) / (π(x)q(y|x))) 3. Иначе остаться в x **Гарантия:** Цепь имеет стационарное распределение π (детальный баланс). Для симметричного q(y|x) = q(x|y): α = min(1, π(y)/π(x)) - алгоритм Метрополиса.

**Детальный баланс** обеспечивает правильность: π(x)·q(y|x)·α(x,y) = π(y)·q(x|y)·α(y,x). Это достаточное (но не необходимое!) условие для π быть стационарным распределением. **Размер шага** критичен: слишком маленький - медленное перемешивание, слишком большой - низкая acceptance rate. Оптимальная acceptance rate для случайного walk в d измерениях ≈ 23%.

Сэмплер Гиббса

При сэмплировании из многомерного распределения π(x₁,...,x_d) часто условные распределения π(x_i | x_{-i}) известны явно, а совместное - нет. **Сэмплер Гиббса** использует это: поочерёдно обновляет каждую координату из её условного распределения. Это специальный случай Метрополиса с acceptance rate = 1 для каждого шага.

**Сэмплер Гиббса:** Для распределения π(x₁,...,x_d): 1. Для k = 1,...,d: Сэмплировать x_k^{new} ~ π(x_k | x_1^{new},...,x_{k-1}^{new}, x_{k+1},...,x_d) 2. Повторить **Гарантия:** Детальный баланс выполнен для каждого условного обновления - π является стационарным распределением. **Условие применимости:** Нужно уметь эффективно сэмплировать из всех условных распределений.

**Блокированный Гиббс:** вместо поочерёдного обновления каждой координаты - обновление блоков переменных. Это снижает автокорреляцию при сильных зависимостях внутри блока. В байесовских иерархических моделях (Stan, PyMC) сэмплер Гиббса используется для конъюгатных компонент.

Hamiltonian Monte Carlo

Стандартный Метрополис «блуждает» случайно, не используя структуру целевого распределения. **Hamiltonian Monte Carlo (HMC)** использует градиент log π для «умного» движения: вместо случайного walk - детерминированная гамильтонова динамика, которая исследует пространство вдоль изолиний распределения. Это основа современных байесовских фреймворков (Stan, PyMC, NumPyro).

**Hamiltonian Monte Carlo:** 1. Добавить импульс p ~ N(0, M) (вспомогательная переменная) 2. Симулировать гамильтонову динамику (leapfrog): H(x,p) = -log π(x) + p²/2M x_{t+ε} = x_t + ε·p_t/M p_{t+ε} = p_t + ε·∇log π(x_{t+ε}) 3. Принять/отклонить по критерию Метрополиса (исправляет ошибку дискретизации) **Преимущество:** Предложения далеко от текущей точки, но всё ещё принимаются с высокой вероятностью.

**No-U-Turn Sampler (NUTS)** - адаптивная версия HMC, используемая в Stan. NUTS автоматически выбирает L (число шагов leapfrog), останавливаясь когда путь «разворачивается» (U-turn). Это устраняет необходимость ручной настройки L - основного гиперпараметра HMC.

Диагностика сходимости

Как убедиться, что MCMC цепь сошлась к π? Мы никогда не знаем точно - только накапливаем свидетельства. Есть несколько диагностик: визуальные (trace plots, autocorrelation plots) и количественные (R-hat Гелмана-Рубина, ESS, BFMI для HMC).

**Ключевые диагностики MCMC:** **R-hat (Gelman-Rubin):** Запустить несколько цепей из разных начальных точек. R̂ = sqrt(Var_between / Var_within + (N-1)/N) R̂ < 1.01 - хорошо, R̂ > 1.1 - не сошлась **ESS (Effective Sample Size):** ESS = N / (1 + 2·∑_{k=1}^∞ ρ_k) где ρ_k - автокорреляция лага k. **BFMI (Bayesian Fraction Missing Information):** специфично для HMC - проверяет «энергию» цепи.

**Практические советы:** - Всегда запускать несколько цепей (минимум 4) из разных начальных точек - Удалять burnin (обычно 50% итераций или до визуальной стабилизации) - Проверять R̂ < 1.01 и ESS > 400 для надёжных оценок - Для HMC: проверять acceptance rate 60-90% и BFMI > 0.2

Ключевые идеи

• **Метрополис-Гастингс** - предложить шаг, принять с prob = min(1, π(y)/π(x)·q(x|y)/q(y|x)); детальный баланс
• **Гиббс** - поочерёдно сэмплировать из условных распределений; acceptance rate = 1; медленен при высоких корреляциях
• **HMC** - использует градиент log π; leapfrog интегрирование; эффективен в высоких размерностях; NUTS - адаптивная версия
• **Диагностика** - R̂ < 1.01, ESS > 400, trace plots; всегда запускать несколько цепей из разных точек

Связанные темы

MCMC объединяет теорию цепей Маркова с вычислительной байесовской статистикой:

• Цепи Маркова — MCMC строит цепь Маркова со стационарным распределением π; теория сходимости CTMC применима
• Мартингальные неравенства — ESS и скорость сходимости MCMC выводятся через мартингальные аргументы
• Stochastic в ML — SGD и Langevin dynamics - MCMC для оптимизации; Score matching - родственно HMC

Вопросы для размышления

• MCMC и SGD оба используют случайность для навигации в пространстве параметров. В чём фундаментальная разница их целей?
• Диффузионные генеративные модели используют Langevin MCMC для сэмплирования. Почему они работают быстрее обычного MCMC?
• Если R-hat = 2.0, что это говорит о топологии целевого распределения? Как переформулировать задачу?

Связанные уроки

• prob-04-bayes