Bayesian reasoning, uncertainty made quantitative.

In the frequentist chapter we estimated a parameter $\theta$ by picking the value that made the data most likely: the maximum-likelihood estimator. The parameter itself had no distribution; it was a fixed unknown, and our estimator was the random thing, re-rolled under hypothetical replications of the experiment. That framing is excellent for repeated, controllable experiments. It is awkward almost everywhere else.

The Bayesian alternative is to treat $\theta$ itself as a random variable — not because it is physically random, but because we are uncertain about it. Before seeing data, we summarise that uncertainty with a prior distribution $p(\theta)$. The data $y$ arrive through a likelihood $p(y \mid \theta)$. Bayes' theorem combines them into a posterior $p(\theta \mid y)$, which is our new, updated description of uncertainty. That is the whole story. Everything else — priors, conjugacy, MCMC, variational inference, Bayesian neural networks — is machinery for doing that update in cases where the integrals get hard.

The payoff is large. Bayesian reasoning handles small samples, hierarchical structure, and missing data gracefully; it produces probabilistic predictions rather than point estimates; it gives honest uncertainty estimates for downstream decisions; and it connects cleanly to information theory, decision theory, and many modern deep-learning ideas — variational autoencoders, Thompson sampling, Bayesian optimisation, active learning, and the uncertainty quantification that production ML systems increasingly require.

The cost is that the integrals are usually intractable, which is why so much of this chapter is about how to do Bayesian inference in practice. The conceptual apparatus is small; the computational ingenuity is where the subject lives.

Start from the definition of conditional probability. For events $A$ and $B$ with $p(B) > 0$,

Equating the two expressions for the joint $p(A, B)$ gives Bayes' theorem:

Applied to inference — replace $A$ with the parameter $\theta$ and $B$ with the observed data $y$ — we get the form every Bayesian argument starts from:

The denominator $p(y) = \int p(y \mid \theta)\, p(\theta)\, d\theta$ is the marginal likelihood or evidence. It is a constant with respect to $\theta$, so for most purposes we work with the unnormalised form

reading the symbol "$\propto$" as "equal up to a normalising constant we will fix later, if we need to." The proportionality is often all you need: it is enough to decide whether $\theta_1$ or $\theta_2$ is more probable, to find the mode, or to run an MCMC sampler.

The disease-testing example

A textbook classic, worth doing once: a disease has prevalence $0.1\%$ in the population. A test has sensitivity $99\%$ (true positive rate) and specificity $99\%$ (true negative rate). You test positive. What is the probability you have the disease?

About $9\%$. The posterior is dominated by the prior: because the disease is rare, almost all positives are false positives. This is the canonical demonstration that Bayes' theorem is not intuitive — human reasoning routinely ignores base rates — and that the normalising constant in the denominator matters a great deal when the prior is strong.

A prior $p(\theta)$ is any probability distribution over the parameter space. You can pick it to reflect genuine knowledge ("previous experiments have given us roughly a Gaussian centred on $0.3$"), to encode mild regularisation ("the parameter is probably small"), or to be as uninformative as possible. Different choices lead to different posteriors — often similar when data are plentiful, sometimes wildly different when data are scarce.

A rough taxonomy

Informative priors use domain knowledge: a Normal$(0.3, 0.05^2)$ prior on a conversion rate because the last six A/B tests landed there. These are appropriate when you have real prior information and are willing to defend it.

Weakly informative priors encode "common sense" without strong commitment — a Normal$(0, 10^2)$ on a regression coefficient, say, which rules out physically impossible magnitudes without preferring any plausible value. The Stan community, following Gelman, recommends these as the default in applied work.

Non-informative priors try to let the data speak for themselves. The uniform prior $p(\theta) \propto 1$ is the simplest, but it is not invariant under reparameterisation — a uniform prior on $\theta$ is not a uniform prior on $\log \theta$. Jeffreys priors fix this by defining $p(\theta) \propto \sqrt{\det F(\theta)}$ where $F$ is the Fisher information matrix. Jeffreys priors are invariant under smooth reparameterisation, which is the mathematically principled non-informative choice.

