When conjugacy fails: MCMC

What do you do when the posterior has no formula?

Bayesian inference needs one object: the posterior $p(\theta \mid y)$, the distribution of the parameters after seeing the data. Everything you might report — a mean, a credible interval, a predictive probability — is an integral against that distribution. When the posterior is one of the handful of conjugate families, the integral is arithmetic and the story ends. When it is not, the shape exists but refuses to be written down, and the integral is what stands between you and every answer you wanted.

The escape is to stop trying to evaluate the integral and start trying to draw from the distribution. If you have samples $\theta_1, \dots, \theta_N$ that behave as though they came from the posterior, then any expectation you care about is just an average of numbers you already hold. You never differentiated or integrated anything. This part builds the machine that produces those samples, and then makes its failure modes visible.

The normalising constant is the wall

A conjugate prior is a prior chosen so that the posterior lands back in the same family as the prior. The choice is convenient but narrow. It works for the Beta–Binomial pair, for the Gaussian–Gaussian pair, for Gamma–Poisson, and for a short list of others. Each pairing has been tabulated precisely because it is rare, and the previous two parts spent their time inside that list.

Step outside it and the algebra collapses. Give a logistic regression a Gaussian prior on its coefficients and the posterior is the product of a smooth prior and a likelihood with $n$ logistic terms. It is a perfectly well-defined density, but it is not a Gaussian, not a Student-$t$, not anything with a name. Put a prior on the parameters of a distribution whose mean is a mixture of two components and the posterior becomes a multi-modal surface that no one-line formula describes. Add a hierarchy — parameters drawn from a population distribution whose own parameters are unknown — and the posterior acquires an integral inside it.

The obstruction is always the same piece of the definition,

The numerator is easy: it is a product you can evaluate at any $\theta$ with a line of code. The denominator is an integral over the whole parameter space, and without a conjugate trick it has no closed form. It is the same number for every $\theta$, so it changes nothing about the shape of the posterior — only its scale. That observation is the crack in the wall, and Metropolis–Hastings drives straight through it.

Monte Carlo turns an integral into an average

Suppose, for a moment, that samples from the posterior were handed to you. The posterior mean is then just the average of the samples, and any other posterior expectation follows the same pattern:

This is the Monte Carlo idea in one line. A credible interval becomes a pair of sample quantiles. A predictive probability becomes a fraction of samples that pass a test. Nothing requires the posterior to have a name, and the error of the approximation shrinks like $1/\sqrt{N}$ — crucially, that rate does not get worse as the number of parameters grows. A grid over ten dimensions with twenty points per axis is already $20^{10}$ cells, while a million samples remain a million samples in any dimension. The price is that the samples are random and correlated, and the $1/\sqrt{N}$ describes an idealised independent draw.

So the whole problem reduces to one question: how do you draw from a distribution you can only evaluate up to an unknown constant? Rejection sampling needs a bounding envelope that is hopeless in high dimensions, and importance sampling collapses when the proposal and the target disagree. The method that survives is to build a Markov chain — a sequence of values where each depends on the last — and arrange for the chain to spend time in each region in proportion to the target density. If a chain does that, its values are a sample from the target for every purpose that matters. Where optimizers hunt for the single highest point of a function, a sampler wanders over the whole surface and reports how much of it is high.

Propose, then accept or reject

Metropolis–Hastings is a recipe for that chain, and it is startlingly short. Standing at the current value $\theta$, propose a nearby value $\theta'$ from a proposal distribution $q(\theta' \mid \theta)$. Evaluate the target at both points. Accept the proposal with probability

and if the proposal is rejected, stay where you are and record the current value again. The most common choice is a symmetric random walk: $q(\theta' \mid \theta) = \mathcal{N}(\theta, s^2 I)$, so the two proposal densities cancel and the acceptance probability is just the ratio of target densities, $\min(1, p(\theta' \mid y)/p(\theta \mid y))$. A move uphill is always accepted; a move downhill is accepted sometimes, with a probability equal to how far downhill it is.

The denominator of Bayes' rule appears in both the numerator and denominator of that ratio, so it cancels. This is the point of the whole construction: the chain needs the posterior only up to a constant, which is precisely the part you can compute. The resulting transition rule satisfies a condition called detailed balance, and detailed balance is what forces the target distribution to be the chain's stationary distribution — the distribution it converges to and then samples from forever.

The one free parameter is the proposal width $s$. It behaves like the step size of an explorer who cannot see the map. Too small and almost every proposal is accepted, but each step is a crawl: the chain mixes slowly and consecutive samples are nearly identical. Too large and proposals jump far into low-density territory, where they are almost always rejected, so the chain sits still. The useful regime is between, where the step size is comparable to the width of the target. For a random walk in more than a couple of dimensions the classical target is an acceptance rate near a quarter, and the demo's readout will show you what each extreme looks like.

