что такое метод монте карло

Интуитивное использование методов Монте-Карло с цепями Маркова

Легко ли это? Я попробовал

Алексей Кузьмин, директор разработки и работы с данными «ДомКлик», преподаватель направления Data Science в Нетологии, перевел статью Rahul Agarwal о том, как работают методы Монте-Карло с цепями Маркова для решения проблем с большим пространством состояний.

Каждый, кто связан с Data Science, хоть раз слышал о методах Монте-Карло с цепями Маркова (MCMC). Иногда тема затрагивается во время изучения байесовской статистики, иногда — при работе с такими инструментами, как Prophet.

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

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

Итак, начнем. MCMC состоит из двух терминов: Монте-Карло и Марковские цепи. Поговорим о каждом из них.

Монте-Карло

В самых простых выражениях методы Монте-Карло можно определить как простую симуляцию.

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

Марковские цепи

Поскольку термин MCMC состоит из двух частей, нужно еще понять, что такое цепи Маркова. Но прежде, чем перейти к Марковским цепям, немного поговорим о Марковских свойствах.

Предположим, есть система из M-возможных состояний, и вы переходите из одного состояния в другое. Пусть пока вас ничего не сбивает с толку. Конкретный пример подобной системы — погода, которая меняется от жаркой к холодной и умеренной. Другой пример — фондовый рынок, который прыгает от «медвежьего» к «бычьему» и состоянию стагнации.

Марковское свойство говорит о том, что для данного процесса, который находится в состоянии X_n в конкретный момент времени, вероятность X_{n + 1}= k (где k — любое из M-состояний, в которое процесс может перейти) зависит только от того, каково это состояние в данный момент. А не о том, как оно достигло нынешнего состояния.
Говоря математическим языком, можем записать это в виде следующей формулы:

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

Процесс, обладающий Марковским свойством, называется Марковским процессом. Из-за чего же важна Марковская цепь? Из-за своего стационарного распределения.

Что такое стационарное распределение

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

У вас есть матрица переходных вероятностей, которая определяет вероятность перехода из состояния X_i в X_j.

Матрица переходных вероятностей, Q

В приведенной матрице переходных вероятностей Q вероятность того, что следующим состоянием будет «бык», учитывая текущее состояние «бык» = 0,9; вероятность того, что следующее состояние будет «медвежьим», если текущее состояние «бык» = 0,075. И так далее.

Что же, давайте начнем с какого-то определенного состояния. Наше состояние будет задаваться вектором [бык, медведь, стагнация]. Если мы начнем с «медвежьего» состояния, вектор будет таким: [0,1,0]. Мы можем вычислить распределение вероятности для следующего состояния, умножив вектор текущего состояния на матрицу переходных вероятностей.

Заметьте, что вероятности дают в сумме 1.

Следующее распределение состояний можно найти по формуле:

И так далее. В конце концов, вы достигнете стационарного состояния, в котором состояние стабилизируется:

Для описанной выше матрицы переходных вероятностей Q стационарное распределение s такое:

Вы можете получить стационарное распределение при помощи следующего кода:

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

Теперь мы можем ответить на вопрос, почему же стационарное распределение так важно.

Стационарное распределение важно, потому что с его помощью можно определять вероятность системы быть в определенном состоянии в случайное время.

Для нашего примера можно сказать, что в 62,5% случаев рынок будет в «бычьем» состоянии, 31,25% — в «медвежьем» и 6,25% — в стагнации.

Интуитивно вы можете рассматривать это как случайное блуждание по цепи.

Случайное блуждание

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

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

Теперь перейдем непосредственно к рассмотрению методов MCMC.

Что такое методы Монте-Карло с цепями Маркова (MCMC)

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

Думайте…

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

На рисунке показана теорема Байеса

Например, нужно сделать сэмпл из апостериорного распределения. Но легко ли вычислить апостериорную компоненту вместе с нормирующей постоянной (evidence)? В большинстве случаев можно найти их в виде произведения правдоподобия и априорной вероятности. Но рассчитать нормирующую постоянную (p (D)) не получится. Почему? Разберемся подробнее.

Предположим, H принимает только 3 значения:

p(D) = p(H=H1).p(D|H=H1) + p(H=H2).p(D|H=H2) + p(H=H3).p(D|H=H3)

В таком случае p(D) подсчитать легко. А что если значение H непрерывно? Получилось бы ли рассчитать это так же легко, особенно если H принимала бесконечные значения? Для этого пришлось бы решить сложный интеграл.

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

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

