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

MCMC изначально использовались для решения многократных интегралов численными методами, например, в байесовской статистике, вычислительной физике^[1], вычислительной биологии^[2] и компьютерной лингвистике^[3][4].

Недавние достижения в области MCMC позволили проводить вычисления в больших иерархических моделях, требующих интегрирования по сотням и тысячам переменных^[5].

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

Методы Монте-Карло с марковскими цепями создают выборки на основе выбранной непрерывной случайной величины с известной функцией плотности распределения. Эти выборки можно использовать для оценки интеграла по этой величине с помощью математического ожидания или дисперсии.

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

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

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

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

• Алгоритм Метрополиса-Гастингса: этот метод генерирует марковскую цепь, используя заданную плотность и фильтрацию новых шагов. Фактически это общая схема, частными случаями которой являются: самый первый и простой MCMC (алгоритм Метрополиса), а также альтернативы, перечисленные ниже.
  • Семплирование по Гиббсу: для этого метода требуется, чтобы были известны все условные распределения целевого распределения. Если вывод из полностью условных распределений не является непосредственным, то используются другие сэмплеры по Гиббсу (см. напр. ^[6][7][8]). Семплирование по Гиббсу популярно, потому что оно не требует предварительной «настройки».
  • Методы Монте-Карло по Ланжевену (Langevin Monte Carlo) и прочие методы, опирающиеся на градиент (и, возможно, вторую производную) логарифма целевой плотности. Цель этих методов – предложить наиболее правдоподобные шаги в направлении более высокой плотности вероятности^[9].
  • Псевдо-маргинальный метод Метрополиса-Гастингса: этот метод заменяет плотность целевого распределения его несмещённой оценкой. Метод применим, когда целевая плотность не задана аналитически (напр. в моделях со скрытыми переменными).
• Выборка по слоям: в основе этого метода лежит принцип, согласно которому выборку некоторого распределения можно построить, равномерно семплируя область под графиком функции плотности этого распределения. Этот метод заменяет однородное семплирование в вертикальном направлении однородным семплированием горизонтального «слоя», определённого текущей вертикальной позицией.
• Алгоритм многократных попыток Метрополиса (Multiple-try Metropolis): это вариация алгоритма Метрополиса-Гастингса, которая позволяет делать в каждой точке многократные попытки. Делая более крупные шаги на каждой итерации, алгоритм помогает избавиться от «проклятия размерностей»^[10][11].

• Метод обратимого прыжка (Reversible-jump) : ещё один вариант алгоритма Метрополиса-Гастингса, допускающий изменения размерности пространства^[12]. Такие методы марковской цепи долгое время использовались в прикладной статистической физике, где в некоторых задачах имело место распределение, называемое «большим каноническим ансамблем» (напр. при переменном количеством молекул в замкнутом сосуде). Метод обратимого прыжка применим при использовании MCMC или семплирования по Гиббсу на непараметрических байесовских моделях, где количество смешивающихся компонент (кластеров) автоматически прогнозируется с помощью данных (напр. процесс Дирихле или «процесс китайского ресторана»).
• Гамильтонов (гибридный) метод Монте-Карло: этот метод пытается избежать случайного блуждания введением дополнительного «вектора момента» и применением гамильтоновой динамики, где функция потенциальной энергии — это целевая плотность. Выборки момента отбрасываются после семплирования. Окончательный результат алгоритма состоит в том, что движение по пространству выборки делается крупными шагами. Таким образом, они менее коррелированы и быстрее сходятся к целевому распределению.

Интерактивные методологии МСМС – это класс методов частиц в среднем поле для получения выборок псевдослучайных чисел из последовательности распределений вероятности с возрастающим уровнем сложности семплирования^[13].

Эти вероятностные модели включают в себя:

• модели состояния пространства путей с возрастающим временным горизонтом
• апостериорные (по отношению к частичным наблюдениям) распределения
• возрастающий уровень ограничений для условных распределений
• убывающие температурные графики, связанные с распределениями Больцмана-Гиббса
• и многое другое