Proposal width, trace plot, acceptance rate, burn-in, R-hat

The target below is a correlated two-component mixture in the plane: two elliptical clumps, tilted in opposite directions, with a moderate dip between them. No normalising constant is used anywhere — the chain only ever evaluates the log of the unnormalised density. Four chains start from the corners of the plot, which is deliberately unhelpful, because watching them find the modes is part of the lesson. The shaded background is the target; the dots are the samples each chain has collected; the solid marker is the most recent position of the first chain.

Press Step a few times and watch a single proposal at a time: the marker hops uphill without hesitation and sometimes declines a downhill move. Press Run and the four chains scatter out of their corners and converge onto the two clumps. The trace plot is the same story in one dimension: a chain that has not yet converged shows a slow drift from its starting value, and a chain that has converged shows a fuzzy band of jitter with no trend. The acceptance-rate readout names the proposal width for you, and the Burn-in slider sets how many early samples are greyed out and excluded from the estimates.

How to tell whether the samples are worth averaging

Burn-in. A chain started at an arbitrary point is not yet sampling from the target; it is travelling towards it. Those early iterations carry a bias that has nothing to do with the posterior, and the standard remedy is to discard them. The grey region of the trace plot is exactly that discarded prefix. Burn-in is a crude instrument — you cannot prove a chain has arrived, only fail to see it still moving — but on a trace plot that has flattened into a stationary fuzz the transient is usually obvious.

Autocorrelation and effective sample size. Successive values of a random-walk chain are correlated, so $N$ samples are worth less than $N$ independent ones. The autocorrelation $\rho_k$ measures how much a value resembles the value $k$ steps earlier, and the effective sample size converts the raw count into the equivalent number of independent draws,

When the proposal width is too small the autocorrelation decays slowly, the denominator is large, and the ESS is a small fraction of $N$. This is why a high acceptance rate is not a virtue in itself: accepting every proposal is worthless if the accepted moves are microscopic. The cure is the same knob that the acceptance rate is reporting on.

R-hat. The trace plot's weakness is that a single chain can look stationary while being stuck in one mode and never visiting another. The fix is to run several chains from over-dispersed starting points and compare them. R-hat pools two variance estimates: the average variance within each chain and the variance between the chains' means. If the chains have all converged to the same distribution, the two agree and

If they are still exploring different regions, the between-chain spread inflates R-hat above one. The conventional rule of thumb is that $\hat{R} < 1.1$ is passable and values near $1.01$ are better. In the demo, R-hat is computed on the retained portion of all four chains; watch it fall towards one as the chains find both clumps, and watch it stall above one if you leave the starting points stranded. This diagnostic is the direct descendant of the Markov-chain filters used in robot state estimation, where a cloud of hypotheses is propagated and compared in exactly this spirit.

Sampling when the formula runs out

The pattern generalises far beyond these two. Every hierarchical model, every latent-variable model and every posterior over neural-network weights is a distribution you sample because you cannot integrate. What changes is the proposal mechanism: Hamiltonian Monte Carlo uses gradients to propose long, informed moves, and Gibbs sampling proposes along coordinate slices. What stays constant is the diagnostic vocabulary of this part — acceptance, trace, autocorrelation, effective sample size and R-hat — and the Monte Carlo reduction of integrals to averages that made the whole enterprise worthwhile.

Further reading

The references below treat MCMC as a practical instrument with failure modes rather than a theorem to admire. If you take away one thing, take away the picture of a chain wandering over a surface it can only evaluate up to a constant, and of diagnostics as the instruments that tell you whether the wandering was honest.

Gelman and colleagues is the standard applied reference, with the convergence chapter that made R-hat a household name. MacKay derives Metropolis–Hastings from detailed balance in a handful of pages. Andrieu, de Freitas, Doucet and Jordan is the machine-learning survey that connects the algorithm to the models that need it, and Betancourt's conceptual introduction is the clearest modern account of why a random walk is the wrong default once gradients are available.

• Andrew Gelman, John Carlin, Hal Stern, David Dunson, Aki Vehtari and Donald Rubin, Bayesian Data Analysis, 3rd ed., chapter 11 — MCMC, convergence and effective sample size.
• David MacKay, Information Theory, Inference, and Learning Algorithms, chapter 29 — Metropolis–Hastings from detailed balance.
• Christophe Andrieu, Nando de Freitas, Arnaud Doucet and Michael Jordan, "An Introduction to MCMC for Machine Learning", Machine Learning, 2003.
• Michael Betancourt, "A Conceptual Introduction to Hamiltonian Monte Carlo", 2017 — why informed proposals beat a random walk.

Cheat sheet