MCMC Concepts — Markov chain Monte Carlo (MCMC) concepts - Statistics Graduate Level

SECTION 1

Historical Context & Motivation

One of the central challenges in Bayesian statistics is computing the posterior distribution. Bayes' theorem gives us the posterior as proportional to the likelihood times the prior, but the normalizing constant—the marginal likelihood—requires integrating over the entire parameter space. For models with more than a handful of parameters, this integral is analytically intractable and numerically prohibitive using standard quadrature methods. Markov chain Monte Carlo (MCMC) methods were developed precisely to overcome this obstacle, enabling researchers to draw dependent samples from posterior distributions without computing the normalizing constant directly.

The intellectual roots of MCMC span physics, mathematics, and statistics. The foundational ideas emerged from physicists who needed to simulate the thermodynamic behavior of many-particle systems at Los Alamos National Laboratory, and the statistical community later recognized the profound generality of these simulation techniques for inference. Today, MCMC is arguably the single most important computational tool in Bayesian analysis, underpinning everything from hierarchical models to machine learning and computational biology.

1949–1953 • The Metropolis Algorithm

Nicholas Metropolis, along with Arianna and Marshall Rosenbluth and Augusta and Edward Teller, published the original Metropolis algorithm for simulating equations of state using random walks on an energy landscape. The paper, targeting statistical mechanics, laid the groundwork for all subsequent MCMC methods.

1970 • Hastings Generalization

W. Keith Hastings generalized Metropolis's symmetric proposal to allow asymmetric proposal distributions, creating the Metropolis-Hastings algorithm. This broadened the applicability of MCMC to a far wider class of target distributions.

1984 • Gibbs Sampler Formalized

Stuart and Donald Geman introduced the Gibbs sampler for image restoration, sampling each variable from its full conditional distribution. This special case of Metropolis-Hastings proved enormously practical for multivariate Bayesian models.

1990 • MCMC Enters Mainstream Statistics

Gelfand and Smith's landmark paper demonstrated the power of Gibbs sampling for Bayesian computation, sparking an explosion of MCMC applications across statistics, epidemiology, ecology, and the social sciences.

2011–Present • Hamiltonian Monte Carlo & Modern Samplers

Hamiltonian Monte Carlo (HMC) and the No-U-Turn Sampler (NUTS) were popularized through Stan, dramatically improving efficiency in high-dimensional parameter spaces. Modern MCMC research focuses on scalability, adaptation, and diagnostics.

The fundamental question MCMC addresses is: how can we characterize a probability distribution from which we cannot directly sample? By constructing a Markov chain whose stationary distribution is the target posterior, MCMC transforms an integration problem into a simulation problem—one that scales to hundreds or even thousands of dimensions.

SECTION 2

Core Principles & Definitions

Understanding MCMC requires grasping several interlocking ideas from probability theory and stochastic processes. At its heart, the method relies on the theory of Markov chains: memoryless stochastic processes whose future state depends only on the current state, not the path taken to reach it. When such a chain is carefully designed so that its long-run equilibrium matches the posterior distribution of interest, we can treat the chain's realized values as (correlated) samples from the posterior.

Markov Property

The probability of transitioning to the next state depends only on the current state, not on any previous states. Formally: P(θ^(t+1) | θ^(t), θ^(t−1), …) = P(θ^(t+1) | θ^(t)).

Stationary Distribution

A distribution π(θ) is stationary for a chain if, once the chain's marginal distribution equals π, it remains π for all subsequent iterations. MCMC constructs chains whose unique stationary distribution is the target posterior p(θ | y).

Ergodicity

An ergodic Markov chain guarantees that time averages converge to expectations under the stationary distribution, regardless of the starting point. This requires the chain to be irreducible (can reach any state) and aperiodic (does not cycle deterministically).

Detailed Balance

A sufficient condition for π to be stationary is the detailed balance (reversibility) condition: π(θ)T(θ → θ′) = π(θ′)T(θ′ → θ) for all states θ, θ′. Most MCMC algorithms are designed to satisfy this.

