Байесовская статистика 101

Любите это или ненавидите, вы больше никогда не будете смотреть на статистику одинаково

Вступление

Байесовская статистика отличается от классической статистики (также известной как частотная статистика) в основном своей интерпретацией вероятности. Первый рассматривает это как «степень уверенности», тогда как второй рассматривает это как «относительную частоту, наблюдаемую во время многих испытаний».

Это различие может показаться слишком абстрактным, но оно имеет большое практическое влияние на методы, разработанные впоследствии.

Частотный подход является наиболее распространенным и, следовательно, повлиял на большее количество статистических методов. Однако многие современные методы основаны на байесовском подходе и могут дать отличные результаты, предоставляя вам альтернативы, когда данных для работы не так много.

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

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

Теорема Байеса

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

P (A): вероятность события A
P (A | B): вероятность события A при условии, что событие B произошло

Событие может быть чем угодно. Например, P (A | B) может означать «вероятность (P) наличия COVID (A) при условии, что ( |) ваш ПЦР-тест оказался положительным (B) ». Чтобы рассчитать эту вероятность с помощью приведенного выше уравнения, нам потребуются:

P (A): вероятность заражения COVID (независимо от результатов теста)
P (B): вероятность положительного результата теста ПЦР (независимо от того, есть ли у вас COVID или нет)
P (B | A): вероятность положительного результата теста ПЦР, если у вас COVID.

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

Байесовская статистика

С этого момента байесовская статистика превращает эту теорему в следующую мантру:

У вас есть предварительное представление о том, как устроен мир. Как только вы получите данные, вы соответствующим образом обновите это убеждение.

Это означает, что вы начинаете решение любой проблемы с определения предшествующего распределения (этот шаг весьма субъективен), а затем используете вероятность наблюдаемых данных для обновления этого предыдущего, создавая апостериорное распределение.

Не ясно? Давайте рассмотрим пример решения проблемы, и вы увидите, как это происходит на практике.

Проблема Монти Холла

Давайте заключим сделку - популярное игровое телешоу, которое стартовало в 60-х годах в США. Первоначальным ведущим которого звали Монти Холл. Знаменитая вероятностная головоломка, основанная на ней, впоследствии стала известной в следующем формате:

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

Вы выбираете Дверь 1. Затем ведущий открывает Дверь 3, обнаруживая, что за ней находится коза. Затем он дает вам выбор: придерживаться двери 1 или переключиться на дверь 2. Что вы делаете?

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

Давайте посмотрим на проблему с байесовской точки зрения.

Начнем с априорного распределения вероятностей 1/3 для каждой двери. Это означает, что каждая дверь априори имеет 1 шанс из 3 оказаться «правильной» дверью (есть 3 двери, и у нас нет оснований полагать, что одна из них более вероятна, чем другая). Априор - это P (A) в теореме Байеса.

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

Дверь 1: хост будет выбирать между дверями 2 и 3 случайным образом, поэтому вероятность того, что он откроет дверь 3 в этом случае, равна 1/2.
Дверь 2: хост должен будет открыть дверь 3 (поскольку машина находится за дверью 2, а вы уже выбрали дверь 1), поэтому вероятность того, что он откроет дверь 3 в этом случае, равна 1
Дверь 3: хост не мог открыть дверь 3, если бы за ним находился автомобиль, поэтому вероятность того, что он откроет дверь 3, в этом случае равна 0.

Это оставляет нам:

Вероятность эквивалентна P (B | A) в теореме Байеса. Теперь давайте посчитаем произведение «Приора» и «Вероятность»:

Вы заметите, что 3 значения (1/6, 1/3 и 0) не дают в сумме 1. Это потому, что нам не хватает последнего элемента теоремы Байеса: P (B). Это сумма этих трех значений. Разделив каждое на эту сумму, мы масштабируем их таким образом, чтобы они прибавляли к 1. В нашей задаче 3 значения складываются до 1/2 (1/6 + 1/3 = 1/2). Чтобы найти апостериор, нам нужно только разделить последний столбец на 1/2:

