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

Алгоритм Метрополиса-Гастингса производит выборку нормального одномерного апостериорного распределения вероятностей.

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

Метрополис-Гастингс и другие алгоритмы MCMC обычно используются для выборки из многомерных распределений, особенно когда количество измерений велико. Для одномерных распределений обычно существуют другие методы (например, адаптивная бракованная выборка ), которые могут напрямую возвращать независимые выборки из распределения, и они свободны от проблемы автокоррелированных выборок, присущей методам MCMC.

История

Алгоритм частично назван в честь Николаса Метрополиса , первого соавтора статьи 1953 года под названием «Уравнение вычислений состояния с помощью быстрых вычислительных машин» , совместно с Арианной В. Розенблут , Маршаллом Розенблутом , Огастой Х. Теллер и Эдвардом Теллером . В течение многих лет алгоритм был известен просто как алгоритм Метрополиса . В статье был предложен алгоритм для случая симметричного распределения предложений, но в 1970 году У.К. Гастингс распространил его на более общий случай. Обобщенный метод в конечном итоге получил оба названия, хотя первое использование термина «алгоритм Метрополиса-Гастингса» неясно.

Некоторые разногласия существуют относительно заслуги в разработке алгоритма Метрополиса. Метрополис, который был знаком с вычислительными аспектами метода, ввел термин «Монте-Карло» в более ранней статье со Станиславом Уламом и возглавил группу теоретического отдела, которая спроектировала и построила компьютер MANIAC I, использованный в экспериментах в 1952. Однако до 2003 года подробного описания разработки алгоритма не было. Незадолго до своей смерти Маршалл Розенблют присутствовал на конференции 2003 года в LANL, посвященной 50-летию публикации 1953 года. На этой конференции Розенблут описал алгоритм и его развитие в презентации под названием «Происхождение алгоритма Монте-Карло для статистической механики». Дальнейшие исторические разъяснения сделаны Губернатисом в журнальной статье 2005 года. рассказывая о 50-летии конференции. Розенблут ясно дает понять, что работу выполнили он и его жена Арианна, и что «Метрополис» не играл никакой роли в разработке, кроме предоставления компьютерного времени.

Это противоречит сообщению Эдварда Теллера, который утверждает в своих мемуарах, что пять авторов статьи 1953 года работали вместе «днями (и ночами)». Напротив, в подробном отчете Розенблюта Теллеру приписывают важное, но раннее предложение «воспользоваться преимуществами статистической механики и взять средние значения по ансамблю вместо того, чтобы следовать подробной кинематике ». Это, по словам Розенблюта, заставило его задуматься об обобщенном подходе Монте-Карло – теме, которую, по его словам, он часто обсуждал с Джоном фон Нейманом . Арианна Розенблут рассказала (Губернатису в 2003 году), что Огаста Теллер начала работу с компьютером, но Арианна сама взяла на себя управление и написала код с нуля. В устной истории, записанной незадолго до его смерти, Розенблут снова приписывает Теллеру постановку исходной проблемы, ему самому ее решение, а Арианне - программирование компьютера.

Описание

Алгоритм Метрополиса – Гастингса может брать выборки из любого распределения вероятностей с плотностью вероятности. $P(x)$, при условии, что мы знаем функцию $f(x)$ пропорционален плотности $P$ и ценности $f(x)$ можно рассчитать. Требование, чтобы $f(x)$ должен быть только пропорционален плотности, а не точно равен ей, что делает алгоритм Метрополиса-Гастингса особенно полезным, поскольку он устраняет необходимость расчета коэффициента нормализации плотности, что на практике часто бывает чрезвычайно сложно.

Алгоритм Метрополиса – Гастингса генерирует последовательность значений выборки таким образом, что по мере создания все большего и большего количества значений выборки распределение значений более точно приближается к желаемому распределению. Эти значения выборки создаются итеративно таким образом, что распределение следующей выборки зависит только от текущего значения выборки, что превращает последовательность выборок в цепь Маркова . В частности, на каждой итерации алгоритм предлагает кандидата на следующее значение выборки на основе текущего значения выборки. Затем с некоторой вероятностью кандидат либо принимается, и в этом случае значение кандидата используется на следующей итерации, либо он отклоняется, и в этом случае значение кандидата отбрасывается, а текущее значение повторно используется на следующей итерации. Вероятность приемки определяется путем сравнения значений функции $f(x)$ текущих и потенциальных значений выборки относительно желаемого распределения.

