Метод Монте-Карло по схеме марковских цепей

Материал из MachineLearning.

(Различия между версиями)

Текущая версия

Статья подготовлена с использованием модели OpenAI GPT‑5.6 Sol с уровнем рассуждений High и проверена участником Д.О. Кистанов 16:57, 18 июля 2026 (MSK)

Промпт приводится полностью в Обсуждение:Метод Монте-Карло по схеме марковских цепей

Содержание

1 История
2 Основная идея
- 2.1 Популярная интерпретация
3 Математические основы
- 3.1 Инвариантность и детальный баланс
- 3.2 Автокорреляция и центральная предельная теорема
4 Основные алгоритмы
5 Диагностика и оценка точности
6 Практическое построение модели и алгоритма
7 Трудности и ограничения
8 Современные направления исследований
9 Применения
10 См. также
11 Примечания
12 Литература
13 Ссылки

Методы Монте-Карло на основе марковских цепей (MCMC, от англ. Markov chain Monte Carlo) — класс вычислительных методов для получения выборок из сложных распределений вероятностей и приближённого вычисления математических ожиданий. В MCMC строят марковскую цепь, стационарным распределением которой служит заданное целевое распределение. После достаточно длительного моделирования состояния цепи используют как зависимую выборку из этого распределения.

MCMC применяется, когда непосредственное генерирование независимых наблюдений затруднено, но плотность целевого распределения можно вычислять хотя бы с точностью до неизвестного нормирующего множителя. Это особенно важно в байесовском выводе, где апостериорная плотность часто известна в виде произведения правдоподобия и априорной плотности, тогда как нормирующая константа представляет собой многомерный интеграл. Другие крупные области применения — статистическая физика, вычислительная химия, обработка изображений, вероятностные графические модели, биоинформатика и машинное обучение.

В отличие от обычного метода Монте-Карло, последовательные значения MCMC, как правило, коррелированы. Поэтому число итераций само по себе не характеризует точность: необходимо оценивать автокорреляцию, эффективный размер выборки и ошибку Монте-Карло. Кроме того, инвариантность целевого распределения ещё не гарантирует быстрого достижения стационарного режима; медленное перемешивание и пропуск мод целевого распределения остаются основными практическими трудностями.

История

Первый широко известный алгоритм этого класса был опубликован в 1953 году группой Н. Метрополиса для моделирования двумерной системы жёстких дисков на компьютере MANIAC.^[1] В 1970 году У. Хастингс сформулировал общий алгоритм с несимметричными предложениями и обсудил оценивание ошибки Монте-Карло.^[1]

В 1984 году С. и Д. Геманы применили гиббсовские распределения и стохастическую релаксацию к байесовскому восстановлению изображений,^[1] а в 1987 году С. Дюэйн с соавторами предложил гибридный, позднее названный гамильтоновским, метод Монте-Карло для решёточной теории поля.^[1] В конце 1980-х и начале 1990-х годов гиббсовское сэмплирование и другие марковские методы вошли в широкую практику байесовской статистики.^[1] Обзор Л. Тирни 1994 года связал алгоритмические конструкции с теорией марковских цепей на общих пространствах состояний.^[1]

С 1990-х годов появились трансразмерные, адаптивные и псевдомаргинальные методы; в 2010-х широкое распространение получили автоматизированные HMC/NUTS и стохастические градиентные алгоритмы. Современная работа сосредоточена не только на создании новых ядер, но и на надёжной диагностике, неасимптотических гарантиях и масштабировании к сложным высокоразмерным моделям.

Основная идея

Пусть требуется вычислить интеграл

$I = E_{\pi} h(X)=\int_{\cal X}h(x)\,\pi(dx),$

где $\pi$ — целевое распределение на пространстве состояний ${\cal X}$ , а $h$ — интегрируемая функция. Если независимую выборку из $\pi$ получить трудно, строят цепь

$X_0,X_1,\ldots,X_N$

с переходным ядром $P(x,dy)$ , сохраняющим $\pi$ :

$\int_{\cal X}\pi(dx)P(x,A)=\pi(A).$

