Computational Bayesian Methods | STAT-101

Where the Framework Runs Out of Algebra

The previous two lessons built the Bayesian machine and showed it running by hand. Every posterior is proportional to the prior times the likelihood, and for a conjugate pair the posterior lands in a known family, so updating is closed-form arithmetic. The catch is that conjugate pairs are the exception. Write the posterior out in full and the denominator is an integral over the whole parameter space:

The Posterior, In Full

p(\theta \mid D) = \dfrac{p(D \mid \theta)\,p(\theta)}{\int p(D \mid \theta)\,p(\theta)\,d\theta}

The numerator is easy: for any candidate θ we can evaluate the prior and the likelihood. The denominator — the marginal likelihood P(D) — is that integral, and outside conjugate models it has no closed form.

For one or two parameters you could grind the integral out numerically on a grid. But real models have dozens, hundreds or thousands of parameters, and a grid over that many dimensions contains more cells than there are atoms worth of computer time. This is the wall that made Bayesian inference impractical for most of the twentieth century, and it is exactly the wall that Markov Chain Monte Carlo (MCMC) walks around.

The MCMC Idea

Instead of evaluating the posterior, MCMC samples from it. The name is two ideas glued together. Monte Carlo: if you can draw samples from a distribution, you can estimate any feature of it — its mean, a 95% interval, the probability of a tail — just by computing that feature on the samples and averaging. Markov chain: a sequence of states in which each one depends only on the state just before it.

The trick is to build a Markov chain whose stationary distribution is the posterior itself. Start it anywhere, let it run, and after a while the states it wanders through are draws from P(θ∣D) — not independent draws, but draws all the same. The decisive advantage is that the chain only ever needs to compare the posterior at two points, so it works with ratios of the posterior. In a ratio the intractable denominator cancels: it is the same constant top and bottom. MCMC never has to compute the very thing that was impossible to compute.

Metropolis–Hastings

The oldest and most general recipe for such a chain is Metropolis–Hastings. From the current state x, it proposes a candidate x′ from some proposal distribution q, then decides whether to move there or stay, using an acceptance probability:

Acceptance Probability

Draw a uniform number on (0, 1); if it is below α, move to x′, otherwise stay at x and record x again. Note p appears only as a ratio, so it need only be known up to a constant — which is the whole point.

When the proposal is symmetric — as likely to suggest a step from x to x′ as the reverse, which a Gaussian centred on the current point always is — the q ratio equals one and the rule collapses to the original Metropolis form:

Symmetric Case

Any uphill move (a more probable x′) is always accepted; a downhill move is accepted with a probability equal to how much less probable it is. That controlled willingness to go downhill is what lets the chain explore instead of just climbing to the mode.

Why does this land on the posterior? The acceptance rule is designed to satisfy detailed balance: the chain is as likely to be found stepping from x to x′ as from x′ back to x, when x is drawn from the posterior. A chain that obeys detailed balance with respect to a distribution has that distribution as its stationary distribution. Run it long enough and it forgets where it started and settles into sampling the posterior.

Tuning: Acceptance versus Mixing

Metropolis–Hastings has one knob that decides whether it works well or barely works at all: the width of the proposal, the typical size of the step it suggests. It sets up a direct trade-off between how often moves are accepted and how well the chain mixes — how quickly it moves across the posterior.

Too small — almost every tiny step is accepted, but the chain shuffles in place. Successive draws are nearly identical, so they are heavily autocorrelated and carry little independent information.
Too large — almost every step lands somewhere far less probable and is rejected, so the chain sits still for long stretches, recording the same value over and over.
A healthy band in between — steps big enough to travel, accepted often enough to keep moving. This is where the sampler actually explores.

Rule of Thumb

For a high-dimensional random-walk Metropolis sampler there is a famous target: an acceptance rate near 0.234 is close to optimal (Roberts, Gelman & Gilks, 1997). Treat it as a diagnostic, not a law — if you are accepting 95% of proposals your steps are too timid, and if you are accepting 2% they are too bold. The right answer for a one-dimensional problem is higher, nearer 0.44.

