С нуля: байесовский вывод, цепь Маркова, Монте-Карло и Метрополис Гастингс, на питоне

Привет и добро пожаловать в первую статью из серии «С нуля», в которой я объясняю и реализую / создаю что-либо с нуля.

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

Как отмечает Хилари Мейсон:

«Если вы хотите использовать новый алгоритм, который вы не очень понимаете, лучший подход - реализовать его самостоятельно, чтобы узнать, как он работает, а затем использовать библиотеку, чтобы извлечь выгоду из надежного кода».

Вот почему я предлагаю объяснить и реализовать с нуля: байесовский вывод (несколько кратко), цепь Маркова, Монте-Карло и Метрополис Гастингс, на Python.

Блокнот и версию в формате pdf можно найти в моем репозитории по адресу: joseph94m

Предпосылки: основные теории вероятностей, исчисление и Python.

1. Введение

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

Тогда я не особо об этом думал. Я подумал: «О, это просто еще одна техника отбора проб», и решил, что буду читать об этом, когда мне это практически понадобится. Эта потребность никогда не возникала, или, возможно, она возникла, и я ошибочно использовал что-то еще для решения своей проблемы.

1.1- Так почему интерес сейчас?

Недавно я увидел несколько дискуссий о MCMC и некоторых его реализациях, в частности об алгоритме Метрополиса-Гастингса и библиотеке PyMC3. Больше всего мое внимание привлекла статья Монте-Карло цепи Маркова в Python: полная реализация в реальном мире. В этой статье Уильям Кёрсен объясняет, как ему удалось изучить этот подход, применив его к реальной проблеме: оценить параметры логистической функции, которая представляет его режимы сна.

Г-н Кёрсен использует реализацию PyMC3 алгоритма Метрополиса-Гастингса для вычисления пространства распределения α и β, получая, таким образом, наиболее вероятную логистическую модель.

1.2- Так почему я говорю обо всем этом?

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

Я подумал, что если я запачкаю руки, я, наконец, смогу это понять. Я буду использовать только numpy для реализации алгоритма и matplotlib для рисования красивых вещей. В качестве альтернативы scipy можно использовать для вычисления функций плотности (о которых я расскажу позже), но я также покажу, как их реализовать с помощью numpy.

1.3- Постановка статьи

В части 1 я представлю байесовский вывод, MCMC-MH и их математические компоненты. Во второй части я объясню алгоритм MH с использованием фиктивных данных. Наконец, в части 3 будет представлено реальное приложение для MCMC-MH.

2- Часть 1: Байесовский вывод, цепь Маркова, Монте-Карло и Метрополис-Гастингс

2.1- Взгляд с высоты птичьего полета на философию вероятностей

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

2.1.1 - Частотное мышление против байесовского мышления

Есть две основные интерпретации вероятностей: байесовская и частотная.

С точки зрения Frequentist, вероятности представляют собой долгосрочные частоты, с которыми происходят события. Специалист по частотам может сказать, что вероятность выпадения решки при подбрасывании монеты в долгосрочной перспективе равна 0,5. Каждый новый эксперимент можно рассматривать как один из бесконечной последовательности возможных повторений одного и того же эксперимента. Основная идея заключается в том, что частотный подход к вероятности не верен. Вероятность возникновения события x из n испытаний приблизительно равна:

а истинная вероятность достигается при n - ›∞. Фраквенционисты никогда не скажут: «Я на 45% (0,45) уверен, что сегодня есть лазанья на обед», поскольку в долгосрочной перспективе этого не происходит. Обычно частотный подход называется объективным подходом, поскольку в нем нет выражения убеждений и / или предшествующих событий.

С другой стороны, в байесовском мышлении вероятности трактуются как выражение веры. Поэтому для байесовца вполне разумно сказать: «Я уверен на 50% (0,5), что сегодня есть лазанья на обед». Комбинируя предыдущие убеждения и текущие события (свидетельства), можно вычислить апостериорное, то есть вероятность того, что сегодня есть лазанья. Идея, лежащая в основе байесовского мышления, состоит в том, чтобы постоянно обновлять убеждения по мере поступления новых доказательств. Поскольку этот подход имеет дело с убеждениями, его обычно называют субъективным взглядом на вероятность.

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

В философии принятия решений байесовский вывод тесно связан с байесовским взглядом на вероятность в том смысле, что он манипулирует априорными, свидетельствами и правдоподобием для вычисления апостериорной. Учитывая какое-то событие B, какова вероятность того, что событие A произойдет? На это отвечает знаменитая формула Байеса:

С участием:

В нашем случае нас больше всего интересует конкретная формулировка формулы Байеса:

То есть мы хотели бы найти наиболее вероятное распределение θ, параметров модели, объясняющей данные, D.

Вычисление некоторых из этих вероятностей может быть утомительным, особенно доказательства P (D). Кроме того, могут возникнуть другие проблемы, например проблемы с обеспечением сопряжения, на которые я не буду углубляться в этой статье. К счастью, некоторые методы, а именно MCMC, позволяют нам брать выборку из апостериорного и рисовать распределения по нашим параметрам, не беспокоясь ни о вычислении доказательств, ни о сопряженности.

2.1.3- Цепь Маркова Монте-Карло

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

Итог: Его можно использовать для вычисления распределения по параметрам с учетом набора наблюдений и предварительного убеждения.

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

MCMC - это класс методов. Metropolis-Hastings - это конкретная реализация MCMC. Он хорошо работает в пространствах большой размерности в отличие от выборки Гиббса и выборки отбраковки.

Этот метод требует простого распределения, называемого распределением предложения (которое я люблю называть моделью перехода) Q (θ ′ / θ), чтобы помочь извлечь выборки из трудноразрешимого апостериорного распределения P (Θ = θ / D).

Metropolis-Hastings использует Q, чтобы случайным образом бродить по распределительному пространству, принимая или отклоняя переходы на новые позиции в зависимости от вероятности выборки. Это «беспамятное» случайное блуждание - часть «цепи Маркова» MCMC.

«Вероятность» каждой новой выборки определяется функцией f. Вот почему f должен быть пропорционален апостериорной части, из которой мы хотим получить выборку. f обычно выбирается как функция плотности вероятности, которая выражает эту пропорциональность.

Чтобы получить новую позицию параметра, просто возьмите текущую, θ, и предложите новую, θ ’, то есть случайную выборку, взятую из Q (θ ′ / θ). Часто это симметричное распределение. Например, нормальное распределение со средним значением θ и некоторым стандартным отклонением σ: Q (θ ′ / θ) = N (θ, σ)

Чтобы решить, следует ли принять или отклонить θ ′, необходимо вычислить следующее соотношение для каждого нового предложенного θ ’:

Используя формулу Байеса, это можно легко переформулировать как:

Что также эквивалентно:

Правило принятия можно сформулировать следующим образом:

Это означает, что если θ ’более вероятно, чем текущее θ, то мы всегда принимаем θ’. Если он менее вероятен, чем текущий θ, то мы можем принять его или отклонить случайным образом с уменьшающейся вероятностью, тем менее вероятно, что это так.

Итак, вкратце, алгоритм Метрополиса-Гастингса делает следующее:

3- Часть 2: Пример фиктивных данных

3.1- Шаг 1: Генерация данных

Мы генерируем 30000 выборок из нормального распределения со средним μ = 10 и стандартным отклонением σ = 3, но мы можем наблюдать только 1000 случайных выборок из них.

3.2- Шаг 2: Чего мы хотим?

Мы хотели бы найти распределение для σ {наблюдаемое}, используя 1000 наблюдаемых образцов. Те из вас, кто является заядлым математиком, скажут, что существует формула для вычисления σ:

Почему мы хотим пробовать и еще много чего? Ну, это всего лишь пример фиктивных данных, настоящая проблема в части 3, где трудно вычислить параметры напрямую. Кроме того, здесь мы не пытаемся найти значение для σ, а, скорее, мы пытаемся вычислить распределение вероятностей для σ.

3.3- Шаг 3: Определите PDF и модель перехода

Из рисунка 1 видно, что данные распределяются нормально. Среднее значение можно легко вычислить, взяв среднее значение из 1000 выборок. Делая это, мы получаем, что μ {наблюдаемый} = 9,8 (хотя на стороне примечания, мы могли бы также предположить, что μ неизвестно и выбрать для него так же, как мы делаем для σ. Тем не менее, я хочу сделать этот начальный пример просто.)

3.3.1- Для модели перехода / распространения предложения

Я не имею в виду конкретное распределение, поэтому выберу простое: Нормальное распределение!

Обратите внимание, что σ ′ не имеет отношения к σ {new} и σ {current}. Он просто указывает стандартное отклонение пространства параметров. Это может быть любое желаемое значение. Это влияет на время сходимости алгоритма и корреляцию между выборками, о которых я расскажу позже.