Improper priors — priors that do not integrate to one, such as $p(\theta) \propto 1$ on an unbounded parameter — are often fine in practice as long as the resulting posterior is proper. But they can lead to paradoxes in model comparison, and are best avoided in hierarchical models.

Empirical Bayes estimates the prior from the data itself — for example, fitting a Gaussian to the observed variation across groups and using it as a prior on each group's parameter. This is philosophically awkward (you are using the data twice) but computationally convenient; in the hierarchical section we will see it emerge naturally as an approximation to the fully Bayesian approach.

In practice, the question "what prior should I use?" is usually better phrased as "what does my model say before seeing data, and does that look reasonable?" Simulating from the prior predictive distribution (Section 14) makes this question concrete and answerable.

The full Bayesian answer to "what have I learned about $\theta$?" is the distribution $p(\theta \mid y)$. In practice you usually want scalar summaries. The standard choices:

The posterior mean $\mathbb{E}[\theta \mid y] = \int \theta\, p(\theta \mid y)\, d\theta$ minimises squared-error loss. It is the Bayes estimator under quadratic loss, and the most common "best single number" summary.

The posterior median minimises absolute-error loss; it is more robust to skewness and heavy tails than the mean.

The posterior mode, or maximum a posteriori (MAP) estimate, maximises $p(\theta \mid y)$. It is particularly convenient because for point estimation it only requires optimisation, not integration:

This is maximum likelihood plus a log-prior penalty — the same object we called "regularised MLE" in the statistics chapter. MAP is fast and often good but it loses the one thing that makes Bayesian inference distinctive: a full distribution. A thin, sharply peaked posterior and a wide, flat one can share the same mode.

Credible intervals

A $95\%$ credible interval is any interval $[a, b]$ with $\int_a^b p(\theta \mid y)\, d\theta = 0.95$. It has a satisfyingly direct interpretation: "given the model and data, there is a $95\%$ probability that $\theta$ is between $a$ and $b$." This is what most people think a frequentist confidence interval means — and what only a Bayesian credible interval actually means.

Two common choices fix the arbitrariness of $[a, b]$. The equal-tailed interval uses the $2.5\%$ and $97.5\%$ posterior quantiles; it is easy to compute from MCMC samples but can be misleading when the posterior is skewed. The highest posterior density (HPD) interval is the shortest interval containing $95\%$ of the mass, and always contains the mode. HPD is the principled choice; equal-tailed is what everyone actually reports.

The posterior predictive distribution

Predictions about future data $\tilde{y}$ marginalise out the parameter:

This is the posterior predictive distribution. It incorporates both observational noise (through $p(\tilde{y} \mid \theta)$) and parameter uncertainty (through $p(\theta \mid y)$) into a single forecast. A frequentist plug-in prediction $p(\tilde{y} \mid \hat{\theta})$ uses only the first; it is overconfident by exactly the variance that averaging over $\theta$ would have added.

A prior $p(\theta)$ is conjugate to a likelihood $p(y \mid \theta)$ if the posterior $p(\theta \mid y)$ is in the same parametric family as the prior. The practical payoff is enormous: there is no integral to do. You just update the parameters of the prior using a simple rule.

Beta–Binomial

The classic example. Observe $y$ successes in $n$ trials with unknown success probability $\theta$. The likelihood is $p(y \mid \theta) = \binom{n}{y} \theta^y (1-\theta)^{n-y}$. With a Beta prior $\theta \sim \mathrm{Beta}(\alpha, \beta)$,

which is $\mathrm{Beta}(\alpha + y,\, \beta + n - y)$. The update is "add successes to $\alpha$, failures to $\beta$" — exactly what intuition says it should be. The prior parameters $\alpha$ and $\beta$ act as pseudo-counts: they behave like $\alpha$ previously observed successes and $\beta$ failures.

Other standard pairs

Every exponential-family likelihood has a natural conjugate prior, and the update is a linear operation on natural parameters. This is more than a coincidence: conjugacy is the exponential-family's way of being well-behaved under Bayesian updating, and it is the reason so many classical Bayesian results have closed forms.

