Gibbs & Metropolis-Hastings Sampling — Gibbs sampling and...

SECTION 1

Historical Context & Motivation

Bayesian inference requires computing the posterior distribution p(θ | y) ∝ p(y | θ) p(θ), but in most realistic models the normalizing constant — the marginal likelihood p(y) — involves an intractable integral over a high-dimensional parameter space. Before the advent of Markov chain Monte Carlo (MCMC) methods, practitioners were largely confined to conjugate models or crude numerical quadrature, both of which scale poorly with dimensionality. The development of MCMC fundamentally changed statistical practice by offering a general-purpose computational strategy: construct a Markov chain whose stationary distribution equals the target posterior, run the chain long enough, and treat the resulting samples as draws from the posterior. Two of the most influential MCMC algorithms — the Metropolis–Hastings algorithm and Gibbs sampling — emerged from different scientific communities and converged into the cornerstone toolkit of modern Bayesian computation.

1953

The Metropolis Algorithm

Nicholas Metropolis, Arianna and Marshall Rosenbluth, and Augusta and Edward Teller published their landmark algorithm for simulating equations of state via Monte Carlo sampling on the MANIAC computer at Los Alamos. The method introduced the accept-reject mechanism based on energy ratios, enabling simulation from the Boltzmann distribution.

1970

Hastings Generalization

W. Keith Hastings extended the Metropolis algorithm to allow asymmetric proposal distributions, producing the general Metropolis–Hastings framework. This broadened the class of problems amenable to MCMC by introducing a correction ratio that accounts for asymmetry in the proposal kernel.

1984

Geman & Geman — Gibbs Sampling

Stuart and Donald Geman introduced the Gibbs sampler in the context of image restoration using Markov random fields. Named after the physicist Josiah Willard Gibbs, the algorithm systematically samples each variable from its full conditional distribution, avoiding the need for a tunable proposal distribution.

1990

Gelfand & Smith — Bayesian Adoption

Alan Gelfand and Adrian Smith demonstrated the power of the Gibbs sampler for mainstream Bayesian inference, catalyzing a revolution in applied statistics. Their paper made it clear that hierarchical models, previously viewed as theoretically elegant but computationally hopeless, could be fit routinely.

2000s+

Modern Variants & Software

Software packages such as BUGS, JAGS, and Stan brought MCMC to a wide audience. Advanced variants — Hamiltonian Monte Carlo, No-U-Turn Sampler, reversible-jump MCMC — extended the reach of MCMC to ever more complex models.

The central question that both the Metropolis–Hastings algorithm and the Gibbs sampler answer is deceptively simple: how can we generate samples from a probability distribution that we can evaluate (up to a normalizing constant) but cannot sample from directly? Understanding the conceptual foundations of these algorithms — their construction, convergence guarantees, and relative strengths — is essential for any practitioner of Bayesian statistics, whether the goal is fitting a simple hierarchical model or deploying a complex state-space model in a research setting.

SECTION 2

Core Principles & Definitions

Before examining the specific algorithms, it is important to ground ourselves in the conceptual machinery that makes MCMC work. Both the Metropolis–Hastings algorithm and the Gibbs sampler are instances of a broader family of methods that exploit the theory of Markov chains — stochastic processes in which the probability of transitioning to the next state depends only on the current state, not on the history of states visited. The key insight is that if we design a chain with a transition kernel satisfying certain regularity conditions (irreducibility, aperiodicity, and positive recurrence), the chain will possess a unique stationary distribution. By engineering that stationary distribution to equal our target posterior, the chain's trajectory becomes an (asymptotically valid) sample from the posterior.

Markov Property

The next state of the chain depends only on the current state, not the full trajectory. Formally, P(θ^(t+1) | θ^(t), θ^(t−1), …) = P(θ^(t+1) | θ^(t)). This memoryless property simplifies analysis and guarantees convergence results.

Stationary Distribution

A distribution π is stationary for a chain with transition kernel T if π(θ') = ∫ T(θ' | θ) π(θ) dθ. Under irreducibility and aperiodicity, this stationary distribution is unique. MCMC designs T so that π = p(θ | y), the posterior.

Detailed Balance

