Bayesian inference

After 3 heads in 3 flips, maximum likelihood estimates the probability of heads as 1 and declares tails impossible. Nobody believes that, because everyone knows something about coins before flipping one. Bayesian inference makes that knowledge part of the calculation. The parameter is treated as a random variable with a prior distribution; the data update it through Bayes’ formula into a posterior; and every question about the parameter, from a point estimate to the probability of the next outcome, is answered from the posterior. This note works through the update, shows that the common regularizers of machine learning are priors, and explains when the Bayesian and the maximum likelihood answers agree.

Notation. $\pi(\theta)$ is the prior density of the parameter, $L(\theta) = p(x \mid \theta)$ the likelihood of the data $x$, $\pi(\theta \mid x)$ the posterior. $\mathrm{Beta}(a, b)$ is the distribution on $[0, 1]$ with density proportional to $p^{a-1} (1-p)^{b-1}$ and mean $a/(a+b)$.

Prior, likelihood, posterior#

Bayes’ formula for densities reads

$$\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 ,$$

or, dropping the normalizing constant that does not depend on $\theta$,

$$\text{posterior} \propto \text{likelihood} \times \text{prior} .$$

The prior says which values of the parameter are plausible before the data; the likelihood, the same function as in maximum likelihood, says how well each value explains the data; the posterior combines them. The denominator $p(x)$, the marginal likelihood or evidence, is the probability of the data averaged over the prior. It plays no role in the shape of the posterior and a central one in comparing models.

Conjugate priors: the beta–binomial#

For a coin with probability of heads $p$, take a $\mathrm{Beta}(a, b)$ prior. After $k$ heads in $n$ flips the likelihood is $p^k (1-p)^{n-k}$, and multiplying,

$$\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) .$$

The posterior is again a beta distribution: the beta family is conjugate to the binomial likelihood. The update is bookkeeping: every head adds 1 to $a$ and every tail adds 1 to $b$. Reading the prior as imagined data has two common conventions, which count from different baselines. In the posterior mean below, $a$ and $b$ act as pseudo-counts: the prior weighs as much as $a + b$ observations, $a$ of them heads. In the density, the exponents are $a - 1$ and $b - 1$, so measured from the uniform $\mathrm{Beta}(1, 1)$, whose exponents are 0, a prior multiplies in the likelihood of $a - 1$ heads and $b - 1$ tails. A $\mathrm{Beta}(2, 2)$ prior, “probably near one half, but not sure”, is therefore worth 4 pseudo-counts in the posterior mean, and is one imagined head and one imagined tail beyond the uniform prior in the density:

Figure

The posterior as the data arrive

The posterior mean

$$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}$$

is a weighted average of the prior mean and the maximum likelihood estimate, with weights given by the number of pseudo-counts and of real observations. After 3 heads in 3 flips it is $5/7 = 0.714$; with the uniform prior it is $(k + 1)/(n + 2) = 0.8$, Laplace’s rule of succession. As $n$ grows the data dominate: after 70 heads in 100 the posterior mean is $0.692$, close to the MLE $0.7$. The same conjugacy holds for the categorical distribution with a symmetric Dirichlet prior, which is exactly add-$\alpha$ smoothing of the n-gram language models seen as a posterior mean.

A credible interval contains the parameter with a stated posterior probability. For $\mathrm{Beta}(72, 32)$ the central 95% interval is $[0.601, 0.777]$, close to the normal-approximation interval around the MLE, $0.7 \pm 1.96 \cdot 0.046 = [0.61, 0.79]$. The interpretations differ: the credible interval is a probability statement about the parameter given these data, the confidence interval of the next note a statement about a procedure over repeated data.

Normal data and shrinkage#

For a normal mean $\mu$ with a normal prior, the posterior is normal too. With prior $N(\mu_0, \tau^2)$ and $n$ observations of $N(\mu, \sigma^2)$ with known $\sigma^2$ and average $\bar x$, precisions (inverse variances) add:

$$\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) .$$

Figure

The posterior mean is a compromise between the prior and the data

The posterior mean is pulled from the data toward the prior mean, shrunk, and it is pulled more when the data are few or noisy. This is the same compromise that the estimators note found can lower the mean squared error: a little bias for a large reduction in variance.

MAP estimation and regularization#

The maximum a posteriori (MAP) estimate is the mode of the posterior:

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

For a linear model with Gaussian noise of variance $\sigma^2$ and a prior $w \sim N(0, \tau^2 I)$ on the weights, the negative log-posterior is

$$\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} ,$$

ridge regression with penalty $\lambda = \sigma^2/\tau^2$: a tighter prior is a stronger penalty. A Laplace prior, with density $\propto e^{-\lvert w_j \rvert / b}$, gives the L1 penalty of the lasso, which sets some weights exactly to zero. Weight decay in neural networks trained by plain gradient descent is the Gaussian prior again. With 20 noisy observations and 10 features in the code below, least squares has a test error of $3.46$, while ridge with $\lambda = 1$ and $\lambda = 3$ gets $1.92$ and $1.75$; with $\lambda = 10$ the prior is too strong and the error rises to $2.72$.