для любого измеримого множества $A$ . Это равенство означает, что если $X_t\sim\pi$ , то и $X_{t+1}\sim\pi$ .

При условиях эргодичности марковской цепи действует аналог закона больших чисел: при $N\longrightarrow\infty$ среднее сходится к $I$ почти наверное,

$\widehat I_N=\frac1N\sum_{t=1}^{N}h(X_t)\longrightarrow I.$

Тем самым задача интегрирования сводится к моделированию переходов цепи и усреднению. Нормирующая константа плотности часто не нужна. Если

$\pi(x)=\frac{\gamma(x)}{Z},\qquad Z=\int\gamma(x)\,dx,$

то отношения $\pi(y)/\pi(x)=\gamma(y)/\gamma(x)$ не содержат неизвестного $Z$ . Именно это свойство делает MCMC удобным для байесовских апостериорных распределений

$p(\theta\mid y)=\frac{p(y\mid\theta)p(\theta)}{p(y)}.$

Математические основы

Инвариантность и детальный баланс

Распространённым достаточным условием инвариантности является условие детального баланса:

$\pi(dx)P(x,dy)=\pi(dy)P(y,dx).$

Оно означает равенство стационарных потоков вероятности между любыми двумя областями пространства состояний. Цепь, удовлетворяющая этому условию, называется обратимой относительно $\pi$ . Обратимость удобна для конструирования и анализа алгоритмов, однако не является необходимой: существуют корректные необратимые MCMC-методы.

Для сходимости распределения $X_t$ к $\pi$ одной инвариантности недостаточно. В конечном пространстве состояний обычно требуют неприводимость и апериодичность. В общем пространстве соответствующую роль играют $\varphi$ -неприводимость, апериодичность и положительная харрисовская рекуррентность. Тогда при стандартных условиях

$\|P^t(x,\cdot)-\pi(\cdot)\|_{\rm TV}\longrightarrow 0,$

где $\|\cdot\|_{\rm TV}$ — Расстояние полной вариации. Условия геометрической эргодичности дают экспоненциальную оценку скорости сходимости и часто позволяют доказать центральную предельную теорему.^[1]

Автокорреляция и центральная предельная теорема

Пусть цепь стационарна, $Y_t=h(X_t)$ , $\gamma_k={\rm Cov}(Y_0,Y_k)$ , а $\rho_k=\gamma_k/\gamma_0$ . При выполнении центральной предельной теоремы для цепи величина $\sqrt N(\widehat I_N-I)$ сходится по распределению:

$\sqrt N\,(\widehat I_N-I)\longrightarrow {\cal N}(0,\sigma_h^2),\qquad \sigma_h^2=\gamma_0+2\sum_{k=1}^{\infty}\gamma_k.$

Интегрированное время автокорреляции и эффективный размер выборки определяют как

$\tau_{\rm int}=1+2\sum_{k=1}^{\infty}\rho_k,\qquad N_{\rm eff}=\frac{N}{\tau_{\rm int}}.$

Тогда асимптотическая дисперсия среднего приблизительно равна $\gamma_0/N_{\rm eff}$ . При положительной автокорреляции тысяча итераций может содержать информацию лишь нескольких десятков независимых наблюдений. При устойчивой отрицательной автокорреляции теоретически возможно $N_{\rm eff}>N$ ; поэтому эффективный размер выборки — мера точности конкретной оценки, а не буквально число независимых состояний.

Основные алгоритмы

Алгоритм Метрополиса — Хастингса

Алгоритм Метрополиса — Хастингса является базовой конструкцией MCMC. Его исходная симметричная версия была предложена Н. Метрополисом, А. Розенблют, М. Розенблют, А. Теллер и Э. Теллером для расчётов в статистической механике.^[1] У. Хастингс обобщил её на несимметричные распределения предложений.^[1]

Пусть текущему состоянию $x$ сопоставлено распределение предложения с плотностью $q(y\mid x)$ . Одна итерация состоит из следующих шагов:

