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

Алгоритм Метрополиса — Гастингса — алгоритм семплирования, использующийся, в основном, для сложных функций распределения. Он отчасти похож на алгоритм выборки с отклонением, однако здесь вспомогательная функция распределения меняется со временем. Алгоритм был впервые опубликован Николасом Метрополисом в 1953 году, и затем обобщён Шаблон:Nobr в 1970 году. Семплирование по Гиббсу является частным случаем алгоритма Метрополиса — Гастингса и более популярно за счёт простоты и скорости, хотя и реже применимо.

Алгоритм Метрополиса — Гастингса позволяет семплировать любую функцию распределения. Он основан на создании цепи Маркова, то есть на каждом шаге алгоритма новое выбранное значение $x^{t+1}$ зависит только от предыдущего $x^t$. Алгоритм использует вспомогательную функцию распределения $Q(x^{\prime }|x^t)$, зависящую от $x^t$, для которой генерировать выборку просто (например, нормальное распределение). На каждом шаге для этой функции генерируется случайное значение $x^{\prime }$. Затем с вероятностью

$u=\frac{P(x^{\prime })Q(x^t|x^{\prime })}{P(x^t)Q(x^{\prime }|x^t)}$

(или с вероятностью 1, если $u>1$), выбранное значение принимается как новое: $x^{t+1}=x^{\prime }$, а иначе оставляется старое: $x^{t+1}=x^t$.

Например, если взять нормальную функцию распределения как вспомогательную функцию, то

$Q(x^{\prime }|x^t)\sim N(x^t,{\sigma }^{2}I).$

Такая функция выдаёт новое значение в зависимости от значения на предыдущем шаге. Изначально алгоритм Метрополиса требовал, чтобы вспомогательная функция была симметрична: $Q(x^{\prime },x^t)=Q(x^t,x^{\prime })$, однако обобщение Гастингса снимает это ограничение.

Алгоритм

Пусть мы уже выбрали случайное значение $x^t$. Для выбора следующего значения сначала получим случайное значение $x^{\prime }$ для функции $Q(x^{\prime }|x^t)$. Затем найдем произведение $a=a_1{a}_2$, где

$a_1=\frac{P(x^{\prime })}{P(x^t)}$

является отношением вероятностей между промежуточным значением и предыдущим, а

$a_2=\frac{Q(x^t|x^{\prime })}{Q(x^{\prime }|x^t)}$

это отношение между вероятностями пойти из $x^{\prime }$ в $x^t$ или обратно. Если $Q$ симметрична, то второй множитель равен 1. Случайное значение на новом шаге выбирается по правилу:

$\begin{array}{cc}\text{If }a\ge 1:& \\ & x^{t+1}=x^{\prime },\end{array}$

$\begin{array}{cc}\text{and if }a<1:& \\ & x^{t+1}=\{\begin{array}{c}{x}^{\prime }\text{ with probability }a\\ x^t\text{ with probability }1-a.\end{array}\end{array}$

Алгоритм стартует из случайного значения $x^0$, и сначала прогоняется «вхолостую» некоторое количество шагов, чтобы «забыть» о начальном значении.

Лучше всего алгоритм работает тогда, когда форма вспомогательной функции близка к форме целевой функции $P$. Однако добиться этого априори зачастую невозможно. Для решения этой проблемы вспомогательную функцию настраивают в ходе подготовительной стадии работы алгоритма. Например, для нормального распределения настраивают его параметр ${\sigma }^2$ так, чтобы доля «принятых» случайных значений (то есть тех, для которых $x^{t+1}=x^{\prime }$) была близка к 60 %. Если ${\sigma }^2$ слишком мала, то значения будут получаться слишком близкими и доля принятых будет высока. Если ${\sigma }^2$ слишком велика, то с большой вероятностью новые значения будут выскакивать в зоны малой вероятности $P$, отчего доля принятых значений окажется низкой.