Цепь Маркова Монте-Карло с нуля: глубокое погружение в Метрополис-Хастинг и выборка Гиббса

Что такое цепь Маркова Монте-Карло

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

Что такое цепь Маркова

Цепь Маркова можно рассматривать как цепной метод выборки, в котором распределение каждой выборки зависит только от предыдущей выбранной выборки. Ниже приведен простой пример цепи Маркова для завтрака:

«У Тома есть только три варианта завтрака за всю его жизнь. А) Хлеб Б) Зерновые и В) Яйца. То, что он ест сегодня, зависит только от того, что он ел вчера. Например, если он съел (C) яйца вчера, у него будет 20% шанс съесть (C) яйцаеще раз и 20% шанс съесть (A) хлеб и 70% сдачу (B) хлопья. Точно так же у него также есть фиксированный образец выбора завтрака, если он ел вчера (А) или (Б). Это формирует цепочечную структуру выборки из трех штатов».

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

Обратите внимание, что мы не рассмотрели определенные условия, которые будут объяснены позже. Вопрос здесь в том, если предположить, что это правда, почему это полезно для выборки из неразрешимого распределения? Это полезно, потому что, если мы можем построить цепь Маркова, стационарное распределение которой является близким приближением к сложному целевому распределению, мы можем затем просто произвести выборку из построенной цепи Маркова как процесс выборки из целевого распределения. Это означает, что мы можем преобразовать задачу выборки из сложного распределения в построение цепи Маркова, стационарное состояние которой является целевым распределением. Это все еще нетривиально, но сейчас у нас может быть более многообещающее направление для работы. Давайте поговорим о том, как построить цепь Маркова.

Из чего состоит цепь Маркова

Прежде чем мы углубимся в то, как построить цепь Маркова для моделирования сложного распределения, мы сначала вернемся к оставшейся проблеме о том, каким определенным условиям должна удовлетворять цепь Маркова. В дискретном случае цепь Маркова определяется так называемым оператором переходаT, который записывает вероятность перехода между различными состояниями. Если мы хотим, чтобы структура вероятности, подобная цепочке, давала стационарное состояние после множества выборок, подобная цепочке структура, то есть оператор перехода, должна удовлетворять следующим условиям:

1. Нередуцируемость. Это означает, что в цепи Маркова каждое состояние может быть достигнуто из любого другого состояния.

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

3. Подробный баланс Это говорит о том, что вероятность перехода из A в B должна быть равна вероятности перехода из B в A. Почему это так? Потому что только при этом условии Цепообразная структура может сходиться к стационарному распределению. Смотрите следующее доказательство взорвали.

Подробное подтверждение баланса

Вероятность перехода из A в B состоит из двух частей: a. A выбирается в целевом распределении P(A) b. вероятность того, что B переходит в данный A T(B|A). Итак, следующее:

Поэтому,

Как построить цепь Маркова для целевого распределения

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

Предположим, что целевое распределение равно π(x); если мы можем построить переходное ядро P, удовлетворяющее условию π(x) P= π(x), мы закончили. Это связано с тем, что в этом случае мы можем просто выполнить выборку, перескакивая между различными распределениями с переходным ядром P и аппроксимируя целевое распределение. Давайте поговорим о переходе Kernel P, раз уж он так важен. Формальное определение переходного ядра:

Ядро марковского перехода на пространстве с мерой (X, χ) — это функция P: X × χтакая что

1. P(x,·) — вероятностная мера для всех x∈ X

2. P(·, A) — измеримая величина для всех A∈ χ

Интуитивно первое условие говорит о том, что P(x, A) можно использовать для описания условного распределения вероятностей следующего состояния A при заданном текущем состоянии x. Это легко понять, поскольку это то же самое, что делает матрица перехода. Второй из них можно игнорировать, поскольку это технический термин, обозначающий, что ядро интегрируемо.

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

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

Обратите внимание на объяснение обложек, что это такое и почему оно работает. Мы не будем углубляться в мотивацию/мыслительный процесс, лежащий в основе дизайна, так как это требует довольно много усилий.

Дизайн ядра алгоритма MH:

где,

p(x, y) – это распределение x при заданном y. (например, p(x, y) = N∼ (x|y,1)). Таким образом, интеграл от p(x, y) по — это сумма вероятности того, что он переместится в следующую точку. Тогда r(x) — это вероятность того, что цепочка останется в точке x. δx(dy) — функция Дирака в dy, которая есть не что иное, как δx(dy) = 1, если x∈ dyи = 0 в противном случае.

Теперь мы определили функцию ядра P, и каждый термин ясно объяснен. Как мы уже говорили ранее, мы не собираемся объяснять, почему он разработан таким образом, мы просто докажем, почему это работает. Мы знаем, что пока ядро перехода удовлетворяет условию π P= π, мы можем производить выборку из этого ядра, чтобы получить стационарное распределение π .

Как проверить, будет ли ядро перехода расслаиваться π P= π. Вернитесь к уравнению подробного баланса: пока P(x, y) удовлетворяет подробному балансу, оно будет иметь стационарное распределение.

Быстрый способ доказать это, обнаруженный Тирни в 1994 году, состоит в том, что, пока p(x, y) удовлетворяет подробному балансу, πявляется инвариантной плотностью P. Доказательство показано ниже:

Итак, опять же, пока p(x, y) удовлетворяет условию подробного баланса, как показано в строках 2–3 в уравнении выше, это работает.

Теперь вы можете видеть, что проблема действительно упростилась до того, как найти распределение p(x, y), то есть распределение вероятностей y при заданном x, которое имеет подробное условие баланса. с целевым распределением π. p(x, y)π(y) = p(y, x)π(x)