Сгенерировать кандидата $Y\sim q(\cdot\mid x)$ .
Вычислить вероятность принятия
$\alpha(x,Y)=\min\left\{1,\frac{\gamma(Y)q(x\mid Y)}{\gamma(x)q(Y\mid x)}\right\}.$
С вероятностью $\alpha(x,Y)$ положить $X_{t+1}=Y$ , иначе сохранить $X_{t+1}=x$ .

Неизвестный нормирующий множитель целевой плотности сокращается в этом отношении. Самопереход при отклонении кандидата является частью марковского ядра и не должен удаляться из выборки. Полученный переход удовлетворяет детальному балансу. Если $q(y\mid x)=q(x\mid y)$ , поправка Хастингса сокращается и остаётся исходное правило Метрополиса.

Распространённый вариант — алгоритм случайного блуждания, где $Y=x+\varepsilon$ . Малый масштаб предложения даёт высокую долю принятия, но короткие перемещения; чрезмерно большой масштаб — частые отказы. Для последовательностей независимых одинаково распределённых компонент в пределе большой размерности оптимальный масштаб гауссовского случайного блуждания пропорционален $2{,}38/\sqrt d$ , а предельная доля принятия близка к $0{,}234$ .^[1] Это асимптотический ориентир, а не универсальная рекомендуемая доля принятия: для других целей, геометрий и типов предложений оптимальные значения иные.

Сэмплирование по Гиббсу

Сэмплирование по Гиббсу применимо, если можно генерировать значения из полных условных распределений. Для вектора $x=(x_1,\ldots,x_d)$ последовательный вариант на каждой итерации обновляет координаты:

$x_j^{(t+1)}\sim\pi\left(x_j\mid x_1^{(t+1)},\ldots,x_{j-1}^{(t+1)},x_{j+1}^{(t)},\ldots,x_d^{(t)}\right).$

Каждое такое обновление можно рассматривать как шаг Метрополиса — Хастингса с вероятностью принятия, равной единице. Координаты можно обновлять по одной или блоками; блочное обновление часто уменьшает зависимость, если компоненты сильно коррелированы. Метод особенно удобен для иерархических и сопряжённых байесовских моделей, но может перемешиваться крайне медленно при сильной апостериорной зависимости или почти детерминированной связи переменных.

Связь гиббсовского сэмплирования с гиббсовскими распределениями была использована С. и Д. Геманами в восстановлении изображений.^[1] Работа А. Гелфанда и А. Смита способствовала широкому распространению метода в байесовской статистике.^[1]

Сэмплирование сечений

Сэмплирование сечений (англ. slice sampling) вводит вспомогательную переменную $u$ и рассматривает область под графиком ненормированной плотности:

$0<u<\gamma(x).$

Попеременно генерируя $u\mid x$ и $x\mid u$ , получают марковскую цепь с предельной плотностью, пропорциональной $\gamma(x)$ . В одномерном случае метод автоматически подстраивает локальный масштаб перемещения, однако поиск и сужение интервала сечения требуют дополнительных вычислений плотности.^[1]

Градиентные методы

Если доступен градиент $\nabla\log\pi(x)$ , предложение можно направить в область большей плотности. В алгоритме Ланжевена с поправкой Метрополиса (MALA) используют, например,

$Y\sim {\cal N}\left(x+\frac{\epsilon^2}{2}M\nabla\log\pi(x),\,\epsilon^2M\right)$

и затем применяют поправку Метрополиса — Хастингса. Матрица $M$ задаёт масштаб и направление. Без шага принятия дискретизированная ланжевеновская динамика обычно имеет смещённое стационарное распределение; поправка устраняет ошибку дискретизации, но не гарантирует хорошего перемешивания для произвольной плотности.

Гамильтоновский метод Монте-Карло (HMC) расширяет пространство переменной положения $q$ искусственным импульсом $p$ и задаёт гамильтониан

$H(q,p)=U(q)+K(p),\qquad U(q)=-\log\gamma(q),\qquad K(p)=\frac12p^T M^{-1}p.$

Численное интегрирование уравнений Гамильтона обратимым сохраняющим объём методом, обычно схемой «чехарда», создаёт дальнее предложение без поведения обычного случайного блуждания. Ошибку численного интегрирования исправляет шаг Метрополиса с вероятностью $\min\{1,\exp(-\Delta H)\}$ . Метод был предложен в вычислительной физике под названием «гибридный Монте-Карло».^[1]