Gibbs Sampling

Metropolis–Hastings is general but wasteful: it throws away every rejected proposal. When the model has more structure, Gibbs sampling exploits it. The idea is to update one parameter at a time by drawing it from its full conditional distribution — its distribution given the current values of all the other parameters. If every full conditional is a distribution you can sample from directly, you never reject anything.

The clean example is the correlated bivariate normal. If (x₁, x₂) is bivariate normal with correlation ρ and unit variances, then each coordinate's full conditional is itself an ordinary normal:

A Full Conditional

x_1 \mid x_2 \sim \mathcal{N}(\rho x_2,\; 1 - \rho^2)

Hold x₂ fixed and x₁ is normal, centred on ρx₂ with variance 1−ρ²; by symmetry the same holds with the roles swapped. Sample x₁ from this, then x₂ from its own conditional, and repeat — every draw is accepted.

Gibbs sampling mixes well and needs no tuning when you can use it, which is why it powered the first generation of practical Bayesian software. Its limitation is exactly its requirement: you can only use it when the full conditionals are known distributions you can draw from. When they are not, you fall back to Metropolis–Hastings, or to a hybrid that uses Gibbs where it can and Metropolis steps elsewhere.

Diagnostics: Did the Chain Converge?

MCMC gives you samples whether or not the chain has actually reached the posterior, so its output cannot be trusted blind. A handful of diagnostics tell you whether to believe it:

Trace plots. Plot each parameter against iteration number. A converged, well-mixing chain looks like a fuzzy horizontal band — a “hairy caterpillar.” Long flat runs, slow drift, or a chain that never settles are all warning signs.
Burn-in. The chain starts wherever you put it, which may be far from the bulk of the posterior. Discard the early iterations — the burn-in — so the estimates are not dragged toward the arbitrary starting point.
Autocorrelation and effective sample size. Consecutive draws are correlated, so 10,000 MCMC samples are worth fewer than 10,000 independent ones. The effective sample size (ESS) reports how many independent draws they are equivalent to; a small ESS means high autocorrelation and shaky estimates.
R-hat across chains. Run several chains from dispersed starting points. If they have converged they should all be exploring the same region, so the variance between chains should match the variance within them. The Gelman–Rubin statistic R̂ measures this, and a value near 1 is the signal that the chains agree.

You Do Not Write These By Hand

In practice nobody codes a sampler from scratch. Probabilistic programming languages — Stan, PyMC and others — let you state the model (priors and likelihood) in a few lines and compile it to an efficient sampler automatically, usually a modern gradient-based method rather than plain random-walk Metropolis. The concepts in this lesson are what those tools are doing under the hood, and they are what the diagnostics above still ask you to check.

Try It

The companion interactive MCMC sampler runs a real Metropolis–Hastings chain and a real Gibbs sampler in the browser. Drag the proposal width and watch acceptance trade off against mixing in the trace plot — small steps crawl, large steps stick, and the band in between explores. It is the fastest way to feel what “mixing” means.

Key Takeaways

The posterior's normalising constant is an integral that is intractable outside conjugate models, so the posterior usually cannot be evaluated directly.
MCMC samples the posterior instead of evaluating it: build a Markov chain whose stationary distribution is the posterior, and estimate any quantity by averaging over its draws. Only posterior ratios are needed, so the missing constant cancels.
Metropolis–Hastings proposes a move and accepts it with probability α = min(1, …); for a symmetric proposal this is min(1, p(x′)/p(x)). Detailed balance is what makes the posterior stationary.
The proposal width trades acceptance against mixing: too small crawls, too large sticks. A high-dimensional acceptance rate near 0.234 is a useful rule of thumb.
Gibbs sampling draws each parameter from its exact full conditional and accepts every draw — excellent when the conditionals are known, unavailable when they are not.
Always diagnose: trace plots, burn-in, autocorrelation and ESS, and R̂ across multiple chains. Probabilistic programming tools (Stan, PyMC) build the sampler for you but do not excuse the checks.