The frequentist model is $y = X\boldsymbol{\beta} + \boldsymbol{\varepsilon}$ with $\boldsymbol{\varepsilon} \sim \mathcal{N}(\mathbf{0}, \sigma^2 I)$. The Bayesian version adds a prior on $\boldsymbol{\beta}$. Take a Gaussian prior $\boldsymbol{\beta} \sim \mathcal{N}(\mathbf{0}, \tau^2 I)$ for concreteness and assume $\sigma^2$ is known.

The posterior is Gaussian (conjugate!), with mean and covariance

The posterior mean is the ridge regression estimator: $(X^\top X + \lambda I)^{-1} X^\top y$ with $\lambda = \sigma^2/\tau^2$. This is the exact correspondence we saw in the statistics chapter, now with the prior showing where the penalty came from. A Laplace prior would have given the lasso. Classical regularisation is Bayesian MAP estimation; the hyperparameter is the prior variance.

Predictive distribution

For a new input $\mathbf{x}^*$, the posterior predictive is

The predictive variance has two components. The $\sigma^2$ term is the irreducible observation noise — you cannot predict it away. The $\mathbf{x}^{*\top}\Sigma_n \mathbf{x}^*$ term is the uncertainty in $\boldsymbol{\beta}$, propagated through the new input. It is largest in regions where the training data are sparse or the directions of $\mathbf{x}^*$ are poorly constrained by $X^\top X$. This is the Bayesian answer to "how confident am I at this input?" — and one reason Bayesian regression is still the first stop for small-data problems with meaningful uncertainty requirements.

You have measurements of some quantity in $J$ groups: conversion rates at $J$ stores, test scores at $J$ schools, batting averages for $J$ players. Two extreme analyses present themselves. No pooling: estimate each group's parameter $\theta_j$ from its own data, treat them as unrelated. Complete pooling: assume $\theta_1 = \cdots = \theta_J$, pool all the data, and estimate one number. The first over-fits; the second under-fits.

A hierarchical (or multilevel) model sits between them. Each group has its own $\theta_j$, but the $\theta_j$ themselves come from a shared distribution:

with priors on $\mu$ and $\tau$. The hyperparameters $\mu$ and $\tau$ are learned from the data just like any other parameters. The posterior for each $\theta_j$ is then pulled — shrunk — toward the overall mean $\mu$, by an amount that depends on how much data that group has and how spread out the groups appear to be.

The 8-schools problem

Gelman's canonical example: an SAT coaching programme is evaluated separately at $8$ schools. Each school reports a point estimate $\hat{\theta}_j$ and standard error $\sigma_j$. The no-pooling analysis is noisy; the complete-pooling analysis washes out any real differences. The hierarchical analysis shrinks each school's estimate toward the grand mean by a factor depending on $\tau$, producing estimates that are simultaneously more accurate (on average) and more honest about uncertainty than either extreme.

Hierarchical models are the most practically useful Bayesian technique — especially when data are unbalanced across groups, when some groups have very little data, or when borrowing strength across groups is itself the modelling goal. Gelman & Hill's Data Analysis Using Regression and Multilevel/Hierarchical Models is essentially a whole book about the technique, and deservedly so.

Centred vs. non-centred parameterisation

A numerical footnote that matters in practice: the obvious parameterisation $\theta_j \sim \mathcal{N}(\mu, \tau^2)$ (called centred) produces posteriors that are often pathologically hard for samplers to explore when $\tau$ is small. The non-centred form writes $\theta_j = \mu + \tau z_j$ with $z_j \sim \mathcal{N}(0, 1)$, and samplers handle it much better. Any serious Stan or PyMC model uses non-centred by default.

Suppose two models, $M_1$ and $M_2$, each with its own parameters $\theta_1$ and $\theta_2$. Bayes' theorem at the model level gives posterior model probabilities

and the ratio between two models is the Bayes factor,

Each numerator and denominator is a marginal likelihood — an integral over the full parameter space — which automatically penalises model complexity through the volume of parameter space assigned prior mass that ends up with low likelihood. In this sense Bayes factors are the principled Bayesian version of "Occam's razor": simpler models that still fit the data win.