Эффективность HMC зависит от шага интегрирования, длины траектории и матрицы масс. Алгоритм NUTS (англ. No-U-Turn Sampler) автоматически прекращает наращивание траектории, когда она начинает разворачиваться, и адаптирует шаг во время настройки.^[1] HMC и NUTS хорошо работают для многих гладких непрерывных моделей высокой размерности, однако непосредственно не обновляют дискретные параметры и чувствительны к сильно меняющемуся масштабу распределения.

Методы для специальных задач

Обратимые скачки (RJMCMC) позволяют переходить между моделями различной размерности; в коэффициент принятия входит якобиан преобразования.^[1]
Псевдомаргинальный MCMC подставляет в отношение Метрополиса — Хастингса неотрицательную несмещённую оценку правдоподобия. Цепь работает в расширенном пространстве, а её маргинальное стационарное распределение остаётся точным целевым распределением; простая подстановка произвольной шумной оценки этим свойством не обладает.^[1]
Параллельное умеривание и обмен репликами моделируют несколько распределений при разных «температурах» и обменивают состояния цепей. Это облегчает переходы между удалёнными модами, но требует выбора температурной сетки и дополнительных вычислений.
Адаптивный MCMC настраивает масштаб, ковариацию или другие параметры предложения по истории цепи. Поскольку переходное ядро при этом меняется, обычная теория однородных цепей неприменима автоматически; для корректности нужны дополнительные условия, например затухающая адаптация и условие удержания.^[1]
Стохастические градиентные методы, в частности SGLD, оценивают градиент по мини-пакету данных. При убывающем шаге SGLD имеет асимптотическое обоснование, а практически используемый постоянный шаг вносит ошибку дискретизации; поэтому не всякий алгоритм, называемый стохастическим градиентным MCMC, является точным при конечном времени вычисления.^[1]

Диагностика и оценка точности

Конечная траектория не сообщает напрямую, достигла ли цепь стационарного режима. Универсальной диагностики, способной доказать сходимость по одному конечному запуску, не существует: согласие диагностик является необходимой практической проверкой, но не математическим доказательством.

Начальная настройка и несколько цепей

Начальные итерации зависят от распределения $X_0$ . Их часто называют периодом прогрева или «выжигания» (англ. warm-up, burn-in) и исключают из итоговых оценок. Однако отбрасывание фиксированной доли траектории само по себе не гарантирует исчезновения начального смещения. Предпочтительно запускать несколько цепей из разумно рассеянных начальных точек, отделять адаптацию от последующего сэмплирования и проверять, исследуют ли цепи одни и те же области.

Классическая диагностика Гелмана — Рубина сравнивает внутрисерийную и межсерийную дисперсии и выражается потенциальным коэффициентом уменьшения масштаба $\widehat R$ .^[1] Современная версия разбивает цепи, применяет ранговую нормализацию и складывание относительно медианы. Она устойчивее при тяжёлых хвостах и различии масштабов; её авторы рекомендуют добиваться $\widehat R<1{,}01$ для всех интересующих скалярных величин.^[1] Значение, близкое к единице, не исключает общего пропуска всеми цепями изолированной моды.

Графические проверки

К стандартным графическим средствам относятся:

график траектории — выявляет дрейф, застревание, редкие переходы и различия между цепями;
ранговый график — показывает, одинаково ли цепи покрывают распределение;
автокорреляционная функция — характеризует зависимость между итерациями;
графики пар параметров и энергии — выявляют сильную зависимость и патологическую геометрию;
для HMC — число расходящихся траекторий, достижение максимальной глубины дерева и статистика энергии.

Гладкий график маргинальной плотности может быть обманчив: цепь, застрявшая в одной моде, способна давать визуально устойчивую оценку внутри этой моды.

Эффективный размер выборки и ошибка Монте-Карло

