Что такое приближения Монте-Карло и как работает алгоритм Метрополиса?

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

Приближения Монте-Карло

Итак, что мне делать, если я не хочу аналитически вычислять какой-то надоедливый интеграл? Введите приближения Монте-Карло.

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

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

Вычисление интегралов с использованием выборок Монте-Карло - это все очень хорошо, но как вообще извлечь выборки из целевого распределения? Нарисовать гауссовские или однородные образцы легко, но что-то более сложное и np.random вас подведет.

Самый простой способ рисования выборок - использовать метод обратной функции CDF, но он основан на наличии доступа к функции обратной функции CDF, которая часто не имеет хорошей аналитической формы и имеет смысл только для одномерных случайных величин.

Алгоритм Метрополиса является одним из строительных блоков многих методов выборки методом Монте-Карло с цепью Маркова (MCMC). Это позволяет нам рисовать образцы, когда все, к чему у вас есть доступ, - это pdf-файл целевого дистрибутива.

В методах MCMC есть предостережение, что мы больше не берем независимые выборки, поэтому у нас не может быть таких же гарантий того, как дисперсия нашей оценки уменьшается с увеличением количества взятых нами выборок. Если выборки отбираются независимо, Центральная предельная теорема говорит нам, что дисперсия нашей оценки будет уменьшаться обратно пропорционально количеству выборок (N). Для MCMC мы можем уменьшить это значение, изменив количество отсчетов от N до N_eff. N_eff (почти) всегда меньше N и связано с тем, насколько коррелированы выборки в цепочке.

Выборка мегаполиса

Шаги алгоритма Метрополиса следующие:

1. Sample a starting point uniformly from the domain of the target distribution or from the prior distribution.
2. Calculate the pdf at that point.
3. Make a proposal for the new sample by taking a step from the current position in accordance with some state transition function.
4. Calculate the new pdf value.
5. Calculate the value of new pdf / old pdf.
6. If the ratio is greater than 1, accept the step.
7. If the ratio is less than 1:
    1. Sample a uniform random variable.
    2. If the ratio greater than the random sample, accept the step.
8. Else reject the proposal, add the current position as a sample and take a new step.

Обратите внимание, что процесс, описанный в 5–8, эквивалентен принятию выборки, основанной на вероятности Бернулли с вероятностью min (1, p (новая) / p (старая)), запомните это на будущее ...

Почему работает сэмплинг в Метрополисе?

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

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

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

Итак, теперь давайте покажем, как алгоритм Metropolis удовлетворяет этому условию ...

Чтобы найти p (i, j), удовлетворяющее подробному балансу, мы сначала предлагаем произвольную вероятность «скачка» q (i, j), а затем получаем p (i, j), принимая только «скачок» с вероятностью α ( i, j). Когда «скачок» отклоняется, состояние остается j = i. Эта идея «принятия» не уникальна для алгоритма Metropolis и существует в большинстве вариантов MCMC.

Это зависит от того, является ли α допустимым распределением вероятностей. Итак, допустимая форма α:

Если скачкообразная вероятность симметрична, то ее можно упростить до:

В противном случае его можно оставить в полном виде, который называется Metropolis-Hasting MCMC.

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

Вы также могли заметить, что поскольку альфа является функцией π (j) / π (i). В результате это означает, что целевое распределение не нужно нормализовать. Это особенно полезно при использовании Метрополиса для байесовской апостериорной оценки, когда трудно вычислить свидетельство.

Замечания по реализации

Распространенная версия алгоритма Метрополиса называется «Случайное блуждание по Метрополису», где предлагаемое состояние представляет собой текущее состояние плюс многомерный гауссовский алгоритм с нулевым средним и ковариационной матрицей σ²I. σ следует выбирать достаточно большим, чтобы отбраковалось достаточно много образцов. Это необходимо для хорошего изучения целевого распределения.

Второе, на что следует обратить внимание, - это концепция выгорания. Выборки, взятые до того, как цепь Маркова достигнет стационарного распределения, должны быть удалены, поскольку они не являются репрезентативными для целевого распределения до тех пор, пока цепь не сойдется. Довольно сложно определить, сколько образцов следует удалить, но в целом удаляется 10–20% образцов (или 10–100 эффективных образцов).

Давайте построим это в Numpy

Здесь для простоты мы реализуем случайное блуждание по Метрополису.

def metropolis(pi, dims, n_samples, burn_in=0.1, var=1):
    theta_ = np.random.randn(dims)*var
    samples = []
    while len(samples) < n_samples:
        theta = theta_ + np.random.randn(dims)*var
ratio = pi(theta)/pi(theta_)
        if np.random.rand(1) < ratio:
            sample = theta
            theta_ = theta
            samples.append(sample)
        else:
            sample = theta_
            samples.append(sample)
    samples = np.array(samples)
    return samples[int(samples*burn_in):,:]

Мы можем видеть, как это работает на сумме двух гауссианов (обратите внимание, что это ненормализованное распределение).

from scipy.stats import multivariate_normal
def make_pdf(mean1, mean2, cov1, cov2):
    pdf1 = multivariate_normal(mean1, cov1)
    pdf2 = multivariate_normal(mean2, cov2)
    def pdf(x):
        return pdf1.pdf(x) * pdf2.pdf(x)
    return pdf
pdf = make_pdf(
    [3, 3],
    [-1, -1],
    np.array([[1,0.1],[0.1,1]], dtype=float),
    np.array([[1,0.1],[0.1,1]], dtype=float),
)
samples = metropolis(pdf, 2, 1_000, 0.1)

На гифке выше показано, как алгоритм проходит распределение, иногда переключаясь между двумя разными режимами распределения. Обратите внимание, что это также выявило одну из слабых сторон алгоритма Метрополиса, он относительно плохо справляется с многомодельными распределениями. Это связано с тем, что для изучения нового режима шаг должен быть достаточно большим, чтобы переходить от одного режима к другому. Это требует либо большого размера шага, либо распределения с близкими друг к другу модами. В этом могут помочь такие модификации, как гамильтониан MCMC, но, как правило, это проблема большинства методов MCMC.

Дополнительное «чтение»

  • Видео на детальном балансе
  • Оригинальная статья об алгоритме Метрополиса-Гастингса
  • Бумага о гамильтоновой динамике (метод MCMC, используемый станом и другими вероятностными языками)