Burn-in & Mixing

Burn-in refers to discarding early samples before the chain has converged to the stationary distribution. Mixing describes how rapidly the chain explores the full support of the target; poor mixing yields highly autocorrelated samples.

✦

KEY TAKEAWAY

Think of MCMC as exploring an unfamiliar city at night. You wander from block to block, deciding at each intersection whether to turn based only on the appeal of the next block—not your entire route history (Markov property). Over a long walk, you naturally spend more time in popular neighborhoods (high-density regions of the posterior). Eventually your trajectory traces out a map of the city's popularity, even though you never saw the full map up front. This is precisely how MCMC converts a random walk into a faithful representation of a target distribution.

SECTION 3

Visual Explanation — Anatomy of an MCMC Sampler

The following diagram illustrates the core loop of the Metropolis-Hastings algorithm, the most general MCMC framework. At each iteration the sampler occupies a current state θ^(t), proposes a candidate state θ* from a proposal distribution q(θ* | θ^(t)), computes an acceptance ratio α, and then either moves to θ* or remains at θ^(t). Repeating this process generates the Markov chain whose empirical distribution converges to the posterior.

The Metropolis-Hastings loop: starting from an initial value θ⁽⁰⁾, the algorithm iteratively proposes candidates, evaluates the acceptance ratio α, and stochastically accepts or rejects moves. The stored sequence {θ^(t)} forms the Markov chain whose empirical distribution approximates the posterior.

Notice that the acceptance ratio α involves only the ratio of posterior densities π(θ*)/π(θ^(t)), which means the normalizing constant cancels. This is the key insight that makes MCMC practical: we only need the unnormalized posterior, i.e., the product of the likelihood and the prior. The proposal distribution q governs where the chain tries to move and thus critically affects the efficiency of the sampler. A proposal that is too narrow produces small steps and high acceptance but slow exploration; a proposal that is too wide produces large jumps that are mostly rejected.

SECTION 4

Mathematical Framework

The mathematical foundation of MCMC rests on constructing a transition kernel T(θ′ | θ) that leaves the target distribution π(θ) = p(θ | y) invariant. Below we formalize the key equations that govern the Metropolis-Hastings algorithm and establish its theoretical validity.

BAYES' THEOREM — POSTERIOR

p(θ | y) = p(y | θ) · p(θ) / p(y) ∝ p(y | θ) · p(θ)

Here p(y | θ) is the likelihood, p(θ) is the prior, and p(y) = ∫ p(y | θ)p(θ)dθ is the marginal likelihood (normalizing constant). MCMC avoids computing p(y) altogether by working with the unnormalized posterior.

METROPOLIS-HASTINGS ACCEPTANCE RATIO

α(θ*, θ⁽ᵗ⁾) = min( 1, [π(θ*) · q(θ⁽ᵗ⁾ | θ*)] / [π(θ⁽ᵗ⁾) · q(θ* | θ⁽ᵗ⁾)] )

π(θ) is the target density (unnormalized posterior), q(θ* | θ^(t)) is the proposal density of jumping from θ^(t) to θ*. When the proposal is symmetric (q(a|b) = q(b|a)), the ratio simplifies to min(1, π(θ*)/π(θ^(t)))—this is the original Metropolis algorithm.

DETAILED BALANCE CONDITION

π(θ) · T(θ′ | θ) = π(θ′) · T(θ | θ′) for all θ, θ′

This reversibility condition ensures that the probability flux from θ to θ′ equals the flux from θ′ to θ under the stationary distribution π. Satisfying detailed balance is sufficient (though not necessary) for π to be the chain's stationary distribution.

ERGODIC THEOREM (LAW OF LARGE NUMBERS FOR MCMC)

(1/N) Σₜ₌₁ᴺ g(θ⁽ᵗ⁾) → E_π[g(θ)] = ∫ g(θ) π(θ) dθ as N → ∞

For any integrable function g(θ), the sample average from the chain converges almost surely to the posterior expectation. This result underpins all MCMC-based estimation: posterior means, variances, quantiles, and predictive distributions are all computed as sample averages over the chain.

