Bayesian Computation & MCMC

The trick underneath every MCMC algorithm: we never need $m(\mathbf{y})$. To sample from $p(\theta \mid \mathbf{y})$ we only need ratios $\pi(\theta')/\pi(\theta)$ (in the MCMC-theory notation introduced in §26.3, where $\pi$ denotes the posterior-as-target), and the unknown $m(\mathbf{y})$ cancels from numerator and denominator. This is the structural reason every MH / HMC / NUTS algorithm moves by proposing a candidate $\theta'$ and comparing its unnormalized density to the current state’s.

Topic 26 develops four algorithms end-to-end: Metropolis–Hastings (§26.2, the featured recipe), Gibbs sampling (§26.3, the conjugate-full-conditional specialization), Hamiltonian Monte Carlo (§26.4, the gradient-aware extension), and the No-U-Turn Sampler (§26.5, the self-tuning HMC default that Stan and PyMC ship). Convergence diagnostics follow in §26.6. The theoretical spine is §26.7’s Markov-chain ergodic theorem, which generalizes Topic 10’s Theorem 10.5 (Kolmogorov/Etemadi SLLN) from iid sequences to dependent Markov chains. Ten specializations and alternatives — adaptive MCMC, reversible-jump, Riemann-manifold HMC, stochastic-gradient MCMC, sequential Monte Carlo, variational inference, parallel tempering, nested sampling, Bayesian nonparametric MCMC, slice sampling — each get a forward-pointing paragraph in §26.10, not a dedicated section. Topic 26’s bar is pedagogical: does this method change how the reader thinks about proposals, acceptance, or convergence?

MCMC is asymptotically exact: under irreducibility, aperiodicity, and positive recurrence, long-run chain averages converge almost surely to posterior expectations (§26.7 Thm 5). The price is computational — chains take time to mix, and diagnosing convergence is an empirical art. Variational inference trades exactness for speed (§26.10 Rem 26); Topic 26 commits to the exactness side of the trade.

Let $\pi$ be a target density on $\mathcal{X} \subseteq \mathbb{R}^d$, known up to a normalizing constant. A proposal kernel is a conditional density $q(x' \mid x)$ that generates a candidate $x'$ given the current state $x$. Given a proposal, the MH acceptance probability is

$\alpha(x, x') \;=\; \min\!\left\{1,\; \frac{\pi(x') \, q(x \mid x')}{\pi(x) \, q(x' \mid x)}\right\}.$

The MH transition kernel is

