Пару месяцев назад, перед началом курса байесовских вычислений, который я сейчас изучаю, я хотел получить более глубокое представление о моделировании цепей Маркова Монте-Карло (MCMC). В то время я запускал свой блокнот Jupyter, импортировал PyMC3 и строил модель, заканчивая все это вызовом оптимизатора Metropolis, не совсем понимая, что происходит.
Внутренняя красота и элегантность метода кроются за несколькими десятками понятий и определений, что превращает все исследование в тяжелое испытание и мешает среднему студенту вроде меня оценить то, что на самом деле происходит под капотом. В этой статье я хотел бы представить более простое введение и практический пример использования MCMC с помощью алгоритма Метрополиса-Гастингса (MC). В частности, я сосредоточусь на выводе параметров для простой модели линейной регрессии. Однако слово «простой» не отражает всей глубины лежащей в его основе идеи. Замена нормального распределения любым другим условным распределением (Пуассона/Гамма/Отрицательное Биномиальное или другими) позволит реализовать почти идентичную реализацию GLM через MCMC, изменив максимум 4 или 5 строк кода. Кроме того, поскольку выборки для вывода берутся из апостериорного распределения параметров, мы вполне естественно можем использовать их для построения доверительных интервалов для параметров (в данном случае правильнее говорить «достоверные интервалы»).
В этой статье будет кратко описано, почему используются методы MCMC, она предложит реализацию алгоритма MH для модели линейной регрессии и завершится визуализацией того, что на самом деле происходит, когда алгоритм ищет набор параметров, которые сгенерировали данные.
На написание этой статьи меня вдохновила лекция на эту тему профессора Ника Херда.
1. Почему МСМС?
Представьте, что вы кошка — ну, это не обязательно, но потерпите — и вы окружены травой. Что-то вдалеке привлекает ваше внимание: это мороженое. На вершине горы. Точнее, не на вершине горы, а на самой вершине группы гор, которая окружена лесом, и горы сливаются друг с другом и скрыты какими-то очень непрозрачными облаками, которые мешают видеть сквозь них. Как добраться до мороженого?
Вы ищете золотую монету ацтеков, которую нашли в джунглях. И затем вы решаете, что каждый раз, когда вы хотите сделать шаг, вы будете дважды подбрасывать монету и следовать приведенным ниже правилам:
- если вы переворачиваете голову-голову, вы делаете шаг вперед.
- если вы перевернете хвост-хвост, вы сделаете шаг назад
- если вы перевернете голову-хвост, вы сделаете шаг вправо
- если вы перевернете хвост, вы сделаете шаг влево
Это разумная стратегия, чтобы добраться до мороженого? Я бы сказал, что это не так. Проблема здесь в том, что стратегия не информирована. То есть у нас нет оценки того, насколько вероятно следующий шаг приблизит нас к обновленной цели. Методы MCMC, и в частности MH, используют функцию правдоподобия, чтобы определить, следует ли принять или отклонить ход, и вводят стохастичность, чтобы избежать принятия всех шагов, которые кажутся приближающими нас к цели. .
В нашем следующем примере, который включает в себя реализацию регрессионной модели, мороженое представляет значение параметров, которые мы пытаемся вывести. Подробная информация представлена ниже.
2. Настройка проблемы
Пусть Y и X будут ответом и предиктором для модели соответственно. Линейная модель утверждает, что условное распределение ответа с учетом такого предиктора является нормальным. То есть:
для подходящих параметров a (наклон), b (смещение) и σ (интенсивность шума).
Наша задача — вывести a, b и σ.
"Ну, нам нужно заранее знать, как они распространяются", вы могли бы сказать. Ну да и нет. Вам нужно установить некоторые «основные правила», которым должна следовать модель. В нашем случае a, b могут быть любыми, положительными или отрицательными, но σ должен быть строго положительным (я никогда не слышал о нормальное распределение с отрицательным стандартным отклонением, не так ли?). Кроме этого, никаких других правил не требуется. Асимптотическая нормальность апостериорного распределения гарантирует, что пока у вас достаточновыборок из апостериорного распределения, они будут нормально распределены несмотря ни на что (если цепь Маркова достигнет своего стационарного распределения).
В этом случае давайте сгенерируем синтетические данные для регрессии, используя параметры a=3, b=20 и σ=5. Это можно сделать в Python с помощью следующего кода:
Синтетические данные можно увидеть на изображении ниже:
В следующем разделе будут определены функции для алгоритма Метрополиса Гастингса и цикла для заданного количества итераций.
3. Алгоритм
Пусть θ=[a,b,σ] — вектор параметров на предыдущем шаге алгоритма, а θ' — предложение для новый набор параметров, MH сравнивает фактор Байеса (произведение правдоподобия и априорности) между двумя конкурирующими гипотезами для параметров (θ' и θ) и масштабирует этот фактор с помощью обратной величины распределения условного предложения. Затем этот коэффициент сравнивается со значением равномерно распределенной случайной величины. Это добавляет модели стохастичность и позволяет исследовать маловероятные векторы параметров и отбрасывать вероятные (редко).
Теперь это звучит немного запутанно, не так ли? Давайте начнем с построения идей шаг за шагом.
3.1 Распространение предложения
Во-первых, мы определяем распределение предложения g(θ’|θ): это распределение для фиксированных параметров на предыдущем временном шаге. В нашем примере a и b могут быть как положительными, так и отрицательными, поэтому естественным выбором будет их выборка из многомерного нормального распределения с центром на предыдущем шаге итерации. Что касается σ, мы можем, например, выбрать его из гамма-распределения. Это действительно зависит от пользователя, определяющего все эти дистрибутивы. Лучший подход (который мы не будем здесь рассматривать) будет включать выборку σ из обратного гамма-распределения.
В результате мы можем определить распределение предложения следующим образом:
И:
Параметр k используется для управления «разбросом» распределений вокруг их средних значений. 𝜔 — это дополнительная настройка гаммы. Чем больше k𝜔, тем больше разброс для гаммы (мы используем формулировку скорости формы гамма-распределения, будьте осторожны, так как scipy заставляет нас использовать вместо этого формулировку шкалы формы). Функция ниже:
Вы можете заметить, что функция proposal содержит два аргумента: prec_theta представляет собой вектор параметров на предыдущем шаге, а search_width представляет регион, в котором мы будем искать предложения. Поиск хорошего набора параметров приводит к компромиссу между исследованием(поиском новых параметров, установленных в неисследованных областях) и эксплуатацией(уточнением поиск ближе к областям, где были найдены хорошие наборы параметров). Поэтому с search_width следует обращаться крайне осторожно. Большое значение search_widthпредотвратит сходимость MCMC к стационарному распределению. Небольшое значение может помешать алгоритму найти оптимум (оптимум) за разумное время (потребуется большее количество выборок и более длительный период приработки). ожидал).
3.2 Функция правдоподобия
Функция правдоподобия выводится из выбора линейной модели. Условное распределение отклика при заданных параметрах является нормальным. Другими словами, мы будем вычислять вероятность для нормального распределения, где среднее значение является произведением предикторов и коэффициентов a и b, а шум σ. В частности, в этом случае мы будем использовать логарифмическое правдоподобие вместо исходного правдоподобия, чтобы улучшить числовую стабильность. В результате произведения просто превратятся в суммы:
3.3 Приор
Нам не нужно тратить на это много слов. Мы действительно можем выбрать все, что нам нравится с точки зрения распределения для априора. Все, что нам нужно сделать, это убедиться, что зона, в которой, вероятно, будут обнаружены предполагаемые параметры, имеет ненулевой априор и что параметр шума σ неотрицательный. Опять же, априор выражается в формате log-pdf.
3.4 Соотношение предложения
Отношение предложений — это вероятность наблюдения старого вектора параметров при наличии нового, деленная на вероятность наблюдения нового вектора параметров при наличии старого. То есть из раздела 3.1 g(θ|θ’)/g(θ’|θ). Опять же, мы будем использовать логарифмическую плотность вероятности, чтобы иметь единообразие шкал по вероятностям и добиться лучшей числовой стабильности.
3.5 Последний штрих
Мы почти на месте. Все, что нам нужно сделать, это собрать все вместе. В псевдокоде:
И в Питоне:
Потрясающий! Мы сделали это! Имейте в виду, что приведенный выше код займет некоторое время для запуска. Вы можете попробовать изменить N, чтобы получить меньше выборок, или отредактировать ширину, чтобы ускорить исследование (даже если вы заполняете больше эксплойтов и можете расходиться в своем поиске).
4. Результаты
Результаты представлены на картинке ниже.
Как мы видим, история образцов достаточно стабильна. Среднее значение и стандартное отклонение для истории выборки составляют (после отбрасывания первой половины истории, как выгорание):
Ну, есть одна проблема, о которой я только упомяну здесь. Выборки сильно коррелированы, поэтому при оценке достоверных интервалов может потребоваться дополнительная осторожность. Одно из решений заключается в прореживании истории путем сохранения лишь небольшого процента параметров (например, сохранение только 1 из 10 принятых смежных предложений и отбрасывание остальных).
5. Как эти результаты соотносятся с традиционной линейной регрессией?
Мы можем выполнить точно такое же вычисление, используя statsmodels.api.OLS
Со следующим результатом:
Среднее значение коэффициентов a и b очень близко к MCMC, которое можно найти в таблице 1.То же самое можно сказать для стандартной ошибки, что еще больше подкрепляет наши представления о достоверности нашей реализации MCMC.
Имейте в виду, что это не лучшее решение, а просто решение. Лучший выбор априорных предложений существует и рекомендуется.
6. Как выглядят итерации?
Я создал GIF итераций MH. К сожалению, в 3D визуализация довольно грязная, поэтому я решил сосредоточиться только на наклоне a и смещении b.
В визуализации:
- красная точка представляет текущее предложение
- серая область вокруг красной точки представляет собой распределение предложения (нормальное) при 3 стандартных отклонениях от среднего. Предложение имеет диагональную ковариационную матрицу, поэтому мы получаем круг вместо эллипса.
- Электрические синие линии представляют отклоненные ходы.
- Красные линии представляют принятые ходы
- Синяя точка представляет собой среднее значение параметров, полученных из
statsmodels.api.OLS.
Medium не позволяет загружать большие GIF, поэтому мне пришлось сжимать объем информации и показывать поведение MH раз в 10 шагов.
Надеюсь, вам понравилась статья! Ваше здоровье!
