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

После 3 орлов в 3 бросках метод максимального правдоподобия оценивает вероятность орла как 1 и объявляет решку невозможной. Никто в это не верит, потому что каждый знает кое-что о монетах ещё до первого броска. Байесовский вывод (Bayesian inference) делает это знание частью вычисления. Параметр считается случайной величиной с априорным распределением (prior); данные обновляют его по формуле Байеса до апостериорного (posterior); и на любой вопрос о параметре, от точечной оценки до вероятности следующего исхода, отвечают по апостериорному распределению. В этой заметке мы проводим это обновление, показываем, что распространённые регуляризаторы машинного обучения — это априорные распределения, и объясняем, когда байесовский ответ совпадает с ответом метода максимального правдоподобия.

Обозначения. $\pi(\theta)$ — априорная плотность параметра, $L(\theta) = p(x \mid \theta)$ — правдоподобие данных $x$, $\pi(\theta \mid x)$ — апостериорная плотность. $\mathrm{Beta}(a, b)$ — распределение на $[0, 1]$ с плотностью, пропорциональной $p^{a-1} (1-p)^{b-1}$, и средним $a/(a+b)$.

Априорное распределение, правдоподобие, апостериорное распределение#

Формула Байеса для плотностей имеет вид

$$\pi(\theta \mid x) = \frac{p(x \mid \theta)\, \pi(\theta)}{p(x)}, \qquad p(x) = \int p(x \mid \theta)\, \pi(\theta)\, d\theta ,$$

или, если отбросить нормирующую константу, которая не зависит от $\theta$,

$$\text{апостериорное} \propto \text{правдоподобие} \times \text{априорное} .$$

Априорное распределение говорит, какие значения параметра вероятны ещё до данных; правдоподобие — та же функция, что и в методе максимального правдоподобия, — говорит, насколько хорошо каждое значение объясняет данные; апостериорное распределение объединяет их. Знаменатель $p(x)$, маргинальное правдоподобие (marginal likelihood), или обоснованность (evidence), — это вероятность данных, усреднённая по априорному распределению. На форму апостериорного распределения он не влияет, но играет центральную роль при сравнении моделей.

Сопряжённые априорные распределения: бета-биномиальная модель#

Для монеты с вероятностью орла $p$ возьмём априорное распределение $\mathrm{Beta}(a, b)$. После $k$ орлов в $n$ бросках правдоподобие равно $p^k (1-p)^{n-k}$, и после перемножения

$$\pi(p \mid k) \propto p^{a + k - 1} (1-p)^{b + n - k - 1} , \qquad p \mid k \sim \mathrm{Beta}(a + k,\ b + n - k) .$$

Апостериорное распределение — снова бета-распределение: бета-семейство сопряжено (conjugate) с биномиальным правдоподобием. Обновление сводится к подсчёту: каждый орёл прибавляет 1 к $a$, каждая решка — 1 к $b$. Прочитать априорное распределение как воображаемые данные можно двумя распространёнными способами, и они ведут отсчёт от разных точек. В апостериорном среднем ниже $a$ и $b$ играют роль псевдосчётчиков (pseudo-counts): априорное распределение весит столько же, сколько $a + b$ наблюдений, из которых $a$ — орлы. В плотности показатели степени равны $a - 1$ и $b - 1$, так что относительно равномерного $\mathrm{Beta}(1, 1)$, у которого они нулевые, априорное распределение домножает на правдоподобие $a - 1$ орлов и $b - 1$ решек. Поэтому априорное $\mathrm{Beta}(2, 2)$ — «вероятно, около половины, но это не точно» — весит 4 псевдосчётчика в апостериорном среднем, а в плотности это один воображаемый орёл и одна воображаемая решка сверх равномерного априорного распределения:

Рисунок

Апостериорное распределение по мере поступления данных

Апостериорное среднее

$$E[p \mid k] = \frac{a + k}{a + b + n} = \frac{a + b}{a + b + n} \cdot \frac{a}{a+b} + \frac{n}{a + b + n} \cdot \frac{k}{n}$$