Разберемся на примере. Допустим, нужна выборка из бета-распределения. Его плотность:

где C — нормирующая постоянная. На самом деле это некоторая функция от α и β, но я хочу показать, что для выборки из бета-распределения она не нужна, поэтому мы примем ее за константу.

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

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

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

Один из таких MCMC-алгоритмов — алгоритм Метрополиса-Гастингса.

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

Интуиция:

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

Так, например, мы хотели бы провести вдвое больше времени на вершине холма высотой 100 метров, чем на соседнем 50-метровом холме. Хорошо, что мы можем сделать это, даже если не знаем абсолютных высот точек на поверхности: все, что нужно знать, это относительные высоты. Например, если вершина холма A в два раза выше вершины холма B, то мы бы хотели провести в A вдвое больше времени, чем на B.

Существуют более сложные схемы предложения новых мест и правил их принятия, но основная идея заключается в следующем:

Алгоритм:

Теперь давайте определим и опишем поставленную задачу более формально. Пусть s = (s1, s2,…, sM) — желаемое стационарное распределение. Мы хотим создать цепь Маркова с таким стационарным распределением. Начнем с произвольной цепи Маркова с M-состояниями с матрицей перехода P, так, что pij представляет вероятность перехода из состояния i в j.

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

Проделывая это для выборки бета-распределения, единственный раз, когда приходится использовать плотность вероятности, — это поиск вероятности принятия решения. Для этого делим sj на si (то есть нормирующая постоянная C сокращается).

Выборка из бета-распределения

Теперь перейдем к проблеме выборки из бета-распределения.

Бета-распределение — непрерывное распределение на [0,1] и может иметь бесконечные значения на [0,1]. Предположим, что произвольная цепь Маркова P с бесконечными состояниями на [0,1] имеет переходную матрицу P такую, что pij = pji = все элементы в матрице.

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

Пора перейти от теории к практике. Напишем на Python наш бета-образец.

Сверим результаты с фактическим бета-распределением.

Как видно, значения очень похожи на бета-распределение. Таким образом, сеть MCMC достигла стационарного состояния

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

Выводы

Это был большой пост. Поздравляю, если дочитали его до конца.

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

Вы также можете использовать MCMC для решения проблем с большим пространством состояний. Например, в задаче о рюкзаке или для расшифровки. Я постараюсь предоставить вам более интересные примеры в следующем посте. Следите за обновлениями.

Источник

Методы Монте-Карло для марковских цепей (MCMC). Введение

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

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

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

Алгоритмы MCMC – это класс методов, предназначенных именно для этого. Они восходят к эпохальной статье Метрополиса и др.; Метрополис разработал первый MCMC-алгоритм, названный в его честь и предназначенный для вычисления равновесного состояния двумерной системы жестких сфер. Фактически, исследователи искали универсальный метод, который позволил бы вычислять ожидаемые значения, встречающиеся в статистической физике.
В этой статье будут рассмотрены основы выборки по MCMC.

Теперь, когда мы понимаем, зачем нам делать выборку, перейдем к сердцу MCMC: марковским цепям. Что такое марковская цепь? Не вдаваясь в технические детали, можно сказать, что марковская цепь – это случайная последовательность состояний в некотором пространстве состояний, где вероятность выбора определенного состояния зависит только от текущего состояния цепи, но не от ее прежней истории: эта цепь лишена памяти. В определенных условиях марковская цепь имеет уникальное стационарное распределение состояний, к которому сходится, преодолев определенное количество состояний. После такого числа состояний состояния в марковской цепи получают инвариантное распределение.

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

Марковская цепь будет обходить дискретное пространство состояний, образуемое тремя состояниями погоды:

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

В строках указаны состояния, в которых в настоящее время может находиться цепь, а в столбцах – состояния, в которые цепи могут перейти. Если взять «временной» шаг марковской цепи за один час, то, при условии, что сейчас солнечно, существует 60% вероятность, что солнечная погода сохранится в течение ближайшего часа. Также есть 30% вероятность, что в следующий час будет облачная погода, и 10%-я вероятность, что сразу же после солнечной погоды пойдет дождь. Это также означает, что значения в каждом ряду в сумме дают единицу.
Давайте немного погоняем нашу марковскую цепь:

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

РОДОНАЧАЛЬНИК ВСЕХ MCMC: АЛГОРИТМ МЕТРОПОЛИСА-ГАСТИНГСА

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