A critical practical consideration is the effective sample size (ESS). Because consecutive MCMC draws are autocorrelated, N draws from the chain carry less information than N independent draws. The ESS quantifies this: ESS = N / (1 + 2 × Σ_k=1^∞ ρ_k), where ρ_k is the lag-k autocorrelation of the chain. High autocorrelation means low ESS, demanding longer runs or better samplers.

SECTION 5

MCMC Algorithm Variants & Comparison

Several important MCMC algorithms have been developed, each tailored to different problem structures. The three most widely used are the Metropolis-Hastings algorithm, the Gibbs sampler, and Hamiltonian Monte Carlo (HMC). Understanding when to apply each algorithm is essential for efficient Bayesian computation.

Three MCMC algorithms sampling from the same correlated bivariate normal target. Random Walk MH (violet) takes many small, correlated steps. The Gibbs sampler (cyan) moves along coordinate axes, struggling with correlation. HMC (amber) uses gradient information to make large, efficient leaps through the parameter space.

Property	Metropolis-Hastings	Gibbs Sampler	HMC / NUTS
Proposal mechanism	User-specified q(θ* \| θ)	Full conditional distributions	Hamiltonian dynamics + leapfrog integrator
Accept/reject?	Yes (α computed each step)	Always accepted (special case of MH with α = 1)	Yes (but acceptance ≈ 65–95%)
Gradient needed?	No	No	Yes (∇ log π)
Scalability (dim.)	Poor (O(d²) scaling)	Moderate (d conditionals per sweep)	Good (O( $d^{5/4}$ ) scaling)
Best suited for	Low-dimensional, general targets	Conjugate models, conditional independence	High-dimensional, continuous parameters

Comparison of the three major MCMC algorithm families

SECTION 6

Worked Example — Metropolis Algorithm for a Normal Posterior

Suppose we observe a single data point y = 5 from a Normal(μ, σ² = 4) likelihood, and we place a Normal(0, 10²) prior on μ. We want to sample from the posterior p(μ | y = 5) using the Metropolis algorithm with a symmetric Normal(μ^(t), 1) proposal. We trace through three iterations starting from μ⁽⁰⁾ = 0.

Metropolis Sampling: 3 Iterations

1
Step 1 — Set up the unnormalized posterior

The unnormalized posterior is π(μ) ∝ exp(−(y − μ)² / (2 × 4)) × exp(−μ² / (2 × 100)). Since the proposal is symmetric (Normal), the acceptance ratio simplifies to α = min(1, π(μ*)/π(μ^(t))). It is convenient to work on the log scale: log π(μ) = −(5 − μ)² / 8 − μ² / 200 + const.

log π(μ) = −(5 − μ)²/8 − μ²/200
2
Step 2 — Iteration 1: propose and evaluate

Current state: μ⁽⁰⁾ = 0. Suppose the proposal yields μ* = 1.2. Compute log π(0) = −(5)²/8 − 0/200 = −3.125. Compute log π(1.2) = −(3.8)²/8 − (1.44)/200 = −1.805 − 0.0072 = −1.8122. The log acceptance ratio is log α = −1.8122 − (−3.125) = 1.3128, so α = min(1, e^1.3128) = 1. Since α = 1, we accept with certainty.

μ⁽¹⁾ = 1.2 (ACCEPTED, α = 1)
3
Step 3 — Iteration 2: propose and evaluate

Current state: μ⁽¹⁾ = 1.2. Suppose μ* = 6.8. Compute log π(6.8) = −(−1.8)²/8 − (46.24)/200 = −0.405 − 0.2312 = −0.6362. Log ratio = −0.6362 − (−1.8122) = 1.176, so α = 1. Accept.

μ⁽²⁾ = 6.8 (ACCEPTED, α = 1)
4
Step 4 — Iteration 3: propose and evaluate