Метод, используемый для предложения новых кандидатов, характеризуется распределением вероятностей $g(x\mid y)$ (иногда пишется $Q(x\mid y)$) нового предложенного образца $x$ учитывая предыдущий образец $y$. Это называется плотностью предложения , функцией предложения или прыгающим распределением . Распространенный выбор для $g(x\mid y)$ представляет собой гауссово распределение с центром в $y$, так что это указывает ближе к $y$ с большей вероятностью будут посещены следующими, превращая последовательность выборок в гауссово случайное блуждание . В оригинальной статье Metropolis et al. (1953), $g(x\mid y)$ предполагалось, что это равномерное распределение, ограниченное некоторым максимальным расстоянием от $y$. Также возможны более сложные функции предложения, такие как гамильтониан Монте-Карло , Ланжевена-Монте-Карло или предварительно обусловленный Кранк-Николсон .

В целях иллюстрации ниже описан алгоритм Метрополиса, частный случай алгоритма Метрополиса – Гастингса, в котором функция предложения симметрична.

Алгоритм Метрополиса (симметричное распределение предложений)

Позволять $f(x)$ быть функцией, которая пропорциональна желаемой функции плотности вероятности $P(x)$ (также известное как целевое распределение) .

1. Инициализация: выберите произвольную точку $x_{t}$ быть первым наблюдением в выборке и выбрать функцию предложения $g(x\mid y)$. В этом разделе $g$ предполагается симметричным; другими словами, оно должно удовлетворять $g(x\mid y)=g(y\mid x)$.
2. Для каждой итерации t :
   • Предложить кандидата $x'$ для следующей выборки, выбрав из распределения $g(x'\mid x_{t})$.
   • Рассчитать коэффициент приемки $\alpha =f(x')/f(x_{t})$, который будет использоваться для принятия решения о принятии или отклонении кандидата . Поскольку f пропорциональна плотности P , мы имеем, что $\alpha =f(x')/f(x_{t})=P(x')/P(x_{t})$.
   • Принять или отклонить :
     • Генерация равномерного случайного числа $u\in [0,1]$.
     • Если $u\leq \alpha$, затем примите кандидата, установив $x_{t+1}=x'$,
     • Если $u>\alpha$, затем отклоните кандидата и установите $x_{t+1}=x_{t}$ вместо.

Этот алгоритм действует путем случайных попыток перемещения по выборочному пространству, иногда соглашаясь на ходы, а иногда оставаясь на месте. Обратите внимание, что приемочный коэффициент $\alpha$ указывает, насколько вероятен новый предложенный образец по отношению к текущему образцу в соответствии с распределением, плотность которого равна $P(x)$. Если мы попытаемся перейти к точке, которая более вероятна, чем существующая точка (т. е. точка в области с более высокой плотностью $P(x)$ соответствующий $\alpha >1\geq u$), мы всегда примем этот переезд. Однако если мы попытаемся перейти к менее вероятной точке, мы иногда отклоним этот шаг, и чем больше относительное падение вероятности, тем больше вероятность того, что мы отклоним новую точку. Таким образом, мы будем стремиться оставаться в регионах с высокой плотностью населения (и возвращать большое количество образцов из них). $P(x)$, при этом лишь изредка посещая регионы с низкой плотностью населения. Интуитивно понятно, почему этот алгоритм работает и возвращает образцы, которые соответствуют желаемому распределению с плотностью. $P(x)$.

По сравнению с таким алгоритмом, как адаптивная отбраковочная выборка который напрямую генерирует независимые выборки из распределения, алгоритмы Метрополиса – Гастингса и другие MCMC имеют ряд недостатков:

• Выборки автокоррелированы . Даже несмотря на то, что в долгосрочной перспективе они правильно следуют $P(x)$, набор соседних выборок будет коррелировать друг с другом и неправильно отражать распределение. Это означает, что эффективный размер выборки может быть значительно меньше, чем количество фактически взятых выборок, что приводит к большим ошибкам.
• Хотя цепь Маркова в конечном итоге сходится к желаемому распределению, исходные образцы могут иметь совсем другое распределение, особенно если отправная точка находится в области низкой плотности. В результате приработки , обычно необходим период где первоначальное количество образцов выбрасывается.

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

В многомерных распределениях классический алгоритм Метрополиса – Гастингса, описанный выше, включает выбор новой многомерной точки выборки. Когда количество измерений велико, найти подходящее распределение прыжков для использования может быть сложно, поскольку разные отдельные измерения ведут себя очень по-разному, а ширина прыжка (см. выше) должна быть «правильной» для всех измерений одновременно, чтобы избегайте слишком медленного перемешивания. Альтернативный подход, который часто работает лучше в таких ситуациях, известный как выборка Гиббса , предполагает выбор новой выборки для каждого измерения отдельно от других, а не выбор выборки для всех измерений одновременно. Таким образом, проблема выборки из потенциально многомерного пространства будет сведена к набору задач по выборке из маломерного пространства. Это особенно применимо, когда многомерное распределение состоит из набора отдельных случайных величин , в которых каждая переменная обусловлена ​​лишь небольшим количеством других переменных, как это имеет место в большинстве типичных иерархических моделей . Затем отдельные переменные выбираются по одной, при этом каждая переменная зависит от самых последних значений всех остальных. Для выбора этих отдельных выборок могут использоваться различные алгоритмы, в зависимости от точной формы многомерного распределения: некоторые возможности - это методы адаптивной выборки с отбраковкой , алгоритм адаптивного отклонения Metropolis, простой одномерный шаг Метрополиса-Гастингса или срезовая выборка .

Формальный вывод

Цель алгоритма Метрополиса – Гастингса — создать набор состояний в соответствии с желаемым распределением. $P(x)$. Для этого алгоритм использует марковский процесс , который асимптотически достигает уникального стационарного распределения. $\pi (x)$ такой, что $\pi (x)=P(x)$.

Марковский процесс однозначно определяется своими вероятностями перехода. $P(x'\mid x)$, вероятность перехода из любого заданного состояния $x$ в любое другое данное состояние $x'$. Имеет уникальное стационарное распределение. $\pi (x)$ когда выполняются следующие два условия:

1. Существование стационарного распределения : должно существовать стационарное распределение. $\pi (x)$. Достаточным, но не необходимым условием является детальный баланс , который требует, чтобы каждый переход $x\to x'$ обратимо: для каждой пары состояний $x,x'$, вероятность находиться в состоянии $x$ и переход в состояние $x'$ должна быть равна вероятности нахождения в состоянии $x'$ и переход в состояние $x$, $\pi (x)P(x'\mid x)=\pi (x')P(x\mid x')$.
2. Уникальность стационарного распределения : стационарное распределение $\pi (x)$ должен быть уникальным. Это гарантируется эргодичностью марковского процесса, который требует, чтобы каждое состояние (1) было апериодическим — система не возвращается в одно и то же состояние через фиксированные промежутки времени; и (2) быть положительно рекуррентными — ожидаемое число шагов для возврата в то же состояние конечно.

Алгоритм Метрополиса – Гастингса включает в себя разработку марковского процесса (путем построения вероятностей перехода), который удовлетворяет двум вышеуказанным условиям, так что его стационарное распределение $\pi (x)$ выбран быть $P(x)$. Вывод алгоритма начинается с условия детального баланса :

$P(x'\mid x)P(x)=P(x\mid x')P(x'),$

который переписывается как

${\frac {P(x'\mid x)}{P(x\mid x')}}={\frac {P(x')}{P(x)}}.$

Подход состоит в том, чтобы разделить переход на два подэтапа; предложение и принятие-отказ. Распространение предложений $g(x'\mid x)$ это условная вероятность предложения состояния $x'$ данный $x$и акцептирующее распределение $A(x',x)$ это вероятность принять предложенное состояние $x'$. Вероятность перехода можно записать как произведение их:

$P(x'\mid x)=g(x'\mid x)A(x',x).$

Подставив это соотношение в предыдущее уравнение, получим

${\frac {A(x',x)}{A(x,x')}}={\frac {P(x')}{P(x)}}{\frac {g(x\mid x')}{g(x'\mid x)}}.$

Следующим шагом в выводе является выбор коэффициента приемки, который удовлетворяет приведенному выше условию. Одним из распространенных вариантов является выбор Метрополиса:

$A(x',x)=\min \left(1,{\frac {P(x')}{P(x)}}{\frac {g(x\mid x')}{g(x'\mid x)}}\right).$

Для этого коэффициента приемки Метрополиса $A$, или $A(x',x)=1$ или $A(x,x')=1$ и в любом случае условие выполняется.

Таким образом, алгоритм Метрополиса – Гастингса можно записать следующим образом:

1. Инициализировать
   1. Выберите начальное состояние $x_{0}$.
   2. Набор $t=0$.
2. Итерировать
   1. Создать случайное состояние-кандидат $x'$ в соответствии с $g(x'\mid x_{t})$.
   2. Рассчитать вероятность принятия $A(x',x_{t})=\min \left(1,{\frac {P(x')}{P(x_{t})}}{\frac {g(x_{t}\mid x')}{g(x'\mid x_{t})}}\right)$.
   3. Принять или отклонить :
      1. генерировать равномерное случайное число $u\in [0,1]$;
      2. если $u\leq A(x',x_{t})$, затем примите новое состояние и установите $x_{t+1}=x'$;
      3. если $u>A(x',x_{t})$, затем отклоните новое состояние и скопируйте старое состояние вперед $x_{t+1}=x_{t}$.
   4. Приращение : установлено $t=t+1$.

При выполнении заданных условий эмпирическое распределение сохраненных состояний $x_{0},\ldots ,x_{T}$ подойдет $P(x)$. Количество итераций ( $T$) требуется для эффективной оценки $P(x)$ зависит от ряда факторов, в том числе от взаимоотношений между $P(x)$ а также распределение предложений и желаемая точность оценки. Для распределения в дискретных пространствах состояний оно должно быть порядка времени автокорреляции марковского процесса.

Важно отметить, что в общей задаче неясно, какое распределение $g(x'\mid x)$ следует использовать или количество итераций, необходимое для правильной оценки; оба являются свободными параметрами метода, которые необходимо адаптировать к конкретной решаемой задаче.

Использование в численном интегрировании

Основная статья: Интеграция Монте-Карло

Алгоритм Метрополиса – Гастингса обычно используется для вычисления интеграла. В частности, рассмотрим пространство $\Omega \subset \mathbb {R}$ и распределение вероятностей $P(x)$ над $\Omega$, $x\in \Omega$. Метрополис – Гастингс может оценить интеграл вида

$P(E)=\int _{\Omega }A(x)P(x)\,dx,$

где $A(x)$ представляет собой (измеримую) функцию интереса.

Например, рассмотрим статистику $E(x)$ и его распределение вероятностей $P(E)$, что является маргинальным распределением . Предположим, что цель состоит в том, чтобы оценить $P(E)$ для $E$ на хвосте $P(E)$. Формально, $P(E)$ можно записать как

$P(E)=\int _{\Omega }P(E\mid x)P(x)\,dx=\int _{\Omega }\delta {\big (}E-E(x){\big )}P(x)\,dx=E{\big (}P(E\mid X){\big )}$

и, таким образом, оценивая $P(E)$ может быть достигнуто путем оценки ожидаемого значения индикаторной функции $A_{E}(x)\equiv \mathbf {1} _{E}(x)$, который равен 1, когда $E(x)\in [E,E+\Delta E]$ и ноль в противном случае.Потому что $E$ находится на хвосте $P(E)$, вероятность нарисовать состояние $x$ с $E(x)$ на хвосте $P(E)$ пропорционально $P(E)$, что мало по определению. Здесь можно использовать алгоритм Метрополиса – Гастингса для более вероятной выборки (редких) состояний и, таким образом, увеличения количества выборок, используемых для оценки. $P(E)$ на хвостах. Это можно сделать, например, с помощью выборочного распределения. $\pi (x)$ в пользу этих государств (например, $\pi (x)\propto e^{aE}$ с $a>0$).

Пошаговые инструкции

Три цепи Маркова, работающие на трехмерной функции Розенброка с использованием алгоритма Метрополиса – Гастингса. Цепочки сходятся и перемешиваются в области, где функция высока. Освещено примерное положение максимума. Красные точки — это те точки, которые остаются после процесса приработки. Более ранние были отброшены.

Предположим, что самое последнее выбранное значение равно $x_{t}$. Следуя алгоритму Метрополиса – Гастингса, мы затем рисуем новое состояние предложения. $x'$ с плотностью вероятности $g(x'\mid x_{t})$ и рассчитать стоимость

$a=a_{1}a_{2},$

где

$a_{1}={\frac {P(x')}{P(x_{t})}}$

это отношение вероятностей (например, байесовское апостериорное) между предложенной выборкой $x'$ и предыдущий образец $x_{t}$, и

$a_{2}={\frac {g(x_{t}\mid x')}{g(x'\mid x_{t})}}$

— соотношение плотности предложения в двух направлениях (от $x_{t}$ к $x'$ и наоборот).Это значение равно 1, если плотность предложений симметрична.Тогда новое состояние $x_{t+1}$ выбирается по следующим правилам.

Если $a\geq 1{:}$

$x_{t+1}=x',$

еще:

$x_{t+1}={\begin{cases}x'&{\text{with probability }}a,\\x_{t}&{\text{with probability }}1-a.\end{cases}}$

Цепь Маркова запускается с произвольного начального значения $x_{0}$, и алгоритм выполняется множество итераций, пока это начальное состояние не будет «забыто». Эти образцы, которые выбрасываются, известны как выгорание . Оставшийся набор принятых значений $x$ представляют выборку из распределения $P(x)$.

Алгоритм работает лучше всего, если плотность предложений соответствует форме целевого распределения. $P(x)$, из которого прямая выборка затруднена, т.е. $g(x'\mid x_{t})\approx P(x')$.Если плотность предложения по Гауссу $g$ используется параметр дисперсии $\sigma ^{2}$ необходимо настроить в период приработки.Обычно это делается путем расчета коэффициента принятия , который представляет собой долю предложенных образцов, принятую в окне последнего $N$ образцы.Желаемая степень принятия зависит от целевого распределения, однако теоретически было показано, что идеальная скорость принятия для одномерного распределения Гаусса составляет около 50%, уменьшаясь примерно до 23% для $N$-мерное гауссово целевое распределение. Эти рекомендации могут хорошо работать при выборке из достаточно регулярных байесовских апостериорных данных, поскольку они часто следуют многомерному нормальному распределению, которое можно установить с помощью теоремы Бернштейна-фон Мизеса .

Если $\sigma ^{2}$ слишком мала, цепочка будет перемешиваться медленно (т. е. скорость приема будет высокой, но последовательные образцы будут медленно перемещаться по пространству, и цепочка будет лишь медленно сходиться к $P(x)$). С другой стороны,если $\sigma ^{2}$ слишком велико, процент принятия будет очень низким, поскольку предложения, скорее всего, попадут в регионы с гораздо более низкой плотностью вероятности, поэтому $a_{1}$ будет очень мал, и снова цепочка будет сходиться очень медленно. Обычно распределение предложений настраивают так, чтобы алгоритмы принимали порядка 30% всех выборок – в соответствии с теоретическими оценками, упомянутыми в предыдущем абзаце.

Байесовский вывод

Основная статья: Байесовский вывод

MCMC можно использовать для получения выборок из апостериорного распределения статистической модели.Вероятность принятия определяется выражением: $P_{acc}(\theta _{i}\to \theta ^{*})=\min \left(1,{\frac {{\mathcal {L}}(y|\theta ^{*})P(\theta ^{*})}{{\mathcal {L}}(y|\theta _{i})P(\theta _{i})}}{\frac {Q(\theta _{i}|\theta ^{*})}{Q(\theta ^{*}|\theta _{i})}}\right),$где ${\mathcal {L}}$ это вероятность , $P(\theta )$ априорная плотность вероятности и $Q$ (условная) вероятность предложения.

См. также

• Генетические алгоритмы
• Методы частиц среднего поля
• Легкий транспорт Метрополис
• Многократная попытка Метрополис
• Параллельный отпуск
• Последовательное Монте-Карло
• Имитация отжига

Дальнейшее чтение

• Бернд А. Берг . Моделирование Монте-Карло марковской цепью и их статистический анализ . Сингапур, World Scientific , 2004.
• Чиб, Сиддхартха; Гринберг, Эдвард (1995). «Понимание алгоритма Метрополиса – Гастингса» . Американский статистик , 49 (4), 327–335.
• Дэвид ДЛ Мин и До Ле Мин. «Понимание алгоритма Гастингса». Коммуникации в статистике – моделирование и вычисления, 44:2 332–349, 2015 г.
• Болстад, Уильям М. (2010) Понимание вычислительной байесовской статистики , John Wiley & Sons ISBN   0-470-04609-0