Approximate Inference Methods

Laplace approximation, importance sampling, MCMC, HMC, diagnostics — when to use which

Outside of conjugate models, the posterior p(θ|X) is intractable — you cannot compute it, only approximate it. The question is how much approximation you can tolerate and at what computational cost.

MAP (maximum a posteriori) is the fastest: find the mode and stop, discarding all information about posterior shape. Laplace adds one matrix inversion to recover a Gaussian approximation around the MAP — fast but wrong if the posterior is multimodal or heavy-tailed. MCMC is the gold standard: given enough time, it converges to the true posterior, but "enough time" is often hours or days for complex models. Most production ML systems use MAP with frequentist standard errors and reserve MCMC for settings where exact uncertainty is the product — clinical decision support, scientific inference, hierarchical models. Knowing when each method is appropriate, and critically how to diagnose whether MCMC has actually converged, is what separates theoretical understanding from practical competence.

Key points

• Laplace approximation: find the MAP, compute the Hessian H = -∇² log p(θ|X) at that point, approximate the posterior as N(θ_MAP, H⁻¹). No sampling — just one optimisation and one Hessian computation. When the posterior is genuinely unimodal and approximately Gaussian (which the Bernstein-von Mises theorem guarantees asymptotically), this is excellent. It fails badly for multimodal posteriors because the Hessian only captures local curvature at one mode — if the posterior has mass elsewhere, the Laplace approximation misses it entirely.
• Importance sampling estimates E_p[f(θ)] using a proposal q: E_p[f(θ)] ≈ Σᵢ wᵢ f(θᵢ) where wᵢ ∝ p(θᵢ|X)/q(θᵢ). The effective sample size ESS ≈ (Σwᵢ)²/Σwᵢ² tells you how many i.i.d. samples the weighted set is worth. In high dimensions, IS collapses catastrophically: the typical sets of p and q have negligible overlap, almost all weights are near zero, and a few lucky samples dominate the estimate. ESS < 5% of N means the IS estimate is unreliable regardless of sample count.
• Metropolis-Hastings: propose θ* ~ q(θ*|θ_current), accept with probability α = min(1, p(θ*|X)q(θ_current|θ*) / p(θ_current|X)q(θ*|θ_current)). The intractable normaliser p(X) cancels in the ratio p(θ*|X)/p(θ_current|X) = p(X|θ*)p(θ*) / p(X|θ)p(θ) — this is the key insight that makes MCMC work at all for unnormalised posteriors. You never need to compute p(X); you only need ratios.
• HMC augments the state with momentum and uses gradient information to make large, correlated proposals that are accepted at high rates. Random-walk MH explores via small random steps — inefficient in high dimensions because it takes many steps to traverse the posterior. HMC uses the gradient of log p(θ|X) to simulate Hamiltonian dynamics, enabling large steps that respect the posterior geometry. NUTS (No-U-Turn Sampler) adapts step size and trajectory length automatically and is the default in Stan and PyMC.
• MCMC diagnostics are non-negotiable. R-hat (Gelman-Rubin): run K independent chains from different starting points. R-hat = √(total variance / within-chain variance). R-hat < 1.01 → chains have converged to the same distribution. R-hat > 1.1 → chains are exploring different regions; you do not yet have samples from the posterior. ESS accounts for within-chain autocorrelation: ESS < 100 per parameter means high Monte Carlo error. Trace plots should look like a "hairy caterpillar" — no trends, no sticking, all chains overlapping.
• HMC divergences are a hard stop, not a warning. A divergence means the leapfrog integrator encountered extreme posterior curvature and went numerically unstable. Any divergences mean the posterior geometry is pathological and posterior estimates are biased. Non-centred reparameterisation fixes the most common cause. Never report results from a sampler with divergences.
• Non-centred parameterisation is the most important practical HMC fix for hierarchical models. Centred: μ_i ~ N(μ, σ²) directly — when σ is small, the posterior forms a funnel: narrow at the tip (small μ_i variations) and wide at the top. HMC's step size must be tiny to navigate the narrow funnel tip, causing slow mixing everywhere else. Non-centred: write μ_i = μ + σ·z_i where z_i ~ N(0,1), sample z_i instead. The funnel geometry disappears. If you see divergences in a hierarchical model, this is the first thing to try.
• Method selection guide: MAP + frequentist CIs for production systems that need speed and scale. Laplace for post-MAP uncertainty estimates where full MCMC is too slow (last-layer BNNs). NUTS/HMC in Stan or PyMC for serious scientific Bayesian analysis with < ~1M parameters. VI (CAVI, BBVI) for large-scale latent variable models where sampling is too slow. Deep ensembles + conformal prediction for large neural networks where full Bayes is infeasible. The choice is about what you need from uncertainty: a fast approximation or an exact distribution.

Recap

• Cost/accuracy ladder: MAP (mode only) → Laplace (Gaussian at MAP) → MCMC (exact, but hours-to-days).
• Laplace = $N(θ_{MAP}, H^{-1})$ — one Hessian; fails for multimodal/heavy-tailed posteriors (only local curvature).
• Importance sampling collapses in high-D: typical sets of p and q barely overlap; ESS < 5% of N → unreliable.
• Metropolis-Hastings ratio cancels $p(X)$ — the key that makes MCMC work on unnormalised posteriors (only ratios needed).
• HMC/NUTS use gradients for large, high-acceptance proposals; NUTS auto-tunes and is the Stan/PyMC default.
• R-hat > 1.01 = not converged — the samples are not from the posterior. HMC divergences are a hard stop, not a warning.
• Non-centred reparameterisation ($μ_i = μ + σ·z_i$) removes the funnel geometry that causes hierarchical-model divergences.

Check your understanding