— взвешенное среднее априорного среднего и оценки максимального правдоподобия (MLE) с весами, которые задаются числом псевдосчётчиков и числом реальных наблюдений. После 3 орлов в 3 бросках оно равно $5/7 = 0.714$; при равномерном априорном распределении — $(k + 1)/(n + 2) = 0.8$, это правило следования Лапласа (Laplace’s rule of succession). С ростом $n$ доминируют данные: после 70 орлов в 100 бросках апостериорное среднее равно $0.692$, что близко к оценке максимального правдоподобия $0.7$. То же сопряжение выполняется для категориального распределения с симметричным априорным распределением Дирихле, и аддитивное сглаживание add-$\alpha$ n-граммных языковых моделей — это в точности апостериорное среднее.

Байесовский доверительный интервал (credible interval) содержит параметр с заданной апостериорной вероятностью. Для $\mathrm{Beta}(72, 32)$ центральный интервал с вероятностью 95% равен $[0.601, 0.777]$, что близко к интервалу нормального приближения вокруг оценки максимального правдоподобия, $0.7 \pm 1.96 \cdot 0.046 = [0.61, 0.79]$. Интерпретации различаются: байесовский доверительный интервал — вероятностное утверждение о параметре при данных наблюдениях, а доверительный интервал из следующей заметки — утверждение о процедуре при повторных выборках.

Нормальные данные и сжатие#

Для среднего $\mu$ нормального распределения с нормальным априорным распределением апостериорное распределение тоже нормально. При априорном распределении $N(\mu_0, \tau^2)$ и $n$ наблюдениях из $N(\mu, \sigma^2)$ с известной дисперсией $\sigma^2$ и выборочным средним $\bar x$ обратные дисперсии (precisions) складываются:

$$\frac{1}{\sigma_n^2} = \frac{1}{\tau^2} + \frac{n}{\sigma^2}, \qquad \mu_n = \sigma_n^2 \Big(\frac{\mu_0}{\tau^2} + \frac{n \bar x}{\sigma^2}\Big) .$$

Рисунок

Апостериорное среднее — компромисс между априорным и данными

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

MAP-оценка и регуляризация#

MAP-оценка (maximum a posteriori, оценка апостериорного максимума) — мода апостериорного распределения:

$$\hat\theta_{\text{MAP}} = \arg\max_\theta \big(\log p(x \mid \theta) + \log \pi(\theta)\big) .$$

Для линейной модели с гауссовским шумом дисперсии $\sigma^2$ и априорным распределением $w \sim N(0, \tau^2 I)$ на весах отрицательный логарифм апостериорной плотности равен

$$\frac{1}{2\sigma^2} \sum_i (y_i - w^\top x_i)^2 + \frac{1}{2\tau^2} \lVert w \rVert^2 + \text{const} ,$$

а это гребневая регрессия со штрафом $\lambda = \sigma^2/\tau^2$: чем уже априорное распределение, тем сильнее штраф. Априорное распределение Лапласа с плотностью $\propto e^{-\lvert w_j \rvert / b}$ даёт штраф L1 лассо, который обращает некоторые веса в точности в ноль. Затухание весов в нейросетях при обычном градиентном спуске — снова гауссовское априорное распределение. На 20 зашумлённых наблюдениях с 10 признаками из кода ниже метод наименьших квадратов даёт ошибку на тесте $3.46$, а гребневая регрессия с $\lambda = 1$ и $\lambda = 3$ — $1.92$ и $1.75$; при $\lambda = 10$ априорное распределение слишком сильное, и ошибка вырастает до $2.72$.

MAP-оценка — одна точка, и остальное апостериорное распределение она отбрасывает. Кроме того, в отличие от самого апостериорного распределения, она зависит от параметризации: логарифм моды плотности $\theta$ в общем случае не равен моде плотности $\log\theta$.

Предсказание по апостериорному распределению#

Обычно интересен не сам параметр, а следующее наблюдение. Апостериорное предсказательное распределение (posterior predictive) усредняет предсказания модели по апостериорному распределению:

$$p(x_{\text{new}} \mid x) = \int p(x_{\text{new}} \mid \theta)\, \pi(\theta \mid x)\, d\theta .$$