Why it often does not work

Bayes factors are extraordinarily sensitive to prior choice. Doubling the prior variance on a parameter halves the prior density at the MLE, halves the marginal likelihood, and halves the Bayes factor. This is unproblematic when priors are strongly informed; it is a disaster when priors are vague. For non-informative priors — the exact situation in which you might want a "neutral" comparison — Bayes factors can be arbitrary.

The Bayesian Information Criterion (BIC) is a large-sample approximation,

with $k$ parameters and $n$ observations. Lower is better; differences of $2$–$6$ indicate positive evidence, $6$–$10$ strong, $>10$ very strong (Raftery's scale). BIC is a rough approximation to $-2\log p(y \mid M)$ but comes for free once you have the MLE. It is widely used and widely criticised.

What modern practitioners use instead

Increasingly, Bayesian model comparison reaches for predictive criteria that compare models on how well they would forecast held-out data:

• WAIC (widely applicable information criterion): approximates leave-one-out predictive log-likelihood, works with most models including those without point estimates.
• PSIS-LOO (Pareto-smoothed importance-sampling leave-one-out): a more accurate and robust version of WAIC; now the default in PyMC, Stan, and brms.
• Posterior predictive checks (Section 14): rather than choosing between models, test whether a given model can generate data that look like the real thing.

These are less theoretically pure than Bayes factors but vastly more usable. They are the tools you will actually reach for in practice.

A sequence $(Y_1, Y_2, \dots)$ is exchangeable if its joint distribution is invariant under permutations — the order of the observations tells us nothing. This is a weaker assumption than i.i.d., which additionally requires independence.

De Finetti's theorem says that every infinitely exchangeable sequence can be written as conditionally i.i.d. given some (possibly infinite-dimensional) parameter $\theta$ drawn from a distribution $\pi$:

Read this carefully. The left side is a joint distribution over observations — no explicit parameter. The right side is the Bayesian decomposition — prior, likelihood, marginal. The theorem says that if you believe in exchangeability, you are already a Bayesian, even if you did not know it. A prior and a likelihood are not optional modelling choices; they are the unique representation of any exchangeable sequence.

The theorem generalises. Partial exchangeability (order within groups does not matter, between groups does) gives hierarchical models. Markov exchangeability gives latent-state models. Every time you make a modelling choice, there is a symmetry assumption behind it, and a de Finetti-style representation that corresponds.

The goal: draw samples $\theta^{(1)}, \theta^{(2)}, \dots$ whose long-run distribution is the posterior $p(\theta \mid y)$. Because we only need the unnormalised density $p(y \mid \theta)\, p(\theta)$ — the constant of proportionality cancels — we can sample from any posterior we can evaluate up to a constant, which is essentially all of them.

Metropolis-Hastings

The simplest general MCMC algorithm. Given a current state $\theta$, propose a new state $\theta'$ from a proposal density $q(\theta' \mid \theta)$. Accept it with probability

Otherwise stay at $\theta$. Under mild conditions, the sequence of states converges to the posterior. A symmetric random-walk proposal $q(\theta' \mid \theta) = q(\theta \mid \theta')$ simplifies the ratio to the posterior ratio alone.

Metropolis-Hastings always works in principle. In practice, high-dimensional posteriors force impossibly small step sizes (or exponentially low acceptance rates), and mixing becomes glacial.

Gibbs sampling

A special case: cycle through parameters one at a time, sampling each from its full conditional $p(\theta_j \mid \theta_{-j}, y)$. When full conditionals are conjugate (the usual case in hierarchical models), each step is an exact draw and acceptance is trivially 1. Gibbs sampling was the original BUGS/JAGS algorithm and is still the right tool for conditionally conjugate models.

Hamiltonian Monte Carlo and NUTS

Modern Bayesian computation runs on Hamiltonian Monte Carlo (HMC) and its self-tuning variant, the No-U-Turn Sampler (NUTS). HMC treats the negative log-posterior as a potential energy, augments the parameters with momentum variables drawn from a Gaussian, and simulates Hamiltonian dynamics to propose moves that travel long distances through parameter space while still respecting the target distribution.

