Markov Chain Monte Carlo

Definition

Markov Chain Monte Carlo (MCMC) is a class of algorithms for drawing samples from a target distribution $\pi(\theta)$ by constructing a Markov chain whose stationary distribution is $\pi(\theta)$ . Because the chain must mix before its output approximates the target, MCMC trades the independence of direct Monte Carlo for the ability to sample from distributions that have no closed-form and resist other approximations. After a burn-in period the chain's draws are used to estimate posterior means, variances, quantiles, and marginal densities.

Key Ideas

A Markov chain with transition kernel $P(\theta, \cdot)$ has $\pi$ as its stationary distribution iff $\pi P = \pi$ . Sufficient (not necessary) condition: detailed balance $\pi(\theta)P(\theta, \theta') = \pi(\theta')P(\theta', \theta)$ — satisfied by construction in the Metropolis-Hastings (MH) accept/reject step.
Ergodicity (irreducibility + aperiodicity + positive recurrence) guarantees that time averages converge to space averages: $N^{-1}\sum_{t=1}^N h(\theta^{(t)}) \xrightarrow{a.s.} E_\pi[h(\theta)]$ .
The Central Limit Theorem for MCMC (Tierney 1994): under a geometric drift condition, $\sqrt{N}(\bar h - E_\pi[h]) \to \mathcal{N}(0, \sigma^2_h)$ where $\sigma^2_h = \sum_{k=-\infty}^\infty \text{Cov}(h(\theta^{(0)}), h(\theta^{(k)}))$ — the asymptotic variance includes all autocorrelations, not just the marginal variance.
Peskun (1973) optimality: among all kernels sharing the same proposal, the standard MH acceptance rule minimises the asymptotic variance of any estimator (it maximises the off-diagonal mass in the transition matrix subject to detailed balance).

Algorithm Taxonomy

Gibbs Sampler

Cycle through each block $\theta_i$ drawing from its full conditional $\pi(\theta_i | \theta_{-i}, y)$ . A special case of MH with acceptance probability 1 (when conditionals are available in closed form). See Gibbs Sampler for full details: convergence proof via fixed-point equation, Rao-Blackwell density estimator, blocked sampling, Adaptive Rejection Sampling (ARS) for log-concave conditionals, data augmentation strategy.

Metropolis-Hastings (MH)

Draw a proposal $\theta^*$ from $q(\theta^* | \theta^{(t)})$ and accept with probability $\alpha = \min[1,\, \pi(\theta^*)q(\theta^{(t)}|\theta^*)/(\pi(\theta^{(t)})q(\theta^*|\theta^{(t)}))]$ . Requires only the ability to evaluate $\pi$ up to a normalising constant. See Metropolis-Hastings Algorithm for full details: five candidate families, Product of Kernels principle, Metropolis-within-Gibbs, prefetching (Strid 2010).

Metropolis-within-Gibbs (Hybrid)

Embed MH steps for non-conjugate blocks inside a Gibbs sweep. The combined kernel has $\pi$ as its invariant distribution by the Product of Kernels principle: if $P_1$ and $P_2$ each leave $\pi$ invariant, so does $P_1 P_2$ . Used in Time-Varying Parameter VAR (TVP-VAR), Markov-Switching GARCH (MS-GARCH) (Das-Yoo 2004), hierarchical Seemingly Unrelated Regression (SUR) (Chib-Greenberg 1995b), and stochastic volatility samplers.

Hamiltonian Monte Carlo (HMC)

Augment $q$ with fictitious momentum $p \sim \mathcal{N}(0, M)$ . Define $H(q,p) = U(q) + K(p)$ with $U = -\log\pi(q)$ and $K = p^T M^{-1} p/2$ . Each HMC step: (1) resample $p$ ; (2) run $L$ leapfrog steps with stepsize $\varepsilon$ to get proposal $(q^*, p^*)$ ; (3) accept with probability $\min(1, \exp(H(q,p)-H(q^*,p^*)))$ . The leapfrog integrator is symplectic — it exactly preserves phase-space volume, eliminating the Jacobian correction, and approximately conserves $H$ so distant proposals are accepted with high probability (~65% optimal). Scaling advantage: the total arithmetic cost per independent sample scales as $O(d^{5/4})$ for HMC versus $O(d^2)$ for random-walk Metropolis — a $d^{3/4}$ improvement (HMC uses $O(d^{1/4})$ gradient evaluations per sample, RWM $O(d)$ density evaluations, each costing $O(d)$ ). Requires differentiable $\log\pi$ ; cannot handle discrete parameters directly. The No-U-Turn Sampler (NUTS, Hoffman-Gelman 2011) automates stepsize and trajectory-length tuning and is the default in Stan. See Hamiltonian Monte Carlo for full details: leapfrog derivation, three key properties (reversibility, volume preservation, energy conservation), mass matrix preconditioning, variations (windowed HMC, partial momentum refreshment, tempered trajectories), and Riemannian extensions.

Reversible Jump MCMC (RJMCMC)

See Reversible Jump MCMC for full details. Extends MH to variable-dimension parameter spaces (Green 1995): proposes moves between models of different sizes via dimension-matching bijections $(\theta,u)\to(\theta',u')$ with Jacobian correction $|\partial(\theta',u')/\partial(\theta,u)|$ . Used for Bayes factor computation when model dimension varies (Han-Carlin 2001 benchmark: RJ $\hat B_{21}=67.83$ for AIDS CD4 model where Gibbs-based methods fail) and for nonparametric regression with unknown dimension (Perron-Mengersen 2001: birth/death/move over spline knot count $r$ ).

Particle Methods (Sequential MC)

Propagate a weighted population of particles forward in time; resample to avoid weight degeneracy. Outside the scope of current wiki sources but referenced in Carvalho-Johannes-Lopes-Polson (2010) particle learning.

Convergence Diagnostics

Gelman-Rubin $\hat R$

Run $M \geq 2$ chains from overdispersed starting points. Compare within-chain variance $W$ and between-chain variance $B$ : $\hat R = \sqrt{(n-1)/n + B/(nW)}$ . Values near 1 indicate convergence; $\hat R > 1.1$ signals poor mixing. Requires multiple parallel chains.

Geweke Spectral Diagnostic (1992)

Split a single chain into an early fraction (first 10%) and a late fraction (last 50%). Compute a $z$ -statistic for the difference in means using spectral variance estimates. Under convergence $z \sim \mathcal{N}(0,1)$ . Reports the Relative Numerical Efficiency (RNE) = $\sigma^2_{\text{iid}} / \sigma^2_{\text{MCMC}}$ ; RNE $<$ 1 indicates autocorrelation loss.

Raftery-Lewis (1992)

Given a desired quantile $q$ , accuracy $r$ , and probability $s$ , derives the minimum number of MCMC draws $N_{\min}$ and required thinning $k$ via a two-state Markov chain approximation. Useful for pre-run planning; implemented in CODA.

Zellner-Min $\text{GSC}^2$ (1995)

Three criteria ( $\text{ARC}^2$ anchored ratio, $\text{DC}^2$ difference, $\text{RC}^2$ ratio) that detect convergence to the wrong distribution — the failure mode standard diagnostics miss. Core identity: $p(\alpha|D) p(\beta|\alpha, D) = p(\beta|D) p(\alpha|\beta, D)$ ; violations signal non-stationarity or multimodality. Applied to bimodal posteriors and misspecified models.

Efficiency and Thinning

Inefficiency factor $\hat\tau = \text{Var}(\hat h_M)/(s^2/M)$ (Chib 2001, eq. 13): ratio of numerical variance to the variance under i.i.d. sampling. Estimated by batch means (batch size chosen so first-order autocorrelation of batch means $< 0.05$ ); inverse is Geweke's Relative Numerical Efficiency. Values $\hat\tau \gg 1$ signal severe autocorrelation and motivate improved blocking or reparameterization.

Thinning is never beneficial (MacEachern-Berliner 1994): for any stationary Markov chain, the variance of the full-chain mean is ≤ the variance of any systematic subsample mean. The only legitimate reason to thin is storage or downstream processing cost. Preferred alternatives: batching and spectral variance estimation. See Gibbs Sampler for the proof and the stationary time-series trap.

Blocking: drawing multiple parameters jointly (block Gibbs) reduces the effective autocorrelation by cutting the number of successive conditioning steps. Critical for high-dimensional problems where single-site updates mix slowly.

Rao-Blackwellisation: replace a Monte Carlo estimate $N^{-1}\sum h(\theta^{(t)})$ with its conditional expectation $N^{-1}\sum E[h(\theta)|\theta_{-i}^{(t)}]$ . Always reduces variance (Casella-George 1992); especially useful for binary indicators in mixture/variable-selection models.

Open Questions

Choosing between Gibbs and MH for a given block is context-dependent: Gibbs requires closed-form conditionals; MH requires tuning of the proposal scale for ~25% acceptance rate in high dimensions.
Diagnosing convergence in multimodal posteriors remains hard; $\text{GSC}^2$ criteria (Zellner-Min 1995) are the most systematic approach in the current wiki but require access to paired conditional evaluations.
Parallel computation (prefetching, particle parallelism) can reduce wall-clock time but complicates the standard ergodicity arguments.