В общем случае любые сэмплеры MCMC можно превратить в интерактивные. Их, в свою очередь, можно использовать как способ параллельного запуска последовательности обычных сэмплеров. Например, интерактивные алгоритмы имитации отжига основаны на независимых перемещениях по алгоритму Метрополиса-Гастингса, которые последовательно взаимодействуют с механизмом выборочного ресемплирования. В противоположность классическим методам MCMC, здесь параметр точности интерактивных сэмплеров зависит только от их количества. Методы взаимодействующих частиц принадлежат классу моделей частиц Фейнмана-Каца^[14][15], также называемому «последовательным методом Монте-Карло» или «методами фильтра частиц» в теории байесовских выводов и в области обработки сигналов^[16]. Интерактивные методы МСМС можно также понимать как циклы в алгоритме генетических частиц с мутациями в виде классических МСМС мутаций.

Использование последовательностей с низким расхождением вместо случайных чисел для простого независимого семплирования Монте-Карло обладает явными преимуществами^[19]. Такая замена используется в методе квази-Монте-Карло (QMC)^[20]. Согласно неравенству Коксма-Хлавка, погрешность интеграции в таком методе гораздо меньше, чем при семплировании независимых одинаково распределённых случайных величин (IID). Это позволяет на порядок уменьшить и погрешность оценки и время сходимости.

Метод Array-RQMC вводит моделирование марковских цепей на основе QMC путём моделирования ${\displaystyle n}$ цепей одновременно. На каждом шаге эмпирическое распределение этих ${\displaystyle n}$ состояний даёт более точное приближение к функции распределения, чем при использовании MCMC^[21]. В эмпирических экспериментах сходимость дисперсии среднего функции состояния иногда имеет скорость порядка ${\displaystyle O(n^{-2})}$ или даже быстрее, вместо ${\displaystyle O(n^{-1})}$ у метода Монте-Карло^[22].

Обычно построить марковскую цепь с необходимыми свойствами нетрудно. Сложнее определить, сколько шагов потребуется для сходимости к стационарному распределению с приемлемой ошибкой^[23]. Хорошая цепь обладает свойством перемешивания: стационарное распределение достигается быстро для любой начальной позиции. Классический эмпирический метод достижения сходимости состоит в том, чтобы запустить несколько независимо смоделированных марковских цепей и проверить, что дисперсии вне цепи и внутри неё примерно равны^[23][24].

Обычно семплирование MCMC может только приближать целевое распределение, поскольку всегда присутствует остаточный эффект начальной позиции. Более сложные алгоритмы, основанные на MCMC, такие как coupling from the past, могут получать точные выборки, что сказывается на количестве вычислений и потраченного времени.

Многие методы случайного блуждания Монте-Карло движутся небольшими шагами, начиная от равновесного распределения без склонности прокладывать путь в одном направлении. Эти методы легко применять и анализировать, но для исследования всего пространства при помощи блуждания требуется длительное время (блуждание часто возвращается на уже пройденные участки).

Дальнейшие рассмотрения сходимости содержатся в Центральной предельной теореме для марковских цепей см.^[25] для обсуждения теории, относящейся к сходимости и стационарности алгоритма Метрополиса-Гастингса.

Программное обеспечение для работы с МСМС семплированием:

• пакеты, которые используют диалекты языка модели BUGS
  • WinBUGS / OpenBUGS / MultiBUGS
  • JAGS
  • NIMBLE
• greta, язык байесовского статистического моделирования
• MCSim
• PyMC3
• pymcmcstat
• R (язык программирования) с пакетами adaptMCMC, atmcmc, BRugs, mcmc, MCMCpack, ramcmc, rjags, rstan, и т. д.
• Stan
• TensorFlow Probability (библиотека вероятностного программирования, встроенная в TensorFlow)
• MCL (кластерный алгоритм для графов и HipMCL (альтернативная версия))^[26]
• emcee (реализация аффинного инвариантного семплера ансамбля MCMC на Python)
• Keanu Java-библиотека многоцелевого вероятностного программирования
• Zeus реализация ансамблевого семплирования по слоям на Python
• Turing.jl пакет для многоцелевого вероятностного программирования в Julia
• Mamba.jl платформа для МСМС метода в Julia