Как и ожидалось, вероятность того, что машина находится за дверью 3, равна 0, поскольку хост ее уже открыл. Мы также видим, что дверь 2 в 2 раза чаще скрывает машину, чем дверь 1.

Если вас особенно интересует эта проблема или результат кажется вам слишком противоречивым, есть целая статья об этом здесь.

Хорошо, я знаю, что это была игрушечная проблема, и что мы использовали очень простой дистрибутив в качестве предыдущего. Что происходит, когда мы сталкиваемся с реальными сложными проблемами с множеством переменных, где априорные значения не так легко определить, а вероятность трудно вычислить? К счастью, в этом нам поможет библиотека Python: PyMC3.

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

Давайте посмотрим, как использовать оба варианта для классических задач, таких как регрессия и классификация.

Регресс

Линейная регрессия «с нуля»

Если вы знакомы с линейной регрессией, вы заметите несколько отличий от традиционных моделей. Главный из них заключается в том, что в байесовской регрессии мы рассматриваем параметры α, β1, β2 и σ2 не как фиксированные значения, а как переменные, следующие за распределениями.

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

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

[IN]:
import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc3 as pm
# Setting a random seed for reproducibility
np.random.seed(23)
# True parameter values
alpha = -0.8
beta = [-1, 1.5]
sigma = 1
# Number of observations
size = 500
# Predictor variables
X1 = np.random.randn(size)
X2 = np.random.randn(size) * 0.8
# Outcome variable
Y = alpha + beta[0] * X1 + beta[1] * X2 + np.random.randn(size) * sigma
# Putting our original data into a dataframe (we'll use it later on)
df = pd.DataFrame(
    data = np.array([X1, X2, Y]),
    index=['X1', 'X2', 'Y']
).T

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

[IN]:
with pm.Model() as model_1:
    # Data
    X_1 = pm.Data('X1', X1)
    X_2 = pm.Data('X2', X2)
    # Priors
    alpha = pm.Normal("alpha", mu=0, sigma=10)
    beta = pm.Normal("beta", mu=0, sigma=10, shape=2)
    sigma = pm.HalfNormal("sigma", sigma=10)
    # Likelihood
    mu = alpha + beta[0] * X_1 + beta[1] * X_2
    Y_obs = pm.Normal("Y_obs", mu=mu, sigma=sigma, observed=Y)
    
    # Posterior
    trace = pm.sample(100, return_inferencedata=False, chains=4)
    az.plot_trace(trace)

Во-первых, обратите внимание на формат построения модели, который является стандартным для PyMC3: вы помещаете все (данные, априорные значения, вероятность и апостериор) в оператор «with», где вы называете свою модель (в нашем случае «model_1»). Позже к этой модели можно будет получить доступ.

Опять же, вам может быть интересно, почему эти приоры. В этом конкретном случае для альфа и бета мы используем нормальные распределения с центром вокруг 0, поскольку у нас нет каких-либо предварительных знаний, которые могли бы указать на сильную связь между X и Y. Что касается сигмы = 10 внутри каждого предшествующего, он устанавливает стандартное отклонение для нормального распределения, поэтому чем выше это число, тем больше будет разброс в наших априорных значениях и, следовательно, тем менее информативными они будут (что может быть хорошо, если вы не очень уверены в этих априорных значениях). Распределение HalfNormal для параметра сигма будет сохранять его положительным.

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

Давайте посмотрим на результат нашей модели:

Каждая «строка» на этом изображении представляет один из наших параметров (альфа, бета, сигма), левая часть представляет апостериорные распределения, а правая часть представляет их сходимость во времени.

Возьмите, например, верхний левый график: он показывает оценочное апостериорное распределение для альфа. Но подождите, на самом деле там 4 дистрибутива, а не один! Это потому, что pm.sample () имеет параметр, называемый последовательностями, установленным на 2 или на количество процессоров в вашей системе, в зависимости от того, что больше. В моем случае это было равно 4, но я решил явно установить его для иллюстраций. Цепочки задают количество независимых цепочек, работающих параллельно, что позволяет нам иметь некоторую статистику сходимости.

