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

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

Статья написана с использованием LLM Gemini Pro 3.1 и проверена участником Nikita Zinoviсh 16:50, 16 июля 2026 (MSD)

Содержание

1 Введение и мотивация
2 Теоретический фундамент: цепи Маркова
3 Строгая постановка задачи
4 Описание алгоритма Метрополиса — Гастингса
5 Математическое доказательство корректности
6 Важные частные случаи
7 См. также
8 Литература

Введение и мотивация

В задачах байесовского вывода, вычислительной статистики и статистической физики фундаментальной проблемой является генерация выборок из сложного многомерного вероятностного распределения. Пусть на измеримом пространстве состояний $X$ задано целевое распределение вероятностей с плотностью $\pi(x)$ . Как правило, эта плотность известна аналитически лишь с точностью до нормировочной константы:

$\pi(x) = \frac{p(x)}{Z}$

где $p(x) > 0$ — вычислимая ненормированная функция (например, произведение априорной плотности и функции правдоподобия), а $Z = \int_X p(x) dx$ — нормировочная константа (статистическая сумма или маргинальное правдоподобие). Для многомерных пространств прямое или численное интегрирование для нахождения $Z$ вычислительно неразрешимо ввиду проклятия размерности.

Алгоритм Метрополиса — Гастингса представляет собой универсальный метод из семейства MCMC (Markov Chain Monte Carlo), позволяющий строить эргодическую марковскую цепь, стационарное распределение которой сходится к $\pi(x)$ . При этом алгоритм оперирует исключительно отношениями плотностей, что делает вычисление многомерного интеграла $Z$ излишним.

Теоретический фундамент: цепи Маркова

Пусть $\{x_t\}_{t=0}^\infty$ — однородная марковская цепь на непрерывном пространстве $X$ . Переходы между состояниями полностью определяются переходным ядром $P(y \mid x)$ , которое описывает условную плотность вероятности перехода системы в состояние $y$ при условии, что текущим состоянием является $x$ .

Распределение $\pi(x)$ называется стационарным (или инвариантным) для данного переходного ядра, если оно сохраняется при действии интегрального оператора перехода:

$\pi(y) = \int_X \pi(x) P(y \mid x) dx$

Достаточным условием стационарности распределения является выполнение уравнения детального равновесия (условия обратимости марковской цепи):

$\pi(x) P(y \mid x) = \pi(y) P(x \mid y) \quad \forall x, y \in X$

Интегрируя обе части этого уравнения по $x$ , нетрудно убедиться, что из детального равновесия автоматически следует условие стационарности. Таким образом, глобальная задача MCMC сводится к синтезу такого марковского ядра $P(y \mid x)$ , для которого требуемое целевое распределение $\pi(x)$ удовлетворяло бы условию детального равновесия.

Строгая постановка задачи

Дано:

Пространство состояний $X \subseteq \mathbb{R}^d$ .
Целевая плотность распределения $\pi(x) \propto p(x)$ , где функция $p(x)$ алгоритмически вычислима в любой точке области определения.

Требуется: Построить конструктивный алгоритм задания переходного ядра $P(y \mid x)$ , строго удовлетворяющего уравнению детального равновесия по отношению к $\pi(x)$ . В силу эргодической теоремы для марковских цепей (при условии выполнения требований неприводимости и апериодичности), распределение элементов сгенерированной последовательности при $t \to \infty$ будет слабо сходиться к целевому распределению $\pi(x)$ , а выборочные средние будут состоятельными оценками математических ожиданий.

Описание алгоритма Метрополиса — Гастингса

Для реализации перехода из текущего состояния $x$ алгоритм использует вспомогательное распределение предложений (proposal distribution) с плотностью $q(y \mid x)$ . Данное распределение отвечает за стохастическую генерацию возможных состояний-кандидатов.

Для компенсации возможной асимметрии распределения $q$ и обеспечения строгого выполнения детального равновесия вводится вероятность принятия (acceptance probability) предложенного шага $\alpha(x, y)$ :

$\alpha(x, y) = \min \left( 1, \frac{\pi(y) q(x \mid y)}{\pi(x) q(y \mid x)} \right) = \min \left( 1, \frac{p(y) q(x \mid y)}{p(x) q(y \mid x)} \right)$

Поскольку плотности $\pi$ входят в формулу в виде отношения, невычислимая нормировочная константа $Z$ тождественно сокращается.

Формальная итерация алгоритма: Пусть в дискретный момент времени $t$ марковская цепь находится в состоянии $x_t = x$ . 1. Сгенерировать состояние-кандидат $y$ из распределения предложений: $y \sim q(\cdot \mid x)$ . 2. Вычислить вероятность принятия $\alpha(x, y)$ . 3. Сгенерировать случайную величину $u$ из стандартного равномерного распределения $U(0, 1)$ . 4. Принять решение о переходе:

- Если $u \leq \alpha(x, y)$ , переход принимается: $x_{t+1} = y$ .
- В противном случае переход отклоняется, и цепь остаётся в прежнем состоянии: $x_{t+1} = x$ .

Математическое доказательство корректности

Для математически строгого доказательства необходимо выписать полное выражение для переходного ядра генерируемой цепи. Ядро $P(y \mid x)$ является композицией непрерывной компоненты (в случае успешного перехода $x \neq y$ ) и сингулярной компоненты (вероятность остаться в точке $x$ при отклонении кандидата).