$K(x, dx') \;=\; \alpha(x, x') \, q(x' \mid x) \, dx' \;+\; r(x) \, \delta_x(dx'), \qquad r(x) \;=\; 1 - \!\int \alpha(x, y) \, q(y \mid x) \, dy,$

where the first term covers accepted moves and the second (the rejection mass $r(x)$) parks the chain at $x$ when the proposal is rejected.

A distribution $\pi$ is invariant for a transition kernel $K$ if $\int \pi(x) K(x, A) \, dx = \pi(A)$ for every measurable $A \subseteq \mathcal{X}$. A sufficient condition for invariance is detailed balance: $\pi(x) k(x, x') = \pi(x') k(x', x)$ for all $x, x'$, where $k$ is the off-diagonal density of $K$.

The MH transition kernel $K$ of Def 1 satisfies detailed balance with respect to $\pi$. Consequently $\pi$ is invariant for $K$.

The canonical 2-D tuning problem. Target $\pi(q) \propto \exp\!\bigl(-\tfrac{1}{2}\bigl[(1-q_1)^2 + 10(q_2 - q_1^2)^2\bigr]\bigr)$ (Rosenbrock banana). Use a symmetric Gaussian proposal $q(x' \mid x) = \mathcal{N}(x, \sigma^2 I)$, so the proposal density ratio in $\alpha$ drops out: $\alpha(x, x') = \min\!\bigl\{1, \pi(x') / \pi(x)\bigr\}$. Start from $(0, 0)$, tune $\sigma$ until the empirical acceptance rate approaches the optimum. The banana’s curved ridge is the difficulty — too-small $\sigma$ leaves the chain sticky at the current point; too-large $\sigma$ overshoots the ridge and rejects almost every proposal.

For a $d$-dimensional Gaussian target with a Gaussian random-walk proposal of scale $\sigma$, Roberts–Gelman–Gilks 1997 proved the asymptotically optimal proposal covariance is $\Sigma = (2.38^2 / d) \, \Sigma_\pi$, where $\Sigma_\pi$ is the target covariance. The induced optimal acceptance rate is 0.234 for $d \ge 5$ and 0.44 for $d = 1$. ROB2004 Ch. 7 gives the full derivation. For Bayesian applications, the natural choice is $\Sigma \propto \hat{\mathcal{I}}(\hat\theta_{\text{MLE}})^{-1} / n$, using Topic 14 Theorem 4 — MLE asymptotic normality tells us the posterior is approximately Gaussian at large $n$, and its covariance is the inverse Fisher information scaled by $1/n$.

When $q(x' \mid x) = q(x \mid x')$ — a symmetric proposal, like a Gaussian random walk — the ratio $q(x \mid x') / q(x' \mid x) = 1$ and $\alpha(x, x') = \min\!\bigl\{1, \pi(x')/\pi(x)\bigr\}$. This is the 1953 Metropolis recipe; Hastings 1970 generalized to asymmetric proposals.

The independence sampler uses a proposal $q(x' \mid x) = q(x')$ that doesn’t depend on the current state — typically a Laplace approximation to the posterior. Acceptance ratio becomes $\pi(x') q(x) / \bigl[\pi(x) q(x')\bigr]$, with the importance-ratio-like structure familiar from importance sampling. Works well when $q$ is a good global approximation; catastrophic when it misses modes.

Target acceptance 20–40% for random-walk Metropolis on smooth unimodal posteriors; larger is over-proposing, smaller is over-rejecting. Tune the proposal covariance to $\Sigma \propto \hat{\mathcal{I}}^{-1}/n$ where available; otherwise use the empirical covariance of a pilot chain. The tuner component below lets you feel the acceptance-rate-vs-scale trade-off directly.

Geometric ergodicity — the rate of convergence to $\pi$ is exponential — requires compatible proposal and target tails (ROB2004 §7.3). Polynomial-tailed targets (Student-$t$ with few degrees of freedom, scale-mixture distributions) violate geometric ergodicity under Gaussian proposals; §26.8 Rem 21 names the standard mitigation (heavier-tailed proposals or reparameterization).

Let $\boldsymbol\theta = (\theta_1, \ldots, \theta_d) \in \mathbb{R}^d$ with joint target $\pi(\boldsymbol\theta)$. The full conditional of coordinate $i$ is

$\pi(\theta_i \mid \boldsymbol\theta_{-i}) \;=\; \frac{\pi(\boldsymbol\theta)}{\int \pi(\boldsymbol\theta) \, d\theta_i},$

where $\boldsymbol\theta_{-i} = (\theta_1, \ldots, \theta_{i-1}, \theta_{i+1}, \ldots, \theta_d)$. The systematic-scan Gibbs sampler iterates: at iteration $t+1$, for $i = 1, \ldots, d$ in sequence, draw $\theta_i^{(t+1)} \sim \pi\bigl(\theta_i \mid \theta_1^{(t+1)}, \ldots, \theta_{i-1}^{(t+1)}, \theta_{i+1}^{(t)}, \ldots, \theta_d^{(t)}\bigr)$.

The systematic-scan Gibbs sampler preserves $\pi$. Equivalently: each coordinate-wise sub-step is an MH step with acceptance probability identically 1.

Let $(X, Y) \sim \mathcal{N}\!\bigl(\mathbf{0}, \Sigma\bigr)$ with \Sigma = \begin{psmallmatrix}1 & \rho \\ \rho & 1\end{psmallmatrix}. Topic 8 Theorem 3 — the conditional-MVN formula — gives the full conditionals in closed form: $X \mid Y = y \sim \mathcal{N}(\rho y, 1 - \rho^2)$ and $Y \mid X = x \sim \mathcal{N}(\rho x, 1 - \rho^2)$. Gibbs alternates between these two draws. For $\rho$ near 1, consecutive draws are highly correlated — axis-aligned updates crawl along the diagonal — and mixing slows. This is the canonical Gibbs failure mode that motivates block updates, reparameterization, or HMC.

When a posterior splits into blocks whose full conditionals are each conjugate, Gibbs becomes the composition of closed-form Topic 25 updates. The five canonical pairs from §25.5 Examples 2–5 become Gibbs building blocks without re-derivation:

• $\theta \mid \mathbf{y}, \sigma^2 \sim \pi(\theta \mid \mathbf{y}, \sigma^2)$ via Normal-Normal conjugacy (Ex 2);
• $\sigma^2 \mid \mathbf{y}, \theta \sim \mathrm{Inv\text{-}Gamma}(\cdot, \cdot)$ via Normal-Normal-IG (Ex 3);
• Poisson rates $\lambda_j$ via Gamma-Poisson (Ex 4);
• Multinomial probabilities via Dirichlet-Multinomial (Ex 5).

Each sub-step is a single conjugate update; Gibbs glues them together. The composition preserves the joint posterior by Thm 2.

Systematic-scan Gibbs updates coordinates in a fixed order $1, 2, \ldots, d$ each sweep. Random-scan Gibbs picks a coordinate uniformly at random each step. Both preserve $\pi$; systematic-scan is the standard practical choice because it guarantees every coordinate is updated each sweep. Random-scan is easier to analyze theoretically (the transition kernel becomes symmetric in a way systematic-scan is not).

When coordinates are strongly correlated, one-at-a-time Gibbs mixes slowly (§26.3 Ex 3 with $\rho \to 1$). Block Gibbs groups correlated coordinates and updates them jointly from the joint full conditional — which must itself be tractable. For bivariate Normal with unknown covariance, drawing $(\theta_1, \theta_2)$ as a block from their joint conditional is a simple MVN draw and mixes far better than alternating $\theta_1 \mid \theta_2$ and $\theta_2 \mid \theta_1$.

Collapsed Gibbs integrates out nuisance parameters analytically before sampling. The canonical example is Latent Dirichlet Allocation: Topic 8 §8.10 Theorem 6 gives the Dirichlet-Multinomial conjugacy that lets us integrate out topic proportions $\boldsymbol\theta_d$ and word distributions $\boldsymbol\phi_k$ in closed form, leaving only the discrete topic assignments to sample via Gibbs. The state space shrinks dramatically and mixing improves — when the integrated-out parameters exist in closed form.

Let $\pi(q)$ be the target density on $\mathbb{R}^d$ ($q$ = position, following NEA2011). Augment with momentum $p \in \mathbb{R}^d$ and joint density

$\tilde\pi(q, p) \;=\; \pi(q) \, \mathcal{N}(p; \mathbf{0}, M),$

where $M \succ 0$ is the mass matrix. The Hamiltonian is $H(q, p) = U(q) + K(p)$ with $U(q) = -\log \pi(q)$ and $K(p) = \tfrac{1}{2} p^\top M^{-1} p$; so $\tilde\pi(q, p) \propto \exp(-H(q, p))$. The leapfrog integrator approximates the Hamiltonian flow over time $L\epsilon$ via $L$ steps of size $\epsilon$:

$p^{(k+1/2)} \;=\; p^{(k)} - \tfrac{\epsilon}{2} \nabla U(q^{(k)}),$

$q^{(k+1)} \;=\; q^{(k)} + \epsilon M^{-1} p^{(k+1/2)},$

$p^{(k+1)} \;=\; p^{(k+1/2)} - \tfrac{\epsilon}{2} \nabla U(q^{(k+1)}).$

The HMC iteration — resample $p \sim \mathcal{N}(\mathbf{0}, M)$, run $L$ leapfrog steps, flip momentum, MH-accept — preserves $\tilde\pi$ on the extended state space $(q, p)$. Consequently $\pi$ is preserved on the marginal $q$-space.

Target $\pi(q) \propto \exp\!\bigl(-\tfrac12\bigl[(1-q_1)^2 + 10(q_2 - q_1^2)^2\bigr]\bigr)$. Gradient $\nabla U$ is smooth and cheap to evaluate. Run HMC with $\epsilon = 0.05$, $L = 25$, mass matrix $M = I$. Starting from $(0, 0)$, HMC flows along the banana ridge — a single trajectory can traverse from one end to the other, whereas random-walk Metropolis with any reasonable scale needs thousands of steps for the same distance. The acceptance rate comfortably exceeds 0.6 under these settings; the 0.234 / 0.44 random-walk rule does not apply.

Symplectic integrators conserve a modified Hamiltonian $\tilde H = H + O(\epsilon^2)$ exactly, so while $H$ oscillates by $O(\epsilon^2)$ per step, the oscillation stays bounded as $L$ grows. Non-symplectic methods (e.g., Runge–Kutta applied to Hamilton’s equations) suffer energy drift that grows linearly with $L$ — the chain would reject almost every proposal at long horizons. NEA2011 §5.4 derives the leapfrog modified Hamiltonian in full.

Standard HMC practice targets acceptance 0.65–0.85 — higher than random-walk’s 0.234 / 0.44 because HMC proposals are near-deterministic and should land close in energy. Too-low acceptance signals $\epsilon$ is too large (energy drift swamps the target density); too-high acceptance signals $\epsilon$ is too small (we’re under-exploring). Divergences — iterations where $|H_{\text{end}} - H_{\text{start}}|$ exceeds a threshold — indicate leapfrog has “fallen off” the energy manifold; in Stan / PyMC these are reported and usually signal a reparameterization opportunity.

The mass matrix $M$ defines the momentum distribution and hence the geometry HMC sees. The theoretically motivated choice is $M = \mathcal{I}(\theta_0)$ — Fisher information at the posterior mode — because then the transformed target is locally standard Normal and HMC’s Gaussian-momentum default is well-conditioned. This connects to Topic 14 Theorem 4 (MLE asymptotic normality; the posterior is approximately Gaussian with covariance $\mathcal{I}^{-1}/n$ at large $n$) and to Topic 25 §25.8’s Remark 16 Laplace-at-MAP approximation, which provides the standard HMC warm-start in practice: initialize at $\hat\theta_{\text{MAP}}$ and set $M = \hat{\mathcal{I}}(\hat\theta_{\text{MAP}})$ from the Laplace Hessian.

Given forward and backward trajectory endpoints $q_+, p_+$ and $q_-, p_-$ (with $q_+ = q_-$ initially), a U-turn is detected when either of

$\langle q_+ - q_-,\; p_+\rangle \;<\; 0 \quad \text{or} \quad \langle q_+ - q_-,\; p_-\rangle \;<\; 0$

holds. Intuitively: the trajectory has turned far enough that its endpoint momenta no longer point outward along the displacement vector. NUTS terminates doubling when a U-turn is detected on any subtree.

The NUTS sampler with dual-averaging step-size adaptation and slice-sampling-based subtree validity preserves $\pi$ on the extended state space $(q, p)$.

The canonical NUTS stress test is Rubin’s 1981 8-schools hierarchical model (GEL2013 §5.5). Eight observed treatment effects $y_j \sim \mathcal{N}(\theta_j, \sigma_j^2)$ with school-level effects $\theta_j \sim \mathcal{N}(\mu, \tau^2)$ and hyperpriors on $(\mu, \tau)$. At small $\tau$, the posterior develops Neal’s funnel — a narrow neck where $\tau$ concentrates near zero. NUTS with Stan’s default settings handles this acceptably; non-centered reparameterization (§26.8 Rem 23, Topic 28) improves mixing further. Full applied workflow in §26.9 Example 11.

During a warmup phase (~1000 iterations), NUTS runs a dual-averaging scheme that targets acceptance around 0.8 by adjusting $\epsilon$ on the fly (HOF2014 §3.2). After warmup, $\epsilon$ is frozen and sampling begins. The adaptation is optimization-free — just Robbins–Monro stochastic approximation on the acceptance deviation — and converges reliably under Hoffman–Gelman’s conditions.

Stan (since 2014), PyMC (since 2016), and NumPyro all ship NUTS as their default sampler. The practical implication for practitioners: specify the model in the DSL, let NUTS sample, check $\hat R$ and ESS. The §26.9 workflow assumes this; implementation internals are deferred to formalML: Probabilistic Programming (§26.10 Rem 33).

Suppose we run $M$ chains of length $N$ from dispersed starting points, on a scalar summary of the chain state. Let $\bar x_m$ and $s_m^2$ be the $m$-th chain’s sample mean and (unbiased) sample variance, and $\bar{\bar x}$ the grand mean across chains. Define

$B \;=\; \frac{N}{M - 1} \sum_{m=1}^M \bigl(\bar x_m - \bar{\bar x}\bigr)^2,\qquad W \;=\; \frac{1}{M} \sum_{m=1}^M s_m^2.$

The Gelman–Rubin statistic is

$\hat R \;=\; \sqrt{\frac{N-1}{N} \,+\, \frac{B}{N\,W}}.$

For a single chain of length $N$ with lag-$t$ autocorrelation $\rho_t$, the effective sample size is

$N_{\text{eff}} \;=\; \frac{N}{1 \,+\, 2\sum_{t=1}^\infty \rho_t}.$

In practice the infinite sum is truncated by Geyer’s initial-positive-sequence estimator: group consecutive lags into pairs $\Gamma_m = \rho_{2m} + \rho_{2m+1}$, accumulate until the first non-positive pair, and stop. For AR(1) with parameter $\varphi$, the theoretical $\tau = (1 + \varphi) / (1 - \varphi)$, so $\varphi = 0.9$ gives $\tau = 19$ and $N_{\text{eff}} \approx N / 19$.

Initialize four chains at positions $\\{-3, -1, 1, 3\\}$ targeting $\mathcal{N}(0, 1)$ with a random-walk Metropolis proposal at the RGG-1997 optimum. At iteration 100, chains are still concentrated near their starts — $\hat R \approx 1.34$ signals mixing has not occurred. By iteration 2000, chains have interpenetrated; $\hat R$ has dropped below 1.01. The dashboard component below lets you scrub through this collapse interactively.

Generate an AR(1) chain $X_{t+1} = \varphi X_t + \sqrt{1 - \varphi^2} \, Z_t$ with $Z_t \sim \mathcal{N}(0, 1)$ iid. This has stationary variance 1, but lag-1 autocorrelation $\varphi$. For $N = 5000$ samples at $\varphi = 0.9$, the Geyer-IPS estimator gives $N_{\text{eff}} \approx 263$ — the 5000 correlated samples carry the information of roughly 263 iid samples. At $\varphi = 0.99$, $N_{\text{eff}} / N$ falls below 1%.