The upshot: HMC/NUTS scales to posteriors with thousands of dimensions, handles strong correlation structure naturally, and requires only gradients of the log-posterior — computed automatically by modern PPLs (Stan, PyMC, Numpyro, Turing.jl) via reverse-mode autodiff, the same technique that trains neural networks.

Variational inference turns posterior inference into optimisation. Pick a family $\mathcal{Q}$ of tractable distributions $q(\theta; \phi)$ parameterised by $\phi$, and find the member closest to the true posterior in KL divergence:

The true posterior is unknown, so we cannot compute the KL directly. But a little algebra shows that minimising this KL is equivalent to maximising a lower bound on the log marginal likelihood — the Evidence Lower BOund, or ELBO:

The two terms have separate interpretations. The first is the expected log-joint, which we want to be large. The second is the negative entropy of $q$, which pulls $q$ toward being broad (and so toward covering more of the posterior mass). Together they balance fit against regularisation.

Mean-field and beyond

The simplest family assumes $q$ factorises across parameters: $q(\theta) = \prod_j q_j(\theta_j)$. This is mean-field VI, and it leads to closed-form coordinate updates that resemble Gibbs sampling but optimise rather than sample. Mean-field chronically underestimates posterior variance because it cannot capture correlations — a warning that is worth taking seriously in practice.

Richer families fix this: structured variational families, normalising flows as $q$, amortised VI where $\phi$ is the output of a neural network, and hybrid approaches that mix VI and MCMC.

The Laplace approximation fits a Gaussian to the posterior by matching the mode and local curvature. Around the MAP estimate $\hat{\theta}$, a second-order Taylor expansion of the log-posterior gives

where $H = -\nabla^2 \log p(\theta \mid y)\big|_{\hat{\theta}}$ is the Hessian of the negative log-posterior. Exponentiating gives

That is the whole method. Run an optimiser to find $\hat{\theta}$; compute the Hessian once at $\hat{\theta}$; report a Gaussian. It reproduces the posterior exactly when the log-posterior is quadratic (as in conjugate Gaussian cases), and is surprisingly accurate when the posterior is unimodal and approximately Gaussian near the mode — which, by the Bernstein-von Mises theorem, it is in large samples for regular models.

Connection to Fisher information

As $n \to \infty$, the Hessian of the negative log-posterior converges to $n F(\theta_0)$, where $F$ is the Fisher information matrix. The Laplace posterior becomes $\mathcal{N}(\hat{\theta}, (n F)^{-1})$ — exactly the asymptotic sampling distribution of the MLE. This is the Bernstein-von Mises theorem in practice: the Bayesian posterior and the frequentist sampling distribution agree in the limit, both governed by Fisher information.

A prior predictive check draws parameters from the prior, simulates data given those parameters, and plots the result. If your prior implies that a conversion rate of $0.7$ or $-0.4$ is routine, you will see that immediately — and fix the prior before the model ever touches data.

A posterior predictive check does the same thing with the fitted posterior: draw $\theta$ from $p(\theta \mid y)$, simulate $\tilde{y}$, and compare its distribution — means, medians, extremes, group-wise summaries, whatever is domain-appropriate — to the observed data.

The comparison does not need to be formal. A scatter of simulated versus observed group means, overlaid histograms of simulated versus observed values, a quantile-quantile plot — any of these will catch gross misspecification. Gelman advocates "model criticism" as the third leg of Bayesian workflow (alongside prior and posterior inference), and in applied work it is where most real debugging happens.

Suppose you must choose an action $a$ from a space $\mathcal{A}$, and the goodness of action $a$ given the true parameter $\theta$ is measured by a loss function $L(\theta, a)$. The posterior expected loss, or risk, of action $a$ is

The Bayes action is $a^* = \arg\min_a \rho(a)$. Different losses give different optimal actions:

• Squared error $L(\theta, a) = (\theta - a)^2$: Bayes action is the posterior mean.
• Absolute error $L(\theta, a) = |\theta - a|$: Bayes action is the posterior median.
• 0-1 loss $L(\theta, a) = \mathbb{1}[\theta \neq a]$: Bayes action is the posterior mode (MAP).
• Asymmetric losses (over-prediction costs $c_+$, under-prediction costs $c_-$): Bayes action is the $c_-/(c_-+c_+)$ quantile.