Эффективный размер следует оценивать не только для параметров, но и для функций, которые входят в научные выводы: вероятностей событий, квантилей и предсказаний. Ранговые bulk-ESS и tail-ESS отдельно оценивают точность в центре и хвостах распределения.^[1] Для векторных функционалов разработан многомерный ESS и правила остановки, основанные на оценке ковариации в многомерной центральной предельной теореме.^[1]

Оценка стандартной ошибки Монте-Карло имеет вид

${\rm MCSE}(\widehat I_N)=\sqrt{\widehat\sigma_h^2/N},$

где $\widehat\sigma_h^2$ оценивают спектральными методами, пакетными средними или методами начальной положительной последовательности. Длину вычисления следует выбирать по требуемой MCSE, а не по фиксированному числу итераций.

Прореживание (сохранение каждого $k$ -го состояния) обычно не увеличивает статистическую информацию при фиксированном числе вычисленных переходов: отброшенные состояния всё ещё содержат информацию. Оно оправдано главным образом ограничениями памяти, стоимостью последующей обработки или необходимостью записывать менее плотную траекторию.^[1]

Практическое построение модели и алгоритма

Надёжное MCMC-исследование обычно включает следующие этапы:

Определение целевого распределения. Проверяют корректность правдоподобия, априорных распределений, ограничений и якобианов преобразований. Несобственное апостериорное распределение нельзя исправить более длинной цепью.
Выбор параметризации. Масштабирование переменных, удаление неидентифицируемости и переход между центрированной и нецентрированной параметризациями могут изменять эффективность на порядки. Для ограниченных параметров часто используют взаимно однозначное преобразование в неограниченное пространство.
Выбор ядра. Для гладких непрерывных моделей часто используют HMC/NUTS; для доступных полных условных распределений — гиббсовские обновления; для дискретных или трансразмерных задач — специальные шаги Метрополиса, блочные или RJMCMC-переходы.
Настройка. Адаптируют масштаб, ковариацию предложения, шаг интегрирования или матрицу масс только по правилам, сохраняющим теоретическую корректность. Результаты периода адаптации обычно не включают в оценки.
Несколько независимых запусков. Начальные точки выбирают так, чтобы проверить разные правдоподобные области пространства, но избегают заведомо недопустимых состояний.
Диагностика. Совместно анализируют $\widehat R$ , bulk-ESS, tail-ESS, MCSE, графики рангов и траекторий, а также специфические предупреждения алгоритма.
Проверка модели. Сходимость вычислений не означает адекватности статистической модели. Для байесовских моделей отдельно выполняют априорные и апостериорные предиктивные проверки и анализ чувствительности.
Воспроизводимость. Фиксируют версию программы, генератор и начальное значение случайных чисел, параметризацию, число цепей, длительность прогрева и сэмплирования, критерии остановки и все предупреждения.

Трудности и ограничения

Медленное перемешивание и мультимодальность

Если моды разделены областью крайне малой вероятности, локальная цепь может практически не переходить между ними. Тогда эмпирические средние выглядят стабильными, хотя отражают лишь одну часть распределения. Возможные средства — репараметризация, блочные и независимые предложения, умеривание, последовательные методы Монте-Карло или специализированные глобальные переходы. Ни один из этих подходов не устраняет проблему автоматически.

Высокая размерность и геометрия

У случайного блуждания типичная длина допустимого шага уменьшается с размерностью, что приводит к диффузионному исследованию пространства. Градиентные методы используют локальную геометрию, но сталкиваются с «воронками», сильной кривизной и масштабами, различающимися на много порядков. В HMC такие особенности проявляются расходящимися траекториями или необходимостью очень малого шага. Часто эффективнее изменить параметризацию модели, чем просто увеличить число итераций.

Вычислительная стоимость и большие данные

В байесовской модели с независимыми наблюдениями полная оценка логарифма правдоподобия и его градиента обычно стоит $O(n)$ на шаг. Подвыборка данных уменьшает стоимость, но наивная замена полной суммы мини-пакетной оценкой, вообще говоря, меняет стационарное распределение. Точные псевдомаргинальные и событийные методы требуют специальных несмещённых оценок, ограничений интенсивности или контрольных вариат.

Дискретные параметры и границы

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

