MCMC
You have a posterior and it is written down completely — an unnormalised density you can evaluate at any point you name. What you cannot do is integrate it. The normalising constant is a high-dimensional integral with no closed form, and every summary you care about — a mean, a credible interval, a marginal, an expectation of some function — is another integral over the same intractable object. So you stop trying to compute and start trying to visit. If you can produce a stream of points whose long-run frequency matches the density, then averages over the stream estimate every expectation you wanted. Markov chain Monte Carlo is the machinery that turns "I can evaluate it" into "I can sample from it", and Metropolis–Hastings is the one idea underneath almost all of it.
The question
When the normaliser is the whole problem
Bayes' rule gives the posterior up to a constant. Write the data as $\mathcal{D}$ and the parameter as $\theta$: the posterior is proportional to the likelihood times the prior, and the missing denominator is the evidence $p(\mathcal{D})=\int p(\mathcal{D}\mid\theta)p(\theta)\,d\theta$. In one dimension that integral is a quadrature problem. In fifty dimensions it is hopeless, because the number of evaluations needed grows exponentially with dimension. Yet the numerator is easy: hand me any $\theta$ and I can tell you $p(\mathcal{D}\mid\theta)p(\theta)$ in microseconds.
The previous part made this painless for a coin, where a Beta prior met a binomial likelihood and the posterior landed back in the Beta family. That is the exception. As soon as the model has a couple of layers — a hierarchical prior, a mixture, a neural likelihood — no conjugate family survives and the normaliser becomes an integral you cannot write down. The conjugate case taught us what a posterior is; this part is about what to do when you cannot compute one.
Here is the shift in viewpoint. I do not need the normalising constant if I can sample from the distribution. A cloud of points drawn from the posterior answers the same questions the density would. The mean is the average of the cloud; the 95% credible interval is its middle 95%; the marginal of one coordinate is its histogram along that axis; and an expectation $\mathbb{E}[f(\theta)]$ is the plain average of $f$ over the points. This is exactly the Monte Carlo principle of Part 8, applied to a distribution rather than a fixed region.
Ordinary sampling needs to know how to draw — inverse-CDF, rejection with an envelope, importance weights. All three need something more than pointwise evaluation: a quantile, an upper bound, a proposal you can sample. The posterior offers none of that. What it does offer is a relative weight between any two points, and that turns out to be enough, provided you are willing to give up independence. Markov chain Monte Carlo draws a correlated sequence whose stationary distribution is the target. The chain has memory, the samples are not independent, and every diagnostic in the subject exists to check that the memory is not lying to you.
Metropolis–Hastings
Propose a move; accept it with the right probability
The algorithm is a walk with a rule. You are standing at a point $x$ and you want to know the target density there, but only up to a constant, so call the unnormalised value $\pi(x)$. Pick a proposal distribution $q(x'\mid x)$ that you can sample from — a Gaussian centred on where you are is the standard choice — and draw a candidate $x'$. Then decide whether to move. The decision is not "move if the density is higher"; that would march uphill and pile every sample on the mode. It is a coin flip whose bias is the ratio of the densities, capped at one.
Look at what the formula asks for, and notice what it does not ask for. Both $\pi$ values appear as a ratio, so any constant that multiplies the target cancels — the intractable evidence integral vanishes. That single observation is why MCMC works at all. If the proposal is symmetric, $q(x\mid x')=q(x'\mid x)$, the ratio simplifies to $\pi(x')/\pi(x)$ and the algorithm is the plain Metropolis rule. When the target is flat at the mode and you propose a worse point, you still accept sometimes: a downhill step is how the chain escapes and how the histogram fills in the tails.
Why does this produce the right distribution? The acceptance rule is constructed to satisfy detailed balance, $\pi(x)\,T(x\to x')=\pi(x')\,T(x'\to x)$, where $T$ is the full transition probability of proposing and accepting. Detailed balance makes $\pi$ stationary, and if the chain is irreducible — it can reach any region with positive density — the stationary distribution is unique. Then a theorem from the Markov-chains part says the empirical distribution of the visited points converges to $\pi$. That is the whole justification, and it is short because it never touches the normaliser either.
The demo below runs one chain on a correlated two-dimensional Gaussian, so the target is a leaning ellipse you can see. Start where the chain starts, then move it a step at a time, or run a batch. The pink path is the chain, the dots are everywhere it has been, and the blue contours are the truth it is trying to fill in.
Blue ellipses are the one-, two- and three-sigma contours of the target. The pink path is the Metropolis chain; each dot is one accepted state. The bold dot is where the chain currently sits.
The same walk can be lifted onto the density itself. The next panel draws the target as a surface — height is probability density — and drops the chain onto it. Every sample sits under the surface at its own density, so a chain that has found the mass shows up as a trail of markers on the flanks of the hill, and a chain that is still climbing shows up as a lonely thread in the low-density corner. Change the step size and watch both the scatter and the surface trail respond together.
Plotly surface of the target density with the chain overlaid as a scatter of markers, positioned at their own density height.
Notice that nothing here required an integral. To run any of it we needed exactly one function: given a point, return the unnormalised log density. That is the entire interface MCMC asks of a model, and it is why the technique reaches into places where closed-form posteriors do not exist — a hierarchical model with hundreds of parameters, a latent-variable model whose likelihood is itself an integral, a Bayesian neural network whose prior is over millions of weights.
Step size, acceptance and mixing
The one knob that decides whether you learn anything
The proposal width h is the algorithm's only real tuning parameter, and it fails in opposite directions. Make h tiny and almost every proposal is accepted, because the density barely changes over a short hop — but the chain creeps. It explores the posterior at a snail's pace, neighbouring samples are nearly identical, and the effective number of independent samples is a fraction of the count. Make h huge and the chain proposes points far out in the tails, whose density ratio is minuscule, so almost everything is rejected and the chain sticks where it is for long stretches. The samples are then few in number and far apart in the way that matters least.
Between the extremes sits a sweet spot where the chain moves fast enough to traverse the distribution in a reasonable number of steps and still accepts often enough to make progress. The acceptance rate is the readout that tells you which regime you are in. A rate near one means sticky and slow; a rate near zero means jumpy and stuck; a rate somewhere in the middle means mixing. Drag the slider in the demo above from left to right and watch the trace flatten, then explode, then settle.
Those numbers are not folklore. In the late 1990s researchers worked out the scaling limit of random-walk Metropolis as dimension grows and found the proposal that optimises the diffusion speed of the chain accepts about 23.4% of the time. In one or two dimensions the optimum drifts up toward roughly one half. The $d^{-1/4}$ rule for the step length is the same result read the other way: as the space grows you must shrink the jump, but only as the fourth root, so the algorithm degrades gracefully rather than catastrophically. The demo here is two-dimensional, so aim for a rate near one half.
The consequence of a bad step size is autocorrelation. Adjacent samples in the stream are correlated because the chain has memory, and the stronger the memory the more samples it takes to convey one fresh piece of information. The trace plot is the visual diagnostic: a healthy chain looks like a fuzzy caterpillar with no slow drifts and no long flat stretches; a sticky chain looks like a stepped line that barely moves; a jumpy chain looks like a spiky mess with the same value repeated. The histogram below is the payoff — only when the trace is healthy does the chain's own distribution line up with the target curve.
Gibbs and Hamiltonian moves
Two ways to move better than a random walk
Metropolis–Hastings is completely general, and generality is expensive: it wanders. A structured model usually hands you something a blind proposal cannot see — the conditional distribution of each coordinate given the others. Gibbs sampling takes that gift and updates one coordinate at a time, drawing each from its exact conditional. There is no accept/reject step at all. Every proposal is a draw from the correct conditional, so every move is accepted, and the chain updates the full vector once per sweep.
The gains come with a catch. Because each move is axis-aligned, Gibbs struggles when coordinates are strongly correlated: it takes a staircase of tiny horizontal and vertical steps to travel along a diagonal ridge, and the chain looks like it is shuffling sideways while barely advancing. The demo below makes the contrast literal. The same correlated Gaussian is walked by Metropolis, which can move diagonally but rejects moves, and by Gibbs, which never rejects but can only step along the axes.
For the Gaussian the conditional is exact and cheap, which is why Gibbs is the natural sampler for the multivariate Gaussian and for every model built out of it — the block updates inside latent-variable models, the sweep in an Ising model, the alternating draws that make a Bayesian mixture tractable. When the conditional is not available in closed form you can still use a Metropolis step inside a Gibbs sweep, which is the hybrid that most real samplers are.
Grey: the Metropolis chain. Pink: the Gibbs chain, drawn through its two half-steps per sweep so the axis-aligned staircase is visible.
Hamiltonian Monte Carlo attacks the random walk itself. Instead of proposing a blind jump, give the chain a simulated physical momentum and let it glide along the gradient of the log density for several small leapfrog steps. Positions and momenta update in an alternating pattern that conserves energy almost exactly, so a proposal can travel a long way along a level set of the density while still being a plausible draw. The gradient buys large, well-directed moves; the length of the trajectory is the product of the step size $\epsilon$ and the number of leapfrog steps $L$.
The three lines are the leapfrog integrator, applied $L$ times and followed by a Metropolis accept/reject on the total energy. In the correlated Gaussian below, the gradient points along the lean of the ellipse, so HMC shoots across the ridge in a handful of moves where Metropolis has to blunder into it. The price is that HMC needs a gradient — cheap for a Gaussian, available by backpropagation for a neural likelihood, and the reason Hamiltonian methods swept through Bayesian statistics once autodiff matured. This is the same gradient calculus that drives nonlinear optimization, repurposed from finding a minimum to exploring a distribution.
Grey: random-walk Metropolis. Pink: Hamiltonian trajectories. The HMC chain crosses the correlated ridge in long, directed arcs instead of blind hops.
Trace plots and burn-in
Reading a chain honestly
A Markov chain started at an arbitrary point has not yet forgotten where it began. The early samples are a transient, drifting from the starting value toward the region of high probability, and during that drift they are not draws from the target at all. Discarding that warm-up is called burn-in, and the trace plot is how you decide how much to throw away. Start the chain far out in a low-density corner and you will see the first hundred iterations of the trace climb steadily; only after it starts wiggling around a stable level is the chain sampling.
Burn-in is necessary but not sufficient, and it is easy to overstate its importance. Throwing away the first $B$ samples cannot create information; it removes a bias at the cost of a few samples. The bigger issue is autocorrelation, which persists after warm-up. Samples that sit near each other in time are correlated, so the $N$ samples you kept may carry only $N_{\text{eff}}$ samples' worth of independent information. The leftmost panel below starts the chain in a corner; move the burn-in slider and watch the histogram fill in. The autocorrelation panel shows the memory directly.
Read the trace for three diseases. A slow drift with no return means the chain is still converging — run longer and burn more, or the estimate is biased. Long flat stretches mean the chain is stuck, a proposal or a multimodality problem. A jagged line that hugs a narrow band means the chain is mixing well and the effective sample count is close to the raw count. None of these is visible in a single summary number; the trace is the honest view, which is why every serious sampler prints one.
The autocorrelation function should decay to near zero within a few lags for a well-tuned chain. If it lingers — heavy tails out to tens or hundreds of lags — you are not getting the sample size you paid for. The cure is a better proposal, a reparameterisation that breaks the correlation between coordinates, or a different algorithm entirely. The diagnostic is the same in every case: plot the trace, plot the autocorrelation, and reduce the effective sample size accordingly before you quote an error bar.
A chain started in a low-density corner. The shaded region is discarded burn-in B; the histogram on the right uses only the kept samples. The bottom panel is the autocorrelation of the kept x sequence.
Where this shows up
When the posterior has no formula
The posterior you cannot integrate
The conjugate part left the posterior as a closed-form Beta; the moment the model has a hierarchy or a non-conjugate likelihood, that closed form is gone and MCMC is the default way to get samples from it anyway.
Conditionals and Gaussians
Gibbs sampling lives off the closed-form conditionals of the multivariate Gaussian: block updates, marginal slices and the precision matrix all become moves on the chain rather than integrals.
Gradients for better proposals
Hamiltonian moves need $\nabla\log\pi$, the same gradient machinery that powers nonlinear optimization. Where a map or a pose graph has an analytic gradient, HMC turns it into directed exploration instead of a blind walk.
Uncertainty over model parameters
When a model's weights are given a posterior rather than a point estimate, that posterior is sampled, not solved. The scoring of such distributions is what model evaluation measures, and the same diagnostics decide whether to believe the uncertainty it reports.
Cheat sheet
Every formula in one place
| Idea | Formula | Reading |
|---|---|---|
| Target | $\pi(\theta)\propto p(\mathcal{D}\mid\theta)\,p(\theta)$ | Known up to a constant; that constant is never computed. |
| MH acceptance | $\alpha=\min\!\left(1,\dfrac{\pi(x')\,q(x\mid x')}{\pi(x)\,q(x'\mid x)}\right)$ | Accept the downhill move sometimes; the ratio cancels the evidence. |
| Symmetric proposal | $\alpha=\min\!\left(1,\pi(x')/\pi(x)\right)$ | Plain Metropolis: only the density ratio matters. |
| Detailed balance | $\pi(x)\,T(x\to x')=\pi(x')\,T(x'\to x)$ | Makes $\pi$ stationary; with irreducibility, unique. |
| Gibbs update | $x_i^{(t+1)}\sim p\!\left(x_i\mid x_{-i}^{(t)}\right)$ | Exact conditionals; every move accepted. |
| Gaussian conditional | $x_1\mid x_2\sim\mathcal{N}\!\left(\mu_1+\tfrac{\Sigma_{12}}{\Sigma_{22}}(x_2-\mu_2),\;\Sigma_{11}-\tfrac{\Sigma_{12}^2}{\Sigma_{22}}\right)$ | The workhorse Gibbs step for correlated clouds. |
| HMC leapfrog | $p\!+\!=\!\tfrac{\epsilon}{2}\nabla\log\pi(x);\;x\!+\!=\!\epsilon p;\;p\!+\!=\!\tfrac{\epsilon}{2}\nabla\log\pi(x)$ | Gradient-guided moves; trajectory length $\epsilon L$. |
| Step-size rule | $\mathbb{E}[\alpha]\approx0.234$ (high $d$), $\approx0.5$ (1–2D), $h\propto d^{-1/4}$ | Tune acceptance to the middle; too high is sticky, too low is stuck. |
| Burn-in | discard first $B$ samples | Remove the transient from the starting point; watch the trace. |
| Autocorrelation | $\rho_k$; $N_{\text{eff}}\approx N/(1+2\sum\rho_k)$ | Decay fast or you have fewer samples than you counted. |
Further reading
Where to go deeper
- Metropolis, Rosenbluth, Rosenbluth, Teller and Teller, "Equation of state calculations by fast computing machines", Journal of Chemical Physics, 1953 — the original algorithm, from statistical physics.
- W. K. Hastings, "Monte Carlo sampling methods using Markov chains and their applications", Biometrika, 1970 — the general, asymmetric-proposal form and the detailed-balance argument.
- Radford Neal, "MCMC using Hamiltonian dynamics", 2011 — the definitive chapter on leapfrog, energy and tuning.
- Christian Robert and George Casella, Monte Carlo Statistical Methods, 2004 — the textbook treatment of Metropolis–Hastings and Gibbs, with the convergence theorems.
- Andrew Gelman et al., Bayesian Data Analysis, chapter 11 — convergence diagnostics, trace plots, $\hat R$ and effective sample size in practice.
- David MacKay, Information Theory, Inference, and Learning Algorithms, chapter 29 — the same ideas written for exactly this audience.