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

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

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

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

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

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

$Q(x^{\prime }|x^t)\sim <!--\; \sim \; -->N(x^t,{\mathrm{\sigma <!--\; \sigma \; -->}}^{2}I).$$Q(x'|x^{t})\sim N(x^{t},\sigma ^{2}I).$

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

АлгоритмПравить

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

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

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

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

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

${\begin{matrix}{\mbox{If }}a\geq 1:&\\&x^{t+1}=x',\end{matrix}}$ 

${\begin{matrix}{\mbox{and if }}a<1:&\\&x^{t+1}=\left\{{\begin{matrix}x'{\mbox{ with probability }}a\\x^{t}{\mbox{ with probability }}1-a.\end{matrix}}\right.\end{matrix}}$ 

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

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