A sufficient (though not necessary) condition for stationarity is detailed balance: π(θ) T(θ' | θ) = π(θ') T(θ | θ'). This reversibility condition ensures that the net probability flow between any two states is zero, guaranteeing π as the stationary distribution.

Ergodicity & Convergence

An ergodic chain visits every region of the state space in proportion to the target density, given enough iterations. The ergodic theorem guarantees that time averages of the chain converge to expectations under π, analogous to the law of large numbers for independent samples.

Burn-in & Mixing

Early iterates may be influenced by the starting value and should be discarded (burn-in). Good mixing means the chain traverses the support of the posterior rapidly. Poor mixing — visible as high autocorrelation — necessitates longer runs or algorithmic improvements.

✦ KEY TAKEAWAY

Think of an MCMC sampler as a hiker exploring a mountainous landscape in dense fog: the hiker can evaluate the elevation (unnormalized density) at any point but cannot see the full terrain. By following a carefully designed random-walk rule — sometimes accepting uphill moves, sometimes rejecting them — the hiker spends time at each location in proportion to its elevation. After a long trek, the histogram of visited locations reveals the full topography of the landscape, just as the chain's samples reveal the posterior distribution.

SECTION 3

Visual Explanation — How the Algorithms Explore Parameter Space

The following diagram contrasts the behavior of the Metropolis–Hastings algorithm with that of the Gibbs sampler on a bivariate posterior. Metropolis–Hastings proposes moves in arbitrary directions, accepting or rejecting each proposal via an acceptance probability. Gibbs sampling, by contrast, cycles through each coordinate in turn, drawing from the full conditional distribution along that axis. The visual difference is striking: Metropolis–Hastings traces a random walk with occasional rejections (shown as stationary dots), while Gibbs sampling moves in axis-aligned steps, always accepting.

Left: Metropolis–Hastings traces a random walk; red dots indicate rejected proposals where the chain stays put. Right: Gibbs sampling moves in axis-aligned steps, alternating between θ₁ and θ₂; every draw is accepted because it comes from the exact full conditional distribution.

Observe the characteristic difference in trajectory shape. The Metropolis–Hastings path wanders freely through parameter space, but the red dots — rejected proposals — represent wasted computation. The Gibbs sampler, in the right panel, takes orthogonal steps: each iteration updates one coordinate while holding the others fixed. This structure eliminates rejections entirely but introduces strong autocorrelation when the posterior exhibits high correlation between parameters, since axis-aligned moves struggle to traverse the long axis of an elongated elliptical contour. This visual contrast encapsulates the fundamental trade-off between the two algorithms: Metropolis–Hastings offers flexibility in proposal design at the cost of potential rejections, while Gibbs sampling guarantees acceptance at the cost of requiring tractable full conditionals and potentially slow mixing under strong correlations.

SECTION 4

Mathematical Framework

The Metropolis–Hastings Algorithm

Suppose we wish to sample from a target density π(θ) that we can evaluate up to a normalizing constant — that is, we know π(θ) ∝ f(θ) for some computable function f. The Metropolis–Hastings algorithm proceeds by iterating two steps at each time t. First, a candidate θ* is drawn from a proposal distribution q(θ* | θ^(t)). Second, the candidate is accepted with a probability α that corrects for any asymmetry in the proposal distribution and for the ratio of target densities.

METROPOLIS–HASTINGS ACCEPTANCE PROBABILITY

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

π(θ) = target density (known up to a constant); q(θ* | θ^(t)) = proposal density; the normalizing constant of π cancels in the ratio. When q is symmetric, q(θ* | θ) = q(θ | θ*), the expression reduces to the original Metropolis ratio α = min{1, π(θ*)/π(θ^(t))}.

After computing α, we draw u ~ Uniform(0, 1). If u ≤ α, we set θ^(t+1) = θ*; otherwise θ^(t+1) = θ^(t). This accept-reject mechanism ensures that the chain satisfies detailed balance with respect to π, and therefore π is the stationary distribution of the chain. The beauty of the construction is that the unknown normalizing constant cancels in the ratio, so we never need to compute p(y).

The Gibbs Sampler

Now consider a parameter vector θ = (θ₁, θ₂, …, θ_d). The Gibbs sampler replaces the joint sampling problem with a sequence of univariate (or lower-dimensional) sampling problems. At each iteration, the algorithm cycles through the d components, drawing each from its full conditional distribution — the distribution of that component given the current values of all others.

GIBBS SAMPLING UPDATE

θⱼ⁽ᵗ⁺¹⁾ ~ p(θⱼ | θ₁⁽ᵗ⁺¹⁾, …, θⱼ₋₁⁽ᵗ⁺¹⁾, θⱼ₊₁⁽ᵗ⁾, …, θ_d⁽ᵗ⁾, y) for j = 1, 2, …, d

Each component θⱼ is drawn from its full conditional distribution, conditioning on the most recently updated values of all preceding components (superscript t+1) and the current values of all subsequent components (superscript t). This deterministic scan strategy is the most common, though random scan variants also exist.

GIBBS AS A SPECIAL CASE OF MH

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

When the proposal distribution for component j is the full conditional p(θⱼ | θ₋ⱼ), one can show that the Metropolis–Hastings acceptance probability equals 1 for every proposal. Thus Gibbs sampling is a special case of Metropolis–Hastings with no rejections. The proof follows from the identity π(θ*) / π(θ) = p(θⱼ* | θ₋ⱼ) / p(θⱼ | θ₋ⱼ) and q(θ | θ*) / q(θ* | θ) = p(θⱼ | θ₋ⱼ) / p(θⱼ* | θ₋ⱼ), whose product is 1.

💡 Why Full Conditionals?

In many Bayesian models — particularly conjugate or conditionally conjugate hierarchical models — the full conditional distributions take recognizable parametric forms (normal, gamma, beta, etc.), making them easy to sample from. This is why the Gibbs sampler became the workhorse of Bayesian software such as BUGS and JAGS: it automates the decomposition of a complex joint posterior into a cycle of simple univariate draws.

SECTION 5

Detailed Comparison — MH vs. Gibbs

While the Gibbs sampler can be viewed as a special case of Metropolis–Hastings, the practical considerations governing their use are quite different. The diagram below illustrates a decision flowchart for choosing between the two algorithms, and the accompanying table summarizes their key properties.

Decision flowchart for selecting between Gibbs sampling and Metropolis–Hastings. When full conditionals are available in closed form, Gibbs is typically preferred for its simplicity and zero rejection rate. When full conditionals are intractable, Metropolis–Hastings offers a flexible alternative with tunable proposal distributions.

Comparison of Metropolis–Hastings and Gibbs Sampling

Feature	Metropolis–Hastings	Gibbs Sampling
Proposal mechanism	User-specified q(θ* \| θ); can be any distribution on the parameter space	Full conditional distributions p(θⱼ \| θ₋ⱼ, y); must be available in tractable form
Acceptance rate	Varies with q; optimal ≈ 23.4% for random walk MH in high dimensions (Roberts et al., 1997)	100% — every draw is accepted because the proposal equals the full conditional
Tuning required	Yes — proposal scale/shape must be tuned for good acceptance and mixing	Minimal — no proposal tuning, but may need blocking or reparameterization
Dimensionality	Can propose moves in all dimensions simultaneously; struggles in very high dimensions without gradient information	Updates one (or a block of) component(s) per step; scales well if conditionals are easy
Correlated posteriors	Can use non-axis-aligned proposals to navigate correlation structure	Axis-aligned moves cause slow mixing under strong posterior correlations

SECTION 6

Worked Example — Gibbs Sampler for a Normal Model

Consider the classical conjugate normal model. We observe y₁, y₂, …, yₙ iid ~ Normal(μ, σ²), and we wish to perform Bayesian inference on both the mean μ and the variance σ² (or equivalently the precision τ = 1/σ²). We place standard conjugate priors: μ ~ Normal(μ₀, σ₀²) and τ ~ Gamma(a, b). The joint posterior p(μ, τ | y) does not simplify to a single recognizable distribution in closed form, but the full conditionals do. This makes the problem ideal for the Gibbs sampler.

Gibbs Sampler for Normal(μ, τ⁻¹) with Conjugate Priors

Step 1 — Derive the Full Conditional for μConditioning on τ (and data y), the posterior for μ is proportional to the product of the likelihood and the prior on μ. Both are normal kernels in μ, so the posterior is normal. After completing the square, we find: p(μ | τ, y) = Normal(μ_n, σ_n²), where σ_n² = (nτ + 1/σ₀²)⁻¹ and μ_n = σ_n² × (nτȳ + μ₀/σ₀²). Here ȳ is the sample mean.

μ | τ, y ~ Normal(μ_n, σ_n²)

Step 2 — Derive the Full Conditional for τConditioning on μ (and data y), the posterior for τ is proportional to the product of the normal likelihood (viewed as a function of τ) and the Gamma prior on τ. This yields a Gamma distribution. Specifically: p(τ | μ, y) = Gamma(a_n, b_n), where a_n = a + n/2 and b_n = b + ½ Σᵢ (yᵢ − μ)².

τ | μ, y ~ Gamma(a + n/2, b + ½ Σ(yᵢ − μ)²)

Step 3 — Initialize the ChainChoose starting values μ⁽⁰⁾ and τ⁽⁰⁾. Common choices include the sample mean and sample precision, or draws from the prior. The chain's long-run behavior is insensitive to starting values, but good initialization reduces the required burn-in period.

Set μ⁽⁰⁾ = ȳ, τ⁽⁰⁾ = 1/s²

Step 4 — Iterate the Gibbs CycleFor t = 0, 1, 2, …, T: (a) draw μ⁽ᵗ⁺¹⁾ from Normal(μ_n, σ_n²) using the current τ⁽ᵗ⁾; then (b) draw τ⁽ᵗ⁺¹⁾ from Gamma(a_n, b_n) using the just-sampled μ⁽ᵗ⁺¹⁾. Notice that step (b) uses the freshly updated μ, not the old value — this is the systematic (deterministic) scan Gibbs sampler.

After discarding burn-in, {μ⁽ᵗ⁾, τ⁽ᵗ⁾} are samples from the joint posterior p(μ, τ | y).

Step 5 — Numerical IllustrationSuppose n = 20, ȳ = 5.3, s² = 2.1, and priors μ₀ = 0, σ₀² = 100, a = 0.01, b = 0.01. At iteration t with τ⁽ᵗ⁾ = 0.5: σ_n² = (20 × 0.5 + 1/100)⁻¹ = (10.01)⁻¹ ≈ 0.0999 and μ_n = 0.0999 × (20 × 0.5 × 5.3 + 0/100) = 0.0999 × 53 ≈ 5.295. So we draw μ⁽ᵗ⁺¹⁾ ~ Normal(5.295, 0.0999). Then we compute b_n = 0.01 + ½ × Σ(yᵢ − μ⁽ᵗ⁺¹⁾)² and draw τ⁽ᵗ⁺¹⁾ from the resulting Gamma.

μ⁽ᵗ⁺¹⁾ ~ Normal(5.295, 0.0999); τ⁽ᵗ⁺¹⁾ ~ Gamma(10.01, b_n)

SECTION 7

Strengths, Limitations, and Practical Considerations

Practical Strengths and Limitations

Consideration	Strength	Limitation
Ease of implementation	Gibbs is simple when full conditionals are known; MH requires only the ability to evaluate the target up to a constant.	MH needs careful tuning of the proposal; Gibbs requires derivation of every full conditional, which may be analytically demanding.
Mixing speed	MH can use non-axis-aligned or gradient-informed proposals; Gibbs can use blocking to update correlated parameters jointly.	Random walk MH mixes slowly in high dimensions; standard Gibbs suffers under strong posterior correlations.
Generality	MH applies to virtually any target density; it is the most general MCMC framework.	Gibbs is restricted to models where all full conditionals are available in tractable form.
Convergence diagnosis	Standard diagnostics (R-hat, trace plots, ESS) apply to both algorithms equally.	Slow mixing can masquerade as convergence; multimodal targets are challenging for both methods without specialized extensions.
Scalability	Gibbs naturally decomposes the problem; data-augmentation variants (e.g., EM-like) are powerful for latent variable models.	Both methods can be prohibitively slow for very large datasets without sub-sampling or variational approximation.

✦ KEY TAKEAWAY

Think of Gibbs sampling as navigating a city on a strict grid system — you can only move north-south or east-west. This works beautifully when your destination lies roughly along the grid axes, but if you need to reach a point on a diagonal, the zig-zagging path is inefficient. Metropolis–Hastings, by contrast, lets you walk in any direction, but you must occasionally check a compass (the acceptance probability) and backtrack if you have wandered too far. Neither method is universally superior; the choice depends on the geometry of the posterior landscape and the availability of analytical full conditionals.

SECTION 8

Connection to Advanced MCMC Methods

The Metropolis–Hastings algorithm and the Gibbs sampler laid the conceptual groundwork for a rich ecosystem of advanced MCMC methods. Understanding where these foundational algorithms sit in the broader landscape helps contextualize their role in modern Bayesian practice. The table below highlights the progression from classical MCMC to its modern descendants.

From Classical to Modern MCMC

Classical MCMC	Modern Extension	Key Innovation
Random Walk Metropolis–Hastings	Hamiltonian Monte Carlo (HMC)	Uses gradient information (Hamiltonian dynamics) to propose distant, high-probability states; dramatically reduces random-walk behavior.
HMC with manual tuning	No-U-Turn Sampler (NUTS)	Automatically tunes the trajectory length in HMC by detecting when the trajectory begins to double back; implemented in Stan.
Standard Gibbs Sampler	Slice Sampling	Introduces an auxiliary variable to define a 'slice' under the density; avoids the need for a closed-form full conditional while maintaining automatic acceptance.
Fixed-dimension MH	Reversible-Jump MCMC	Extends MH to spaces of varying dimension, enabling Bayesian model selection (e.g., choosing the number of mixture components).
Single-chain MCMC	Parallel Tempering	Runs multiple chains at different 'temperatures'; hot chains explore broadly and swap states with cold chains, enabling escape from local modes.

Despite the proliferation of more sophisticated samplers, the Gibbs sampler and Metropolis–Hastings remain indispensable. Many modern algorithms use them as sub-routines — for instance, a Metropolis-within-Gibbs strategy applies MH steps for components whose full conditionals are not available in closed form, while using standard Gibbs updates for the remaining components. This hybrid approach exemplifies the modular nature of MCMC: once you understand the two foundational algorithms, composing and extending them is largely a matter of engineering. Furthermore, the conceptual clarity of MH and Gibbs provides the theoretical lens through which we analyze more advanced methods — their convergence properties, detailed balance conditions, and mixing diagnostics all trace back to the principles introduced in this lesson.

SECTION 9

Practice Problems

PROBLEM 1 — CONCEPTUAL

Explain why the normalizing constant of the posterior distribution (i.e., the marginal likelihood p(y)) is irrelevant when implementing the Metropolis–Hastings algorithm. What specific feature of the MH acceptance ratio makes this cancellation possible?

PROBLEM 2 — BASIC CALCULATION

Suppose we are running a Metropolis–Hastings algorithm with a symmetric proposal distribution to sample from a target π(θ) ∝ exp(−θ²/2). The current state is θ⁽ᵗ⁾ = 1.0, and the proposed state is θ* = 2.5. Compute the acceptance probability α. Should the proposal be accepted if the drawn uniform random variable is u = 0.15?

PROBLEM 3 — INTERMEDIATE

Consider a bivariate normal posterior p(μ, τ | y) where μ | τ, y ~ Normal(3.0, (10τ)⁻¹) and τ | μ, y ~ Gamma(6, β(μ)) with β(μ) = 2 + 5(μ − 3)². Starting from (μ⁽⁰⁾, τ⁽⁰⁾) = (0, 1), write out the first complete iteration of the Gibbs sampler — specifying the exact distributions from which each draw is made — and explain why these full conditionals guarantee convergence to the joint posterior.

PROBLEM 4 — APPLIED

A researcher fits a Bayesian hierarchical model using a Gibbs sampler with 10,000 iterations (after discarding 2,000 as burn-in). The trace plot for one parameter shows clear trending behavior (the chain slowly drifts upward through the entire run), and the effective sample size (ESS) is only 35. Diagnose the problem, explain its implications for posterior inference, and propose at least two concrete remedies.

PROBLEM 5 — CRITICAL THINKING

Prove that the Gibbs sampler is a special case of the Metropolis–Hastings algorithm by showing that the MH acceptance probability equals 1 when the proposal distribution for component θⱼ is the full conditional p(θⱼ | θ₋ⱼ). Discuss the implications of this result for the design of hybrid MCMC algorithms that combine Gibbs and MH steps.

SUMMARY

Lesson Summary

This lesson introduced the two foundational Markov chain Monte Carlo algorithms in Bayesian statistics. The Metropolis–Hastings algorithm is the most general MCMC framework: it proposes candidate states from a user-chosen proposal distribution and accepts or rejects them via an acceptance probability that corrects for any mismatch between the proposal and the target posterior. The key identity — that the unknown normalizing constant cancels in the acceptance ratio — is what makes Bayesian computation feasible for complex models. The Gibbs sampler is a special case of MH that decomposes the joint sampling problem into a cycle of draws from full conditional distributions, achieving a 100% acceptance rate at every step. Both algorithms rely on the theory of Markov chains — specifically detailed balance and ergodicity — to guarantee convergence to the target posterior.

In practice, Gibbs sampling excels in conditionally conjugate models where full conditionals are recognizable distributions, while Metropolis–Hastings serves as the fallback for any parameter whose full conditional is intractable. The Metropolis-within-Gibbs hybrid combines the best of both worlds. Challenges such as slow mixing under strong posterior correlations and convergence diagnostics remain central practical concerns, and modern successors like Hamiltonian Monte Carlo address these by incorporating gradient information. Mastery of Gibbs and MH provides the conceptual foundation needed to understand, implement, and critique the entire MCMC toolkit.