Gibbs Sampling

Abstract: Gibbs sampling is a Markov Chain Monte Carlo (MCMC) algorithm used to obtain a sequence of observations that approximates a specified multivariate probability distribution when direct sampling is difficult. By iteratively sampling from the conditional distribution of each variable given the current values of all others, the algorithm constructs a Markov chain whose stationary distribution matches the target joint distribution. It is widely used in Bayesian inference, latent variable models, and statistical physics.

Introduction

Gibbs sampling, introduced by Geman and Geman in 1984 for image reconstruction and popularized by Gelfand and Smith in 1990 for Bayesian computation, is a cornerstone algorithm in modern statistical computing. It belongs to the broader class of Markov Chain Monte Carlo methods, which are designed to sample from complex, high-dimensional distributions where analytical integration or direct sampling is intractable.

The core intuition behind Gibbs sampling is elegant: rather than sampling directly from a difficult joint distribution \(p(x_1, x_2, \dots, x_k)\), one can sample from a sequence of simpler conditional distributions \(p(x_i \mid x_{\neq i})\). Under mild regularity conditions, the resulting Markov chain converges to the target joint distribution as its stationary distribution.

Mathematical Formulation

Let \(\mathbf{X} = (X_1, X_2, \dots, X_k)\) be a random vector with joint probability distribution \(\pi(\mathbf{x})\) defined on \(\mathcal{X} \subseteq \mathbb{R}^k\). Assume \(\pi\) is known up to a normalizing constant, and that the full conditional distributions \(\pi(x_i \mid x_1, \dots, x_{i-1}, x_{i+1}, \dots, x_k)\) are available and tractable for sampling.

At iteration \(t\), the Gibbs sampler updates each component sequentially:

\[ X_i^{(t+1)} \sim \pi\left(X_i \mid X_1^{(t+1)}, \dots, X_{i-1}^{(t+1)}, X_{i+1}^{(t)}, \dots, X_k^{(t)}\right) \]

Notice that the most recently updated components are used immediately in subsequent conditional draws within the same iteration. This deterministic scan ensures the Markov property and ergodicity under standard conditions (irreducibility and aperiodicity).

Algorithm

The standard Gibbs sampling procedure can be expressed algorithmically as follows:

// Initialize all variables for i in 1..k: X[i]random_sample() // Iterate for N steps for t in 1..N: for i in 1..k: X[i]sample(X[i] | X[-i]) yield X // Store or analyze current state

Key implementation notes:

  • The order of variable updates can be deterministic (as shown) or randomized (random scan Gibbs), with minimal impact on convergence in practice.
  • Block Gibbs sampling updates subsets of variables jointly when conditionals are more tractable in groups.
  • Collapsed Gibbs sampling integrates out certain variables analytically before sampling, reducing autocorrelation.

Convergence & Diagnostics

Like all MCMC methods, Gibbs sampling produces dependent samples. Convergence to the stationary distribution is asymptotic, and practical implementations require:

  1. Burn-in period: Discarding initial samples before the chain reaches stationarity.
  2. Thinning: Retaining every \(k\)-th sample to reduce autocorrelation (though modern practice often prefers storing all samples and accounting for dependence analytically).
  3. Diagnostics: Trace plots, autocorrelation plots, Gelman-Rubin statistic (for multiple chains), and effective sample size (ESS) estimation.
⚠️ Practical Warning High correlation between variables can cause Gibbs sampling to mix poorly, resulting in long runtimes or false convergence. In such cases, reparameterization, Hamiltonian Monte Carlo, or Metropolis-adjusted updates are often preferred.

Applications

Gibbs sampling is ubiquitous across quantitative disciplines:

  • Bayesian Inference: Posterior sampling when conjugate priels yield tractable conditionals.
  • Latent Variable Models: Latent Dirichlet Allocation (LDA), Hidden Markov Models, Gaussian Mixture Models.
  • Statistical Physics: Sampling configurations in Ising models and Potts models for phase transition analysis.
  • Missing Data Imputation: Iteratively sampling missing values given observed data and model parameters.
  • Computer Vision: Image segmentation, stereo matching, and restoration via conditional pixel sampling.

Comparison with Other MCMC Methods

Method Accept/Reject Step Best Suited For
Gibbs Sampling None (100% acceptance) High-dimensional spaces with known conditionals
Metropolis-Hastings Yes When conditionals are intractable or proposal distributions are tuned
Hamiltonian Monte Carlo Yes (gradient-guided) Strongly correlated posteriors, continuous parameters

References

  1. [1] Geman, S., & Geman, D. (1984). Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images. IEEE Transactions on Pattern Analysis and Machine Intelligence, 6(6), 721-741.
  2. [2] Gelfand, A. E., & Smith, A. F. M. (1990). Sampling-based approaches to calculating marginal densities. Journal of the American Statistical Association, 85(410), 398-409.
  3. [3] Robert, C. P., & Casella, G. (2004). Monte Carlo Statistical Methods (2nd ed.). Springer.
  4. [4] Bishop, C. M. (2006). Pattern Recognition and Machine Learning. Springer. (See Section 11.4 for LDA implementation)