Вероятность того, что кандидат будет отклонён (и система останется в $x$ ), выражается как интеграл:

$r(x) = 1 - \int_X q(y \mid x) \alpha(x, y) dy$

С учётом этого, полное переходное ядро $P(y \mid x)$ выражается через дельта-функцию Дирака:

$P(y \mid x) = q(y \mid x) \alpha(x, y) + r(x) \delta(y - x)$

Проверим выполнение условия детального равновесия: $\pi(x) P(y \mid x) = \pi(y) P(x \mid y)$ .

Случай 1: $x = y$ . Уравнение обратимости принимает вид $\pi(x) P(x \mid x) = \pi(x) P(x \mid x)$ и выполняется тривиально.

Случай 2: $x \neq y$ . Сингулярная часть ядра обращается в нуль, и задача сводится к доказательству тождества для непрерывной компоненты:

$\pi(x) q(y \mid x) \alpha(x, y) = \pi(y) q(x \mid y) \alpha(y, x)$

Подставим явный вид функции $\alpha(x, y)$ в левую часть:

$\pi(x) q(y \mid x) \min \left( 1, \frac{\pi(y) q(x \mid y)}{\pi(x) q(y \mid x)} \right)$

Поскольку величины $\pi(x)$ и $q(y \mid x)$ строго неотрицательны, их можно внести под знак функции минимума:

$\min \big( \pi(x) q(y \mid x), \, \pi(y) q(x \mid y) \big)$

Проделаем аналогичные алгебраические преобразования с правой частью уравнения $\pi(y) q(x \mid y) \alpha(y, x)$ :

$\pi(y) q(x \mid y) \min \left( 1, \frac{\pi(x) q(y \mid x)}{\pi(y) q(x \mid y)} \right) = \min \big( \pi(y) q(x \mid y), \, \pi(x) q(y \mid x) \big)$

Полученные выражения абсолютно идентичны в силу коммутативности функции минимума. Следовательно, переходное ядро $P(y \mid x)$ строго удовлетворяет уравнению детального равновесия. Это завершает доказательство того факта, что $\pi(x)$ является инвариантной мерой для сконструированной цепи Маркова.

Важные частные случаи

Классический алгоритм Метрополиса

В случае, когда распределение предложений является симметричным, то есть выполняется равенство $q(y \mid x) = q(x \mid y)$ для любых $x, y \in X$ , множитель $q$ в отношении сокращается. Формула вероятности принятия принимает свой канонический вид:

$\alpha(x, y) = \min \left( 1, \frac{p(y)}{p(x)} \right)$

Физический (и вероятностный) смысл этого условия прозрачен: если предложенное состояние $y$ обладает большей или равной целевой плотностью по сравнению с текущим, такой переход принимается достоверно ( $\alpha = 1$ ). Менее вероятные состояния принимаются с вероятностью, равной отношению их плотностей.

Алгоритм случайного блуждания (Random Walk Metropolis)

Наиболее распространённый на практике подход заключается в генерации кандидата путём добавления случайного центрированного возмущения: $y = x + \varepsilon$ , где вектор $\varepsilon \sim g(\varepsilon)$ , а плотность шума $g$ симметрична относительно начала координат. Стандартным выбором является Многомерное нормальное распределение: $q(y \mid x) = \mathcal{N}(y \mid x, \Sigma)$ . Поскольку $q(y \mid x) = g(y - x) = g(x - y) = q(x \mid y)$ , данный метод является частным случаем классического алгоритму Метрополиса. Фундаментальной задачей при его использовании является настройка ковариационной матрицы $\Sigma$ для обеспечения оптимальной скорости сходимости цепи.

Независимый сэмплер (Independence Sampler)

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

$\alpha(x, y) = \min \left( 1, \frac{p(y) q(x)}{p(x) q(y)} \right) = \min \left( 1, \frac{w(y)}{w(x)} \right)$

где $w(x) = \frac{p(x)}{q(x)}$ — весовая функция. Эффективность независимого сэмплера критически зависит от того, насколько точно априорно заданное распределение $q(x)$ аппроксимирует целевую функцию $p(x)$ (особенно в областях с «тяжёлыми хвостами»).

См. также

Литература

Hastings, W. K. Monte Carlo sampling methods using Markov chains and their applications // Biometrika. — 1970. — Vol. 57, No. 1. — P. 97–109.
Robert, C. P., Casella, G. Monte Carlo Statistical Methods. — Springer, 2004. — 645 p. (Фундаментальная монография со строгим математическим обоснованием алгоритма).
Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., Rubin, D. B. Bayesian Data Analysis. — CRC Press, 2013. — 675 p.

Источник — «http://machinelearning.ru/wiki/index.php?title=%D0%90%D0%BB%D0%B3%D0%BE%D1%80%D0%B8%D1%82%D0%BC_%D0%9C%D0%B5%D1%82%D1%80%D0%BE%D0%BF%D0%BE%D0%BB%D0%B8%D1%81%D0%B0-%D0%93%D0%B0%D1%81%D1%82%D0%B8%D0%BD%D0%B3%D1%81%D0%B0»