Для монеты вероятность того, что следующий бросок даст орла, равна апостериорному среднему $p$: $0.714$ после 3 орлов в 3 бросках и никогда не 1. Усреднение по параметру вместо подстановки одной оценки делает предсказания честными в том, как мало сказали данные; для нормальных данных предсказательный интервал шире интервала с подставленной оценкой (plug-in), и намного шире при малом $n$. Для нейросетей вместо этого усреднения на практике используют ансамбли и дропаут Монте-Карло (Monte Carlo dropout), но это эвристики для оценки неопределённости, а не сэмплы из апостериорного распределения: независимо обученные сети — не выборка из апостериорного распределения, а дропаут Монте-Карло — вариационное приближение к нему лишь при определённых предположениях об априорном распределении и о сети (Gal and Ghahramani, 2016).

Когда у апостериорного распределения нет формулы#

Сопряжённые пары — исключение. Для логистической регрессии, смеси распределений или нейросети нормирующий интеграл $p(x)$ не выражается в замкнутом виде, и в дело вступают три семейства методов:

• Сэмплирование. Марковские методы Монте-Карло (Markov chain Monte Carlo, MCMC) порождают сэмплы из апостериорного распределения, используя только ненормированное произведение правдоподобия и априорной плотности.
• Вариационный вывод (variational inference). Апостериорное распределение приближается ближайшим членом удобного для вычислений семейства через максимизацию нижней оценки $\log p(x)$; на этом методе основан вариационный автоэнкодер.
• Приближение Лапласа (Laplace approximation). Апостериорное распределение заменяется нормальным с центром в MAP-оценке и ковариацией, обратной к кривизне логарифма апостериорной плотности.

Априорные распределения и большие выборки#

Выбор априорного распределения — самая критикуемая часть метода: два аналитика с разными априорными распределениями получают по одним и тем же данным разные апостериорные. Стандартных ответов три. Слабоинформативные априорные распределения (weakly informative priors) кодируют только масштаб и разумные пределы значений (у монеты не два орла, вес не бывает порядка $10^6$). Анализ чувствительности показывает, насколько заключения меняются вместе с априорным распределением. И наконец, при достаточном объёме данных априорное распределение перестаёт иметь значение: при условиях регулярности теорема Бернштейна–фон Мизеса (Bernstein–von Mises theorem) утверждает, что апостериорное распределение приближается к нормальному распределению $N\big(\hat\theta_{\text{MLE}}, \frac{1}{nI(\hat\theta)}\big)$ — тому же нормальному приближению, которое даёт метод максимального правдоподобия. На больших выборках байесовские и частотные (frequentist) интервалы почти совпадают; на малых выборках, при большом числе параметров или при сильном априорном знании — нет, и тогда дополнительная информация входит именно через априорное распределение.

Что изменилось#

	Метод максимального правдоподобия	Байесовский вывод
Параметр	фиксированный, неизвестный	случайный, с априорным распределением
Результат	точка $\hat\theta$	распределение $\pi(\theta \mid x)$
3 орла в 3 бросках	$\hat p = 1$	апостериорное среднее $0.714$ при априорном $\mathrm{Beta}(2,2)$
Регуляризация	добавляется вручную	априорное распределение: L2 = гауссовское, L1 = лапласовское
Предсказание	подставить $\hat\theta$	усреднить по апостериорному распределению
Большие выборки	$\hat\theta \approx N(\theta, 1/(nI))$	апостериорное $\approx N(\hat\theta, 1/(nI))$

Ограничения#

• Априорное распределение — модельное допущение. Оно может кодировать реальное знание или реальную предвзятость; на малых данных оно определяет ответ, и его нужно явно указывать и проверять.
• Точные апостериорные распределения редки. Вне сопряжённых семейств вывод требует сэмплирования или приближений, которые дороги и могут незаметно давать сбой.
• Апостериорное распределение хорошо ровно настолько, насколько хороша модель. Если правдоподобие неверно, апостериорное распределение уверенно концентрируется на наилучшем неверном ответе — точно так же, как метод максимального правдоподобия.

Что нужно запомнить#