Хорошо, честно, но сделать это по-прежнему нетривиально, верно? Предположим, что мы знаем π (в большинстве случаев, вероятно, нет), будем ли мы решать это аналитически? Нет, предлагаемое решение очень простое.

Мы предположили, что q(x, y) будет очень простым базовым распределением, из которого легко выбрать образец, как в случае нормального распределения. Скажем, q(x, y) = N ~(y| x,1). Ясно, что это наивное предложение не удовлетворяет условию детального баланса. Чтобы сбалансировать его, мы вводим вероятность принятия α(x, y) в левой части, чтобы повторно сбалансировать уравнение. Теперь мы имеем:

Скажем, если LHS в уравнении 4 велико, α(x, y) будет меньше 1. Если RHS велико, α(x, y) будет больше 1. Как только мы это получим, детали функции сбалансированы. π — инвариантное распределение.

Теперь мы можем взять q и α вместе как p в исходном переходном ядре, подставив их в предложение P, мы будем иметь:

Это обратимо и имеет инвариантную плотность π, как мы построили таким образом. Давайте просто остановимся на мгновение, чтобы посмотреть, что мы сделали.

Мы пытались найти распределение p(x, y), которое детализировало бы баланс с π. Это нетривиально, если мы решили это аналитически. Вместо этого алгоритм MH умело разработал способ принудительного соблюдения баланса деталей, предложив простое распределение q(x, y) и отношение α (x, y) и их умножение составляют p(x,y). Сq(x, y) очень легко иметь дело, поскольку это то, что мы предложили. α(x, y) можно рассчитать по коэффициенту распределения. Нам просто нужно немного потрудиться, чтобы найти α на каждом этапе (поскольку α изменится, как только вы перейдете к следующей точке), баланс деталей всегда сохраняется.

Каково на практике значение α(x, y)? Название α(x, y) портит ответ: мы используем это как приемлемую скорость в процессе выборки. Давайте проработаем это:

Предположим, что здесь нет α(x, y) и есть q(x, y ), который является полностью обратимым (подробные условия баланса). Мы просто берем каждую точку этой выборки из q(x,y), потому что из приведенного выше уравнения она точно дает π после определенного количества выборок.

Теперь q(x, y) не является идеальным распределением, поэтому построенная на нем цепь Маркова не идеальна. Мы хотим контролировать частоту плохой выборки с помощью дополнительного параметра α(x, y). Если α меньше 1, это означает, что LHS π(x)q(x, y) больше RHS π(y)q(y, x), что означает вероятность того, что текущая точка будет выбрана, на самом деле больше, чем ваша следующая точка, поэтому не переходите к этой точке легко. Давайте перейдем туда только с вероятностью α. С другой стороны, если α больше 1, что означает, что вероятность того, что следующая точка будет выбрана выше, чем ваша текущая точка, наверняка вы захотите туда попасть. Поэтому на практике имеем следующую норму приемки.

Кроме того, поскольку уровень принятия является отношением, скаляр целевого распределения не имеет значения. Это большое преимущество алгоритма MH в том, что вам нужно знать целевое распределение только до постоянного скаляра, это отлично подходит для оценки апостериорного распределения.

Итак, в целом, ниже показано, как работать с алгоритмом MH.

Или посмотрите это YouTube Video, чтобы узнать, как это работает на практике.

Выборка Гиббса

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

Предположим, что мы знаем условное распределение π(·|x1) и π(·|x2) , и мы хотим выбрать из π(x1, x2).

Мы собираемся выбрать из π(·|xi), построив цепь Маркова, т.е. переходное ядро P(x2, у2 x1). Мы можем сделать это так же, как алгоритмы MH, т.е. предложить распределение и использовать скорость принятия для удовлетворения баланса деталей. Поскольку предполагается, что условное распределение относительно легко выбрать, мы собираемся использовать точное условное распределение в качестве предполагаемого распределения, т. е. π(x , y) = q(x, y) Что из этого следует? Это означает, что коэффициент приемлемости α наверняка будет равен 1. Ниже приводится доказательство.

Предположим, что мы делаем это для бинарного распределения дисперсии, зависящего от x2, y2, мы используем MH для выборки x1(текущий) y1 (следующий образец).

Следующий шаг, мы повторим это, при условии x1, y1, для выборки x2 y2 , что также даст коэффициент принятия 1.

Подытожим идею выборки Гиббса:

1. Построение переходного ядра P(x2, y2|x1) и P (x1, y1|x2), которые условные распределения π( .|x1) и π(.|x2) в качестве целевых распределений.

2. Произведение двух ядер условного перехода может дать стационарное распределение π(x1, x2), которое является совместным распределением что нас интересует. (Доказательство ниже)

Произведение двух ядер условного перехода имеет целевое распределение в виде инвариантной плотности. Итак, что нам нужно сделать сейчас, так это взять образец из продукта переходных ядер? Но что значит сэмплировать из продукта ядер? По свойству переходного ядра:

если P и Q являются переходными ядрами, произведение Переходное ядро PQ соответствует первому выполнению перехода относительно P, а затем Q

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

Вот как выполняется выборка Гиббса

Или загляните в это YouTube Video.

Заключение

Алгоритмы алгоритмов Метрополиса-Хастинга и выборки Гиббса довольно просты в реализации, и вы можете найти множество материалов, объясняющих, как работают эти алгоритмы. Этот блог посвящен математической части, которая поддерживает эти два алгоритма. Я надеюсь, что вы смогли следовать и нашли это полезным.