3.3.2- Для PDF

Поскольку f должен быть пропорционален апостериорному , мы выбираем f как следующую функцию плотности вероятности (PDF) для каждой точки данных di в набор данных D :

3.4- Шаг 4: Определите, когда мы принимаем или отклоняем σ {новый}

3.4.1- Формула принятия

Если это отношение не больше 1, то мы сравниваем его с равномерно сгенерированным случайным числом в замкнутом множестве [0,1]. Если отношение больше случайного числа, мы принимаем σ {new}, в противном случае мы отклоняем его. Это гарантирует, что даже если образец менее вероятен, чем текущий, мы все равно можем захотеть его попробовать. (Аналогично моделированному отжигу)

3.4.2- Вероятность

Общая вероятность для набора наблюдения D равна:

3.4.3- Приор P (μ, σ)

У нас нет никаких предпочтений в отношении значений, которые могут принимать σ {new} и σ {current}. Единственное, что стоит отметить, это то, что они должны быть положительными. Почему? Интуитивно стандартное отклонение измеряет дисперсию. Дисперсия - это расстояние, и расстояния не могут быть отрицательными.

Математически:

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

3.4.4- Форма окончательного принятия

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

Наше условие приемлемости из уравнения (1) выглядит следующим образом:

Эта форма может быть уменьшена еще больше, если извлечь квадратный корень и умножить из журнала, но неважно, что сейчас математики уже достаточно!

3.4.5. Реализация этой разглагольствования

Реализация довольно простая, правда ?!

3.6- Шаг 6: Запустите алгоритм с начальными параметрами и соберите принятые и отклоненные образцы.

Алгоритм принял 8317 выборок (которые могут отличаться при каждом новом запуске). Последние 10 отсчетов содержат следующие значения σ:

[2.87920187, 3.10388928, 2.94469786, 3.04094103, 2.95522153, 3.09328088, 3.07361275, 3.08588388, 3.12881964, 3.03651136]

Давайте посмотрим, как алгоритм работал до этих значений:

Итак, начиная с начального значения σ, равного 0,1, алгоритм довольно быстро сходился к ожидаемому значению 3. Тем не менее, это всего лишь выборка в одномерном пространстве… так что это не очень удивительно.

Тем не менее, мы будем рассматривать начальные 25% значений σ как «выгорание», поэтому мы их опускаем.

3.6.2. Визуализируем след σ и гистограмму следа

Наиболее вероятное значение σ составляет около 3,1. Это немного больше, чем исходное значение 3,0. Разница связана с тем, что мы наблюдаем только 3,33% от первоначальной популяции (1000 из 30 000).

3.6.3. Прогнозы. Как наша модель сможет спрогнозировать исходную численность населения в 30 000 человек?

Во-первых, мы усредняем последние 75% принятых выборок σ и генерируем 30000 случайных индивидов из нормального распределения с μ = 9,8 и σ = 3,05 (среднее значение последних 75% принятых выборок), что на самом деле лучше, чем наиболее вероятное значение 3,1.

И вуаля:

А теперь перейдем к делу!

4- Часть 3: Пример из реальной жизни

Солнечное пятно - это область на поверхности Солнца (сфера), температура которой ниже, чем температура окружающей среды. Эти пониженные температуры вызваны концентрацией потока магнитного поля, который подавляет конвекцию за счет эффекта, подобного вихретоковому торможению. Солнечные пятна обычно появляются парами с противоположной магнитной полярностью. Их количество варьируется в зависимости от примерно 11-летнего солнечного цикла.

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

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

4.2- Кажется, мы могли бы смоделировать это явление с помощью гамма-распределения, с новым циклом, сбрасываемым каждые 12 лет.

Гамма-распределение Γ - это двухпараметрическое семейство непрерывных распределений вероятностей. Параметры - форма a и масштаб b. Случайная величина X, имеющая гамма-распределение, обозначается как X ~ Γ (a, b), а в нашем случае X - это количество солнечных пятен. Два параметра a и b - это неизвестные, для которых мы хотели бы вычислить распределения.

Например, в первом цикле количество солнечных пятен начинается с самого высокого значения около 300 в конце 1749 года и падает до самого низкого значения через 6 лет, в течение 1755 года. Затем число снова возрастает до максимального значения в 1761 и 1762 годах до этого. снова падение в 1766 году и так далее. . .

Давайте в этом убедимся, построив гистограмму количества солнечных пятен:

4.3- Действительно, кажется, что частота отсчетов следует гамма-распределению.

Гамма-распределение для PDF имеет такое f, что:

где Γ - гамма-функция: Γ (a) = (a-1)! (не путать с гамма-распределением!)

Следуя той же процедуре, что и в примере с фиктивными данными, мы можем записать логарифмическую вероятность из этого PDF-файла (см. Код ниже). В качестве альтернативы можно использовать функцию scipy.stats.gamma (a, b) .pdf (x) для его вычисления. Однако имейте в виду, что реализация scipy на несколько порядков медленнее, чем реализованная мной.

Поскольку a и b должны быть строго положительными, мы обеспечиваем это в предыдущем случае:

Запустите код и соберите образцы:

Начиная с a = 4 и b = 10, алгоритм принял 8561 пару выборок, последнее значение для a - 0,98848982, а последнее значение для b - 84,99360422, что довольно далеко от начальных значений.

Как и в примере с фиктивными данными, давайте посмотрим, как алгоритм работал до этих значений:

Как видно из рисунков 10, 11 и 12, алгоритм быстро сходится к зоне [a = 1, b = 85].

Совет: когда алгоритм начинает сильно отклонять выборки, это означает, что мы достигли зоны насыщения вероятности. Обычно это можно интерпретировать как достижение оптимального пространства параметров, из которого мы можем производить выборку, т.е. у алгоритма очень мало причин для принятия новых значений. Это отмечено на рисунках 11 и 12, где алгоритм больше не принимает значения за пределами небольшого диапазона.

4.3.1 - Мы считаем начальные 50% значений a и b «выгоранием», поэтому мы их опускаем. Давайте визуализируем кривые и b и гистограмму кривых.

4.4- Время прогноза

Сначала мы усредняем последние 50% принятых выборок a и b и генерируем случайных людей из Γ-распределения. a {средний} = 0,9866200759935773 и b {средний} = 83,70749712447888.

И прогнозы:

4- Оценка

4.1- Оценка распространения предложения

Как указать параметры для распределения Q? Должны ли мы отойти от текущей выборки θ или оставаться относительно близко? На эти вопросы можно ответить, измерив автокорреляцию между принятыми выборками. Мы не хотим, чтобы удаленные выборки были слишком коррелированы, поскольку мы пытаемся реализовать цепь Маркова, т.е. выборка должна зависеть только от своей предыдущей выборки, а график автокорреляции должен показывать быстрое экспоненциальное уменьшение между корреляцией выборки i и i-1, i-2,… ik

Автокорреляция определяется путем вычисления следующей функции для каждого лага k:

Графики ниже показывают автокорреляцию для a, b для k от 1 до 100. Запаздывание k = 0 означает, что мы измеряем корреляцию выборки с самим собой, поэтому мы ожидаем, что она будет равна 1. чем выше k, тем ниже должна быть корреляция.

В нашем случае нам повезло с достаточно низкой корреляцией. В общем, мы можем захотеть настроить параметры распределения предложения Q автоматически. Распространенный метод - постоянно корректировать параметры предложения так, чтобы отклонялось более 50% предложений. В качестве альтернативы можно использовать расширенную версию MCMC, называемую гамильтонианом Монте-Карло, которая уменьшает корреляцию между последовательными выборочными состояниями и быстрее достигает стационарного распределения.

6. Заключение

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

Я надеюсь, что всем, кто прочитал эту статью, она показалась интересной и познавательной. Если будут положительные отзывы, появятся и другие статьи из этой серии «С нуля», в которых я объясню и реализую вещи с нуля (очевидно!), Поэтому, если вам это нравится, пожалуйста, предложите, о чем вы хотите, чтобы я рассказал дальше!

Приветствуются любые вопросы, и я постараюсь на них ответить! Отзывы более чем приветствуются, так как это моя первая статья, и я хочу улучшить ее.

Использованная литература:

Питер Дрисколл, «Сравнение методов наименьших квадратов и байесовской аппроксимации с наборами данных о лучевых скоростях»

Карсон Чоу, «MCMC и подгонка моделей к данным»

Джон Х. Уильямсон, «Основы данных - вероятности»

Саймон Роджерс, «Первый курс машинного обучения»

Благодарности:

Соавтор медиума, Вера Чернова, за статью: «Как встроить Jupyter Notebook в посты Medium за три шага: 1–2–3!»

Друзьям и коллегам по анализу данных Марку Аскью и Константиносу Иоаннидису за отзывы об этой статье перед публикацией.