Современные направления исследований

Исследования MCMC развиваются в нескольких направлениях.

Несмещённые оценки за конечное время. Обычное MCMC-среднее имеет начальное смещение при любом фиксированном числе шагов. Сцепление двух цепей и телескопические суммы позволяют строить несмещённые оценки, которые можно независимо вычислять параллельно, хотя их эффективность зависит от времени встречи цепей.^[1]
Необратимые событийные процессы. Кусочно-детерминированные марковские процессы, такие как Zig-Zag, заменяют диффузионное случайное блуждание направленным движением со случайными событиями. Для некоторых моделей возможны точные схемы подвыборки без ошибки дискретизации.^[1]
Локально сбалансированные предложения. Предложения Баркера используют градиент иначе, чем MALA, и могут быть устойчивее к ошибочному масштабу и неоднородной геометрии; теоретические и эмпирические преимущества зависят от класса целевых распределений.^[1]
Локальная адаптация HMC. Опубликованный в 2026 году WALNUTS расширяет NUTS, меняя шаг схемы «чехарда» внутри одной траектории по контролю ошибки энергии. В экспериментах авторов метод был устойчивее стандартного NUTS на многомасштабных распределениях; это новый результат, а не универсально установленное превосходство для всех задач.^[1]

Одновременно развивается теория неасимптотических оценок времени смешивания, автоматического выбора параметризации, параллельных взаимодействующих цепей и сочетания MCMC с обучаемыми транспортными отображениями. Для таких методов особенно важно различать точные алгоритмы, сохраняющие заданное распределение, и приближённые алгоритмы с контролируемым либо неконтролируемым смещением.

Применения

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

См. также

Примечания

Литература

Brooks, S.; Gelman, A.; Jones, G. L.; Meng, X.-L. Handbook of Markov Chain Monte Carlo. — Boca Raton: Chapman & Hall/CRC, 2011. — ISBN 978-1-4200-7941-8
Gelman, A.; Carlin, J. B.; Stern, H. S.; Dunson, D. B.; Vehtari, A.; Rubin, D. B. Bayesian Data Analysis. — 3-е изд.. — Boca Raton: Chapman & Hall/CRC, 2013. — ISBN 978-1-4398-4095-5
Meyn, S. P.; Tweedie, R. L. Markov Chains and Stochastic Stability. — 2-е изд.. — Cambridge: Cambridge University Press, 2009. — ISBN 978-0-521-73182-9
Robert, C. P.; Casella, G. Monte Carlo Statistical Methods. — 2-е изд.. — New York: Springer, 2004. — ISBN 978-0-387-21239-5

Ссылки

Источник — «http://www.machinelearning.ru/wiki/index.php?title=%D0%9C%D0%B5%D1%82%D0%BE%D0%B4_%D0%9C%D0%BE%D0%BD%D1%82%D0%B5-%D0%9A%D0%B0%D1%80%D0%BB%D0%BE_%D0%BF%D0%BE_%D1%81%D1%85%D0%B5%D0%BC%D0%B5_%D0%BC%D0%B0%D1%80%D0%BA%D0%BE%D0%B2%D1%81%D0%BA%D0%B8%D1%85_%D1%86%D0%B5%D0%BF%D0%B5%D0%B9»

Категории: Методы Монте-Карло | Байесовское обучение | Энциклопедия анализа данных

@@ Строка 1: / Строка 1: @@
-{{well|Статья подготовлена с использованием модели [https://openai.com/index/gpt-5-6/ OpenAI GPT‑5.6 Sol] с уровнем рассуждений High и проверена участником [[Участник:Denis Kistanov|Д.О. Кистанов]] 16:57, 18 июля 2026 (MSD)
+{{well|Статья подготовлена с использованием модели [https://openai.com/index/gpt-5-6/ OpenAI GPT‑5.6 Sol] с уровнем рассуждений High и проверена участником [[Участник:Denis Kistanov|Д.О. Кистанов]] 16:57, 18 июля 2026 (MSK)
 Промпт приводится полностью в [[Обсуждение:Метод Монте-Карло по схеме марковских цепей]]
 }}