The MAP estimate is a single point and discards the rest of the posterior. It also depends on the parametrization, unlike the posterior itself: the logarithm of the mode of the density of $\theta$ is in general not the mode of the density of $\log\theta$.

Predicting with the posterior#

The quantity of interest is usually not the parameter but the next observation. The posterior predictive distribution averages the model’s predictions over the posterior:

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

For the coin, the probability that the next flip is heads is the posterior mean of $p$: $0.714$ after 3 heads in 3, never 1. Averaging over the parameter instead of plugging in one estimate makes predictions honest about how little the data said; for normal data the predictive interval is wider than the plug-in one, and much wider when $n$ is small. For neural networks, ensembles and Monte Carlo dropout are used as practical stand-ins for this average, but they are uncertainty heuristics rather than posterior samples: independently trained networks are not draws from a posterior, and Monte Carlo dropout is a variational approximation to one only under specific assumptions on the prior and the network (Gal and Ghahramani, 2016).

When the posterior has no formula#

Conjugate pairs are the exception. For a logistic regression, a mixture or a neural network the normalizing integral $p(x)$ has no closed form, and three families of methods take over:

• Sampling. Markov chain Monte Carlo draws samples from the posterior using only the unnormalized product of likelihood and prior.
• Variational inference. Approximate the posterior by the closest member of a tractable family, maximizing a lower bound on $\log p(x)$, the method behind the variational autoencoder.
• The Laplace approximation. Replace the posterior by a normal distribution centred at the MAP estimate with covariance given by the inverse curvature of the log-posterior.

Priors and large samples#

The choice of prior is the most criticized part of the method: two analysts with different priors get different posteriors from the same data. Three answers are standard. Weakly informative priors encode only scale and plausibility (a coin is not two-headed, a weight is not $10^6$). Sensitivity analysis reports how much the conclusions move with the prior. And with enough data the prior stops mattering: under regularity conditions the Bernstein–von Mises theorem says that the posterior approaches the normal distribution $N\big(\hat\theta_{\text{MLE}}, \frac{1}{nI(\hat\theta)}\big)$, the same normal approximation that maximum likelihood gives. For large samples, Bayesian and frequentist intervals nearly coincide; for small samples, with many parameters or with strong prior knowledge, they do not, and the prior is then exactly where the extra information enters.

What changed#

	Maximum likelihood	Bayesian inference
The parameter	fixed, unknown	random, with a prior
Result	a point $\hat\theta$	a distribution $\pi(\theta \mid x)$
3 heads in 3 flips	$\hat p = 1$	posterior mean $0.714$ with a $\mathrm{Beta}(2,2)$ prior
Regularization	added by hand	a prior: L2 = Gaussian, L1 = Laplace
Prediction	plug in $\hat\theta$	average over the posterior
Large samples	$\hat\theta \approx N(\theta, 1/(nI))$	posterior $\approx N(\hat\theta, 1/(nI))$

Limitations#

• The prior is a modelling choice. It can encode real knowledge or real bias; for small data it shapes the answer, and it must be stated and checked.
• Exact posteriors are rare. Outside conjugate families, inference needs sampling or approximation, which are expensive and can fail silently.
• A posterior is only as good as the model. If the likelihood is wrong, the posterior concentrates confidently on the best wrong answer, exactly as maximum likelihood does.

What you should remember#

• Posterior ∝ likelihood × prior; the marginal likelihood $p(x)$ normalizes it.
• Beta prior + binomial data: $\mathrm{Beta}(a + k, b + n - k)$; in the posterior mean the prior acts as $a$ and $b$ pseudo-counts, of total weight $a + b$; the posterior mean is a weighted average of the prior mean and $k/n$.
• Normal prior + normal data: precisions add, the posterior mean is shrunk toward the prior mean.
• MAP = argmax of log-likelihood + log-prior; Gaussian prior = L2 (ridge, weight decay), Laplace prior = L1 (lasso).
• Predict with the posterior predictive, an average over the parameter, not a plug-in.
• With many observations the posterior is approximately $N(\hat\theta_{\text{MLE}}, 1/(nI))$ and the prior washes out.

Code#

The beta–binomial updates and intervals, the normal update, and ridge regression as a MAP estimate:

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

Next#

A credible interval states how probable the parameter values are, given a prior. The frequentist alternative avoids priors and makes a statement about the procedure instead: an interval constructed so that, over repeated samples, it contains the true value 95% of the time. The next note builds such intervals from the central limit theorem and, when no formula is available, from the data alone by the bootstrap.

References#

• A. Gelman, J. B. Carlin, H. S. Stern, D. B. Dunson, A. Vehtari, D. B. Rubin. Bayesian Data Analysis, 3rd ed. CRC Press, 2013. Chapters 2–3 (single- and multi-parameter models, conjugacy), Chapter 4 (normal approximation, Bernstein–von Mises).
• C. M. Bishop. Pattern Recognition and Machine Learning. Springer, 2006. Section 2.1 (the beta distribution), Section 3.3 (Bayesian linear regression), Section 4.4 (the Laplace approximation).
• J. K. Blitzstein, J. Hwang. Introduction to Probability, 2nd ed. CRC Press, 2019. Section 8.3 (beta–binomial conjugacy).
• 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.