Current state: μ⁽²⁾ = 6.8. Suppose μ* = 10.5. Compute log π(10.5) = −(−5.5)²/8 − (110.25)/200 = −3.78125 − 0.55125 = −4.3325. Log ratio = −4.3325 − (−0.6362) = −3.6963. So α = e^−3.70 ≈ 0.025. Draw u ~ Uniform(0,1). Suppose u = 0.41. Since u > α, we reject the proposal.

μ⁽³⁾ = 6.8 (REJECTED, stays at previous value)
5
Step 5 — Interpretation

After many more iterations, the histogram of sampled μ values will approximate the posterior distribution. For this conjugate model, the exact posterior is Normal(μ_post, σ²_post) where μ_post = (5/4 + 0/100) / (1/4 + 1/100) = 1.25/0.26 ≈ 4.808 and σ²_post = 1/0.26 ≈ 3.846. Running the chain for thousands of iterations and discarding burn-in would yield a histogram centered near 4.81.

Exact posterior: Normal(4.808, 3.846). MCMC chain converges to this distribution.

SECTION 7

Convergence Diagnostics & Practical Considerations

A crucial aspect of MCMC practice is assessing whether the chain has converged to its stationary distribution and is mixing well. Unlike optimization, where a single point estimate is sought, MCMC requires that the entire distribution of samples faithfully represents the posterior. Several diagnostic tools have been developed to detect common pathologies such as non-convergence, poor mixing, and multimodality.

Diagnostic	What It Measures	Rule of Thumb / Threshold
Trace plots	Visual inspection of parameter values across iterations. Look for 'hairy caterpillar' appearance: rapid oscillation around a stable mean.	No trends, drifts, or stuck periods visible
R̂ (Gelman-Rubin)	Ratio of between-chain to within-chain variance across multiple parallel chains. Values near 1 indicate convergence.	R̂ < 1.01 (modern recommendation; historically < 1.1)
Effective Sample Size (ESS)	Estimates the number of independent draws equivalent to the autocorrelated chain. Low ESS indicates high autocorrelation.	ESS > 400 per parameter for reliable posterior summaries
Autocorrelation function (ACF)	Plots autocorrelation at successive lags. Ideally drops to zero within a few lags.	Rapid decay to zero; persistent correlation signals poor mixing
Divergent transitions (HMC)	Indicates the leapfrog integrator has failed to track Hamiltonian dynamics, often due to high curvature in the posterior.	Zero divergences is ideal; any divergences require reparameterization

Common MCMC convergence diagnostics

✦

KEY TAKEAWAY

MCMC diagnostics are like quality-control tests on a manufacturing line. Just as a factory checks that products meet specifications before shipping, a statistician must check that the MCMC chain has converged and mixed properly before using the samples for inference. Never trust MCMC output without running multiple chains and checking R̂, ESS, and trace plots. Failing to diagnose convergence issues can lead to severely biased posterior summaries and misleading conclusions.

SECTION 8

Connections to Advanced Theory & Modern Extensions

The foundational MCMC concepts covered so far connect directly to several active areas of research in computational statistics and machine learning. Understanding these connections provides both deeper theoretical insight and practical awareness of when standard MCMC may be insufficient.

Core MCMC Concept	Advanced Extension	Key Idea
Random Walk MH	Adaptive MCMC	Automatically tune proposal covariance during sampling using the chain's history (e.g., Haario et al.'s Adaptive Metropolis). Requires diminishing adaptation to preserve ergodicity.
Single-chain sampling	Parallel tempering	Run multiple chains at different 'temperatures' (flattened posteriors) and swap states. Helps chains escape local modes in multimodal posteriors.
Full-data likelihood	Stochastic gradient MCMC	Use mini-batch gradient estimates (SGLD, SGHMC) for scalability to massive datasets, combining stochastic optimization with Bayesian sampling.
MCMC for posterior sampling	Variational inference (VI)	An alternative to MCMC that approximates the posterior with a simpler family of distributions by optimizing a divergence measure. Faster but biased; MCMC is asymptotically exact.
Gibbs with conjugacy	Non-conjugate Gibbs via auxiliary variables	Introduce latent variables (e.g., Pólya-Gamma augmentation for logistic regression) to create conjugate full conditionals even in non-conjugate models.

