Источник: Nuances of Programming
Слева : моделированное необработанное совместное распределение коэффициентов
Справа : моделированное совместное распределение коэффициентов без отбраковки
В предыдущей статье я дал краткое введение в байесовскую статистику и рассказал, как байесовский анализ объединяет априорные убеждения и экспериментальные данные, чтобы получить апостериорное распределение интересующего параметра. Задача, которую я использовал в статье, была выбрана потому, что в ней апостериорное распределение легко находится аналитически. Но, увы, чаще всего этот подход не работает, и необходимо численное решение. Алгоритм Монте-Карло по схеме цепей Маркова (Markov Chain Monte Carlo, далее MCMC) является одним из методов моделирования апостериорного распределения некоторого параметра.
Если у вас есть базовые знания байесовской статистики, смело читайте дальше. Если нет, взгляните на эту статью — в ней есть всё, что вам нужно знать, прежде чем приступить к изучению MCMC.
Алгоритм Метрополиса-Гастингса
Хорошим введением в MCMC является алгоритм Метрополиса-Гастингса . В нём 5 шагов. Давайте определим некоторые параметры и функции перед тем, как углубиться в сам алгоритм.
θ — случайный параметр, для которого моделируется распределение
Y — случайная переменная, зависящая от θ
y — реализация случайной переменной Y (выборка)
f_ θ (θ) — априорное распределение (обновляющееся на каждом шаге значением предыдущего периода)
f_Y|θ(y| θ) — правдоподобие
f_ θ|Y(θ| y) — апостериорное распределение
J_ θ(θ|θ_t-1) — скачкообразное распределение, используемое для формирования предсказанных новых значений θ
Теперь, когда мы определили параметры, переменные и функции, пройдём по шагам алгоритм Метрополиса-Гастингса.
Шаг 1 :
Выбираем начальное значение: θ_t -1 . При t = 1 , θ_t -1 — наше начальное априорное значение для θ .
Шаг 2 :
Формируем предсказание для периода t - θ*_t - путём выборки из скачкообразного распределения J(θ|θ_t -1 ). Метод Монте-Карло в MCMC проявляется в том, что случайное число моделируется из скачкообразного распределения. Цепи Маркова проявляются в том, что θ*_t зависит только от значения θ_t -1 . Процесс “не обладает памятью”, следовательно является цепью Маркова.
Шаг 3
Находим отношение апостериорного предсказания к апостериорному распределению априорного значения.
Затем перегруппируем для симметричных скачкообразных значений.
Шаг 4
Формируем случайное число a , принадлежащее к равномерному непрерывному распределению(0 ,1 ).
Если r > a , θ_t = θ*_t (мы принимаем предполагаемое значение). Вероятность min(r, 1 ).
Если r ≤ a , θ_t = θ_t -1 (мы отвергаем предсказанное значение и используем априорное значение). Вероятность 1 — min(r, 1 ).
Если r ≥ 1 (θ*_t более вероятно, чем θ_t -1 ), всегда будет принято. Иногда мы принимаем θ*_t , потому что не хотим попасть в ловушку в локальном минимуме.
Шаг 5
Повторяем, пока не достигнем максимального количества итераций.
Не так уж и страшно, когда вы знаете, как это работает, верно?
Обычно в качестве скачкообразного распределения используется нормальное распределение, потому что центрированное в 0 оно может создавать предсказания как меньшие, так и большие относительно предыдущего предсказания. Среднеквадратичное отклонение скачкообразного распределения является важным параметром настройки. Если априорное распределение находится далеко от области плотности апостериорного распределения на числовой оси, а среднеквадратичное отклонение скачкообразного распределения невелико, пройдёт немало времени, прежде чем у вас получится формировать предсказания из верной области числовой оси. Если среднеквадратичное отклонение скачкообразного распределения слишком велико по отношению к реальному распределению параметра, в моделированном апостериорном распределении будет много пробелов, потому что в этих областях не сформировано ни одного предсказания. К счастью, пакеты для MCMC содержат адаптивные среднеквадратичные отклонениях для скачкообразных распределений, поэтому на практике вы можете не беспокоиться по поводу этого параметра.
После формирования апостериорного распределения стандартной практикой является отбраковка области. Это связано с тем, что алгоритму может потребоваться слишком много итераций для поиска вероятных значений в пространстве параметров. Чтобы использовать 10% отбраковки, просто отбросьте первые 10% моделированных значений. После удаления отбракованной части получим моделированное апостериорное распределение, которое можно использовать для байесовского вывода — построения байесовских доверительных интервалов.
Применение модифицированного алгоритма Метрополиса-Гастингса к логистической регрессии
Теперь, когда мы поняли алгоритм Метрополиса-Гастингса, мы можем модифицировать его, чтобы он соответствовал коэффициентам логистической регрессионной модели. Расчётными параметрами будут 𝛽𝑗 для 𝑗=1,2,…𝑘 , где 𝑘 — число коэффициентов в модели. Основное отличие состоит в том, что вместо одного параметра будут рассчитываться несколько. Я справился с этим, формируя предсказание для каждого 𝛽𝑗 каждого цикла, при этом рандомизируя порядок каждого цикла 𝛽𝑗 .
Чтобы построить алгоритм MCMC для логистической регрессионной модели, нужно определить 4 функции. Это позволит нам вычислить соотношение апостериорных распределений для предсказанных 𝛽𝑗 на каждом шаге алгоритма MCMC.
Первая функция
inv_logit , отменяющая логит-преобразование. Логит-преобразование превращает вероятность в логит-функцию.
Выходные данные логистической регрессии представлены в формате логит-функции. Функция inv_logit нужна нам для переложения логит-функции в вероятности при расчёте логарифмического правдоподобия для вектора 𝛽 с учётом данных.
Вторая функция
normal_log_prior , вычисляющая логарифм априорного распределения вектора 𝛽 с учётом априорных убеждений о математических ожиданиях и среднеквадратичных отклонениях отдельных 𝛽𝑗 . Это натуральный логарифм плотности многомерного нормального распределения, где каждый элемент равен 𝑁∼(априорное математическое ожидание j, априорное среднеквадратичное отклонение j ) при 𝛽 . Натуральный логарифм используется для предотвращения исчезновения значащих разрядов.
Третья функция
log_likelihood , вычисляющая логарифм правдоподобия вектора 𝛽 с учётом данных. Поскольку мы моделируем задачу второго класса с помощью логистической регрессии, вероятность однократного наблюдения вычисляется по формуле Бернулли:
Правдоподобие вектора 𝛽 с учётом данных является произведением отдельных вероятностей Бернулли.
Следовательно, логарифмическое правдоподобие:
Натуральный логарифм используется для предотвращения исчезновения значащих разрядов.
Четвёртая функция
log_posterior рассчитывает апостериорное распределение вектора 𝛽 с учётом данных и наших априорных убеждений. апостериорное распределение ∝ априорное распределение∗правдоподобие, таким образом 𝑙𝑛(апостериорного распределения) ∝ 𝑙𝑛(априорного распределения)+𝑙𝑛(правдоподобие) .
Определив эти функции, я смог реализовать модифицированный алгоритм Метрополиса-Гастингса. Кроме того, я добавил функции для удаления отбракованной области, вычисления байесовских доверительных интервалов и прогнозирования. Чтобы сделать прогноз, пользователь должен указать, хочет ли он использовать среднее значение выборки, математическое ожидание или модальное значение моделированных апостериоров для каждого 𝛽𝑗 в качестве коэффициентов при вычислении 𝛽 .
После определения класса протестируем его на реальных данных.
Применение MCMC логистической регрессии к реальным данным
Чтобы протестировать модель логистической регрессии MCMC, я использовал этот набор данных , содержащий данные о 86 сладостях. Я выбрал простую модель, содержащую выделение отрезка на оси и “процент цены” (“Процентиль цены за единицу товара в сравнении с остальными товарами в наборе.”), для предсказания вероятности того, что сладость окажется шоколадкой.
Вот модель:
Как видите, моделированное совместное распределение 𝛽0 и 𝛽1 формирует эллиптическую область (как и совместные распределения в обобщённых линейных моделях). На графике совместного распределения с отбраковкой коэффициенты “идут” от начала координат, немного проходят по кругу и затем достигают эллиптической области. Эти коэффициенты, идущие от априорных значений в эллиптическую область, мы можем определить как “истинное” моделированное совместное распределение.
В 95% случаев байесовский доверительный интервал для 𝛽0 равен [-3,67; -1,32], а для 𝛽1 равен [2,50; 6,72]. Моделированное среднее значение выборки 𝛽0 равно -2,44, выборки 𝛽1 — 4,48. Более высокая цена положительно влияет на вероятность того, что сладость окажется шоколадом.
Сравнивая графики процентиля цены и вероятности оказаться шоколадом, видим, что модель имеет красивую S-образную форму, демонстрируемую в учебниках по логистической регрессии. Не похоже, что она проделывает отличную работу по классификации сладости как шоколада (слишком много ложных негативных значений в сравнении с истинными положительными), поэтому, вероятно, нужна более сложная модель.*
*Безусловно стоит использовать проверки и тестовые наборы данных для оценки производительности модели при построении реальных прогнозирующих моделей. Процесс обучения, проверки и тестирование выходят за рамки этой статьи. Здесь мы рассматриваем только применение MCMC.
Сравниваем с логистической регрессией, рассчитанной с помощью оценки максимального правдоподобия
Теперь, когда мы настроили совместную работу логистической регрессии и модифицированного алгоритма Метрополиса-Гастингса, сравним результаты со стандартной оценкой максимального правдоподобия (ОМП). Результаты очень похожи!
from statsmodels.discrete.discrete_model import Logit
# добавляем выделение отрезка на оси, так как статистический модуль
# этого не делает my_data['Intercept' ] = 1
# логистическая регрессия с использованием ОМП mle_mod = Logit(my_data[target], my_data[['Intercept' ] + vars_of_interest])
mle_mod_fit = mle_mod.fit(disp=False )
# печатаем print(mle_mod_fit.summary())
Если я всё сделал верно, вы поняли идею MCMC, в частности, алгоритм Метрополиса-Гастингса, а также применение к таким практическим задачам, как подбор логистической регрессионной модели.
Спасибо! Если вы интересуетесь MCMC, байесовским статистическим моделированием и машинным обучением, взгляните на pymc3 .
Репозиторий для этого проекта можно найти здесь .
Читайте также:
Читайте нас в Telegram , VK и Яндекс.Дзен
Перевод статьи John Clements : Intro to Markov Chain Monte Carlo