• Апостериорная плотность ∝ правдоподобие × априорная плотность; маргинальное правдоподобие $p(x)$ её нормирует.
• Бета-априорное распределение + биномиальные данные: $\mathrm{Beta}(a + k, b + n - k)$; в апостериорном среднем априорное распределение действует как $a$ и $b$ псевдосчётчиков общим весом $a + b$; апостериорное среднее — взвешенное среднее априорного среднего и $k/n$.
• Нормальное априорное распределение + нормальные данные: обратные дисперсии складываются, апостериорное среднее сжимается к априорному среднему.
• MAP-оценка = argmax суммы логарифма правдоподобия и логарифма априорной плотности; гауссовское априорное распределение = L2 (гребневая регрессия, затухание весов), лапласовское = L1 (лассо).
• Предсказывать нужно по апостериорному предсказательному распределению — среднему по параметру, а не подстановкой оценки.
• При большом числе наблюдений апостериорное распределение приближённо равно $N(\hat\theta_{\text{MLE}}, 1/(nI))$, и влияние априорного распределения исчезает.

Код#

Бета-биномиальные обновления и интервалы, обновление для нормального распределения и гребневая регрессия как MAP-оценка:

import numpy as np
from scipy import stats, optimize

# the coin with a Beta(2, 2) prior: posterior means and 95% credible intervals as the flips arrive
for a, b in ((2, 2), (5, 2), (9, 5), (72, 32)):
    post = stats.beta(a, b)
    print((a, b), round(post.mean(), 3), post.ppf([0.025, 0.975]).round(3))
# (2, 2) 0.5 [0.094 0.906] | (5, 2) 0.714 [0.359 0.957] | (9, 5) 0.643 [0.386 0.861] | (72, 32) 0.692 [0.601 0.777]
print(round(0.7 - 1.96 * 0.0458, 3), round(0.7 + 1.96 * 0.0458, 3))   # 0.61 0.79: the normal interval around the MLE

# normal mean with prior N(0, 1), four observations with mean 2 and noise variance 1: precisions add
precision = 1 / 1 + 4 / 1
print(round(4 * 2 / precision, 2), round(1 / precision, 2))          # 1.6 0.2

# MAP with a Gaussian prior on the weights is ridge regression
rng = np.random.default_rng(0)
X, w = rng.normal(size=(20, 10)), rng.normal(size=10)
y = X @ w + rng.normal(size=20)
ridge = np.linalg.solve(X.T @ X + np.eye(10), X.T @ y)
neg_log_post = lambda v: 0.5 * ((y - X @ v) ** 2).sum() + 0.5 * (v ** 2).sum()   # noise variance 1, prior N(0, I)
print(float(np.abs(ridge - optimize.minimize(neg_log_post, np.zeros(10)).x).max()) < 1e-5)   # True

# on new data the prior helps: test error of least squares (lambda = 0) and ridge
Xt = rng.normal(size=(10_000, 10))
yt = Xt @ w + rng.normal(size=10_000)
for lam in (0.0, 1.0, 3.0, 10.0):
    wl = np.linalg.solve(X.T @ X + lam * np.eye(10), X.T @ y)
    print(lam, round(float(((yt - Xt @ wl) ** 2).mean()), 3))
# 0.0 3.463 | 1.0 1.917 | 3.0 1.746 | 10.0 2.719

Дальше#

Байесовский доверительный интервал говорит, насколько вероятны значения параметра при заданном априорном распределении. Частотная альтернатива обходится без априорных распределений и делает утверждение о процедуре: интервал строится так, чтобы на повторных выборках он содержал истинное значение в 95% случаев. Следующая заметка строит такие интервалы по центральной предельной теореме, а когда формулы нет — по одним только данным с помощью бутстрэпа.

Литература#

• A. Gelman, J. B. Carlin, H. S. Stern, D. B. Dunson, A. Vehtari, D. B. Rubin. Bayesian Data Analysis, 3rd ed. CRC Press, 2013. Главы 2–3 (модели с одним и несколькими параметрами, сопряжённость), глава 4 (нормальное приближение, теорема Бернштейна–фон Мизеса).
• C. M. Bishop. Pattern Recognition and Machine Learning. Springer, 2006. Раздел 2.1 (бета-распределение), раздел 3.3 (байесовская линейная регрессия), раздел 4.4 (приближение Лапласа).
• J. K. Blitzstein, J. Hwang. Introduction to Probability, 2nd ed. CRC Press, 2019. Раздел 8.3 (бета-биномиальная сопряжённость).
• R. Tibshirani. Regression Shrinkage and Selection via the Lasso. Journal of the Royal Statistical Society B, 58(1), 1996.
• Y. Gal, Z. Ghahramani. Dropout as a Bayesian Approximation: Representing Model Uncertainty in Deep Learning. ICML, 2016.