From core MCMC to advanced extensions

An especially important distinction for practitioners is between MCMC and variational inference. While MCMC provides asymptotically exact samples from the posterior, variational inference trades exactness for speed by casting inference as optimization over a parametric family of approximate distributions. In modern Bayesian deep learning and large-scale applications, variational methods often serve as fast approximations, while MCMC remains the gold standard when accuracy is paramount—particularly for model comparison via marginal likelihoods and for detecting posterior multimodality.

💡

Looking Ahead

Recent developments in normalizing flows and neural MCMC use deep learning to learn efficient proposal distributions or transport maps, potentially combining the scalability of VI with the asymptotic exactness of MCMC. This remains an exciting frontier at the intersection of Bayesian statistics and machine learning.

SECTION 9

Practice Problems

Question 1: CONCEPTUAL

Explain why MCMC methods do not require computing the normalizing constant p(y) in Bayes' theorem. Specifically, identify where in the Metropolis-Hastings acceptance ratio the normalizing constant cancels and why this cancellation is the key that makes MCMC practical for complex models.

Show answer

The Metropolis-Hastings acceptance ratio involves the ratio π(θ*)/π(θ^(t)), where π(θ) = p(y | θ)p(θ)/p(y). Since p(y) appears in both the numerator and denominator, it cancels: π(θ*)/π(θ^(t)) = [p(y|θ*)p(θ*)] / [p(y|θ^(t))p(θ^(t))]. This means we only need to evaluate the unnormalized posterior—the product of the likelihood and the prior—at any point. Computing p(y) = ∫ p(y|θ)p(θ)dθ requires integrating over the entire parameter space, which is intractable for most models of practical interest. By avoiding this integral, MCMC enables Bayesian inference in models with hundreds or thousands of parameters.

Question 2: BASIC CALCULATION

Consider a Metropolis algorithm targeting a posterior π(θ) ∝ exp(−θ²/2) (a standard normal) with a symmetric proposal. At the current state θ^(t) = 1.5, the proposal generates θ* = 0.3. Compute the acceptance probability α. Would you accept or reject if u = 0.7?

Show answer

Since the proposal is symmetric, α = min(1, π(θ*)/π(θ^(t))) = min(1, exp(−0.3²/2) / exp(−1.5²/2)) = min(1, exp(−0.045)/exp(−1.125)) = min(1, exp(1.08)) = min(1, 2.945) = 1. The acceptance probability is 1, meaning we accept the proposal with certainty regardless of the uniform draw. Since u = 0.7 < 1 = α, we accept: θ^(t+1) = 0.3. This makes intuitive sense: the proposal moves from the tail (θ = 1.5) toward the mode (θ = 0), a region of higher density, so it is always accepted.

Question 3: INTERMEDIATE

A researcher runs a Gibbs sampler on a two-parameter model (μ, σ²) with strong posterior correlation between μ and σ². After 10,000 iterations with 2,000 burn-in, the trace plot for μ shows slow, wandering behavior and the ESS is only 85. (a) Diagnose the problem. (b) Propose two concrete strategies to improve the sampler's performance.

Show answer