Вы также можете увидеть понятие цепочек на графиках с правой стороны, которые имеют по одной строке на цепочку. Но что показывают эти линии? Вы видите 100 в pm.sample ()? Это установка количества отрисовываемых сэмплов. Таким образом, эти строки в основном показывают предполагаемое наиболее вероятное значение этого параметра в каждой из 100 выборок для каждой из 4 цепочек. Для альфы вы можете видеть, что она составляет около -0,8 (хотя с небольшим уклоном в сторону более низких значений).

На средних левом и правом графиках у вас есть одинаковая информация для бета-версии 1 и бета-версии 2 (вы можете отличить их по разным цветам).

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

У нас есть апостериорные распределения и их графики, но вам может быть интересно, как с их помощью делать прогнозы. PyMC3 имеет функцию под названием fast_sample_posterior_predictive, которая позволяет нам делать именно это:

[IN]:
with model_1:
pm.set_data({
        'X1': [0.6, 0.02],
        'X2': [-1.3, 0.3]
    })
    y_test = pm.fast_sample_posterior_predictive(trace, samples=100)
print(y_test['Y_obs'].mean(axis=0))
[OUT]:
[-3.27941019 -0.31231568]

Результатом являются прогнозы для двух наблюдений, которые мы включили в наш пример, созданные путем взятия 100 выборок с использованием наших апостериорных распределений.

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

[IN}:
pm.find_MAP(model=model_1)
[OUT]:
{'alpha': array(-0.85492716),
 'beta': array([-1.05603871,  1.48859634]),
 'sigma_log__': array(0.04380838),
 'sigma': array(1.04478214)}

Даже несмотря на то, что PyMC3 может потребоваться некоторое время для запуска, это определенно делает нашу жизнь намного проще. Но что именно он делает?

PyMC3 не работает точно так же, как мы построили наш пример Монти Холла, в том смысле, что он не вычисляет точные апостериорные распределения. Вместо этого он делает умный ход, фактически выборка данных из апостериорных данных для их оценки (вот почему вы видите функцию «выборки» в «апостериорной» части нашего кода).

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

Линейная регрессия с использованием GLM

Теперь давайте посмотрим, как мы можем использовать класс GLM из PyMC3, чтобы получить апостериор:

[IN]:
from pymc3.glm import GLM
# Creating our model
with pm.Model() as model_glm:
    GLM.from_formula('Y ~ X1 + X2', df)
    trace = pm.sample(100)
    pm.traceplot(trace)
    plt.show()
[OUT]:

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

Интерпретация результатов такая же, как и в предыдущем примере.

Классификация

Логистическая регрессия «с нуля»

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

[IN]:
# Creating a binary variable
df['Y_bin']=df['Y']>df['Y'].mean()

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

[IN]:
with pm.Model() as model_log:
    # Priors
    alpha = pm.Normal("alpha", mu=0, sigma=10)
    beta = pm.Normal("beta", mu=0, sigma=10, shape=2)
    # Likelihood    
    p = pm.Deterministic('p', pm.math.sigmoid(alpha + beta[0] * X1 + beta[1] * X2))
    Y_obs = pm.Bernoulli("Y_obs", p, observed=df['Y_bin'])
    
    # Posterior
    trace = pm.sample(100, return_inferencedata=False, chains=4)
    az.plot_trace(trace)
[OUT]:

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

Логистическая регрессия с использованием GLM

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

[IN]:
from pymc3.glm.families import Binomial
with pm.Model() as model_glm_logistic:
    GLM.from_formula('Y_bin ~ X1 + X2', df, family=Binomial())
    trace = pm.sample(100)
    pm.traceplot(trace)
    plt.show()

Разница между этой версией (логистическая регрессия) и первой GLM (линейная регрессия) в основном заключается в параметре family = Binomial (), по крайней мере, с точки зрения построения кода, очевидно.

Заключение

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

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