This last case is the correct answer to "what should I predict if false positives and false negatives cost different amounts?" — and is a clean example of how Bayesian decision theory connects posterior inference to downstream economics.

The language of decision theory — loss functions, risk, admissibility, minimax — is also the native language of classical statistics, machine learning, and reinforcement learning. Bayesian decision theory is the bridge: it shows that Bayes estimators are optimal under any loss, that Bayesian decision rules are admissible, and that frequentist minimax procedures often have a Bayesian justification under some (often improper) prior.

The Bayesian influence on ML is broader than "Bayesian ML" would suggest. Here is a sampler of where the ideas in this chapter actually land in practice.

• Regularisation is MAP estimation. Weight decay is a Gaussian prior; L1 is a Laplace prior; dropout can be interpreted as a variational approximation. The loss function of a neural network has a Bayesian reading even if you never invoked it.
• Variational autoencoders. The ELBO is the variational lower bound; the encoder is an amortised variational posterior; the reparameterisation trick makes the expectation differentiable. A VAE is a Bayesian model trained by stochastic VI.
• Gaussian processes. A distribution over functions, defined by a mean and covariance kernel. The posterior over function values conditional on training data is a closed-form Gaussian — the purest Bayesian regression method, still the best tool for small-data regression with calibrated uncertainty.
• Bayesian neural networks. A prior on weights, a likelihood from the architecture, and a posterior approximated by Laplace, VI, or HMC. Recent work (Daxberger et al., "Laplace Redux") has revived the Laplace approximation as a cheap way to add uncertainty to pre-trained models.
• MC Dropout and deep ensembles. Practical approximations to Bayesian uncertainty in deep networks. MC Dropout reinterprets standard dropout as a variational approximation; deep ensembles approximate the posterior by training multiple independent networks. Both work surprisingly well in practice.
• Thompson sampling. In multi-armed bandits and reinforcement learning, sampling actions from the posterior over reward parameters gives near-optimal exploration. Thompson sampling is the default contextual-bandit algorithm at many large tech companies, for Bayesian reasons that were figured out in 1933.
• Bayesian optimisation. Optimising expensive black-box functions by maintaining a Gaussian-process posterior over the objective and choosing the next evaluation point to maximise an acquisition function (expected improvement, UCB, Thompson sampling). The dominant framework for hyperparameter tuning at scale — Google's Vizier, Optuna's TPE, and most AutoML systems are Bayesian at heart.
• Active learning. Pick the next training point to maximise expected information gain about the parameters, measured by posterior entropy reduction or mutual information. The same mathematics as experimental design, now a standard ingredient in data-efficient ML.
• Probabilistic programming languages. Stan, PyMC, Numpyro, Turing.jl, Pyro. Each treats a model as a program that declares priors and a likelihood, and provides MCMC or VI as backends. The whole Bayesian workflow — priors, posterior, checks, predictions — runs inside a single coherent syntax.
• Large-language-model calibration. Interpreting LLM token probabilities as Bayesian beliefs, calibrating them against ground truth, and using predictive uncertainty to decide when to defer to a human or another tool. Bayesian intuition is reappearing in production LLM stacks precisely because calibrated uncertainty matters more at scale.
• Hierarchical models in recommender systems and ads. Partial pooling across users, products, and contexts — a direct lift of Gelman-Hill hierarchical modelling into industrial A/B testing and CTR prediction pipelines.
• Diffusion models. Train a denoising neural network; at inference, run reverse-time SDE dynamics that can be read as sampling from a learned posterior over data given noise. The connection between score matching, denoising diffusion, and Bayesian inference is a lively current research direction.

A conspicuous pattern runs through this list: Bayesian ideas keep getting re-discovered under new names because they work. If you take away one thing from this chapter, let it be the habit of asking, for any ML method you encounter, "what is the implicit prior, and what is the implicit likelihood?" The answers are usually illuminating.