(a) The Gibbs sampler updates each parameter conditionally on the other, moving along coordinate axes. When μ and σ² are highly correlated, the conditional distributions are narrow slices of the joint posterior, forcing the sampler to take many small steps—a phenomenon called slow mixing. The low ESS of 85 from 8,000 post-burn-in samples indicates extreme autocorrelation (autocorrelation time ≈ 94). (b) Strategy 1: Reparameterize the model to reduce posterior correlation. For instance, centering the data (subtracting the sample mean) can decorrelate μ and σ². Strategy 2: Switch to Hamiltonian Monte Carlo (e.g., via Stan's NUTS sampler), which uses gradient information to propose moves along the correlation structure rather than along coordinate axes, yielding much higher ESS per iteration.

Question 4: APPLIED

You are fitting a Bayesian logistic regression model with 50 predictors to a clinical trial dataset of 500 patients. You have access to Stan (which implements NUTS/HMC) and JAGS (which uses Gibbs sampling with Metropolis-within-Gibbs steps). (a) Which software would you recommend and why? (b) After running 4 chains of 2,000 iterations each in Stan, you observe 23 divergent transitions. What should you do?

Show answer

(a) Stan with NUTS/HMC is strongly preferred for this problem. Logistic regression has no conjugate full conditionals, so JAGS would resort to Metropolis-within-Gibbs (adaptive random walk updates for each coefficient), leading to poor mixing in 50 dimensions. HMC exploits the gradient of the log-posterior to make efficient, large moves through the 50-dimensional parameter space, typically yielding far higher ESS per second. (b) Divergent transitions in HMC indicate that the leapfrog integrator is encountering regions of high curvature in the posterior geometry. First, increase adapt_delta (e.g., from 0.8 to 0.95 or 0.99), which forces a smaller step size. If divergences persist, reparameterize: use a non-centered parameterization if there are hierarchical priors, or consider a QR decomposition of the design matrix to reduce collinearity. Finally, check for data issues such as complete separation, which creates a flat or improper posterior that HMC cannot handle.

Question 5: CRITICAL THINKING

Prove that the Metropolis-Hastings acceptance probability α(θ*, θ) = min(1, [π(θ*)q(θ | θ*)] / [π(θ)q(θ* | θ)]) satisfies the detailed balance condition with respect to the target distribution π(θ). Discuss what happens when the chain is reducible (i.e., can get trapped in a subset of the parameter space) and whether detailed balance alone guarantees convergence.

Show answer

The transition kernel of MH is T(θ′ | θ) = q(θ′ | θ)α(θ′, θ) for θ′ ≠ θ, plus a point mass at θ for rejection. To verify detailed balance, consider θ ≠ θ′. Without loss of generality, assume π(θ)q(θ′|θ) ≥ π(θ′)q(θ|θ′). Then α(θ′,θ) = π(θ′)q(θ|θ′)/[π(θ)q(θ′|θ)] and α(θ,θ′) = 1. We have π(θ)T(θ′|θ) = π(θ)q(θ′|θ) × [π(θ′)q(θ|θ′)/(π(θ)q(θ′|θ))] = π(θ′)q(θ|θ′). And π(θ′)T(θ|θ′) = π(θ′)q(θ|θ′) × 1 = π(θ′)q(θ|θ′). These are equal, confirming detailed balance. However, detailed balance alone is not sufficient for convergence. It guarantees that π is a stationary distribution, but uniqueness and convergence require ergodicity: the chain must be irreducible (can reach any set of positive π-measure) and aperiodic. If the chain is reducible—say the proposal can never reach certain regions—then the chain may converge to a different stationary distribution supported on a subset, yielding biased inference. This is why choosing a proposal with sufficiently broad support is critical in practice.

SUMMARY

Lesson Summary

Summary

Markov chain Monte Carlo (MCMC) methods solve the fundamental computational challenge of Bayesian statistics: sampling from posterior distributions whose normalizing constants are intractable. By constructing a Markov chain whose stationary distribution equals the target posterior, MCMC converts an integration problem into a simulation problem. The Metropolis-Hastings algorithm proposes candidate states and accepts them with a probability that ensures detailed balance, while the Gibbs sampler cycles through full conditional distributions, and Hamiltonian Monte Carlo leverages gradient information for efficient exploration of high-dimensional spaces.

Rigorous practice demands attention to convergence diagnostics: trace plots, R̂, and effective sample size (ESS) must all be checked before interpreting results. Key practical decisions include selecting the appropriate algorithm, tuning proposal distributions, choosing burn-in length, and reparameterizing models to improve mixing. With these tools, MCMC enables Bayesian inference for virtually any probabilistic model, from simple conjugate settings to complex hierarchical and nonparametric frameworks.