Если бы мы смогли спроектировать переходное ядро таким образом, чтобы следующее состояние уже было выведено из , то этим можно было бы и ограничиться, так как наша марковская цепь… непосредственно делала бы выборку из . К сожалению, чтобы добиться этого, нам нужна возможность делать выборку из , чего мы делать не можем – иначе вы бы этого не читали, верно?

Обходной путь – разбить переходное ядро на две части: шаг проб (proposal step) и шаг приема/отбрасывания (acceptance/rejection step). На шаге проб фигурирует вспомогательное распределение , из которого выбираются возможные следующие состояния цепочки. Мы не только можем делать выборку из этого распределения, но и в силах произвольно выбирать само распределение. Однако, при проектировании следует стремиться прийти к такой конфигурации, в которой образцы, взятые из этого распределения, минимально коррелировали бы с актуальным состоянием и одновременно имели хорошие шансы пройти этап приема. Вышеуказанный этап приема/отбрасывания – вторая часть переходного ядра; на данном этапе корректируются ошибки, содержавшиеся в пробных состояниях, выбранных из . Здесь вычисляется вероятность успешного приема и принимается проба с такой вероятностью в качестве следующего состояния в цепи. Получение следующего состояния из тогда выполняется следующим образом: сначала пробное состояние берется из . Затем оно принимается в качестве следующего состояния с вероятностью или отбрасывается с вероятностью , и в последнем случае актуальное состояние копируется и используется в качестве следующего.
Следовательно, имеем

Достаточным условием для того, чтобы марковская цепь имела в качестве стационарного распределения является следующее: переходное ядро должно подчиняться детальному равновесию или, как пишут в физической литературе, микроскопической обратимости:

.

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

.

А вот здесь начинается магия: известно нам только до константы, но это не имеет значения, поскольку данная неизвестная константа обнуляет выражение для ! Именно это свойство paccpacc обеспечивает работу алгоритмов, основанных на алгоритме Метрополиса-Гастингса, на ненормированных распределениях. Часто используются симметричные вспомогательные распределения с , и в таком случае алгоритм Метрополиса-Гастингса редуцируется до оригинального (менее общего) алгоритма Метрополиса, разработанного в 1953 году. В оригинальном алгоритме

.

В таком случае полное переходное ядро Метрополиса-Гастингса можно записать как

.

РЕАЛИЗУЕМ АЛГОРИТМ МЕТРОПОЛИСА-ГАСТИНГСА НА PYTHON

Отлично, теперь, когда мы разобрались в принципе работы алгоритма Метрополиса-Гастингса, давайте перейдем к его реализации. Сначала установим логарифмическую вероятность того распределения, из которого собираемся делать выборку – без нормировочных констант; предполагается, что мы их не знаем. Далее работаем со стандартным нормальным распределением:

Далее выбираем симметричное вспомогательное распределение. В целом, производительность алгоритма Метрополиса-Гастингса можно повысить, если включить во вспомогательное распределение уже известную вам информацию о том распределении, из которого вы хотите сделать выборку. Упрощенный подход выглядит так: берем текущее состояние x и выбираем пробу из , то есть, задаем некоторую величину шага Δ и переходим влево или вправо от нашего текущего состояния не более чем на :

Наконец, вычисляем вероятность того, что предложение будет принято:

Теперь собираем все это вместе в по-настоящему краткую реализацию этапа построения выборки для алгоритма Метрополиса-Гастингса:

ТЕСТИРУЕМ НАШ АЛГОРИТМ МЕТРОПОЛИСА-ГАСТИНГСА И ИССЛЕДУЕМ ЕГО ПОВЕДЕНИЕ

Наверное, теперь вам не терпится увидеть все это в действии. Этим и займемся, примем несколько информированных решений об аргументах stepsize и n_total :

Все хорошо. Итак, это сработало? Мы добились принятия проб примерно в 71% случаев, и у нас есть цепочка состояний. Следует отброс ить несколько первых состояний, в которых цепь еще не сошлась к своему стационарному распределению. Давайте проверим, на самом ли деле у выбранных нами состояний – нормальное распределение:

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

Мееееедленно идем к цели…

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

Настоятельно рекомендую вам самим поэкспериментировать с кодом, приведенным здесь – так вы освоитесь с поведением алгоритма в различных обстоятельствах и глубже его поймете. Попробуйте несимметричное вспомогательное распределение! Что будет, если не настроить критерий приема как следует? Что произойдет, если попытаться сделать выборку из бимодального распределения? Можете ли придумать способ автоматической настройки размера шага? Какие здесь есть подводные камни? Ответьте на эти вопросы самостоятельно!

Источник