Gibbs Sampling

Gibbs Sampling

Gibbs sampling is an MCMC algorithm for a joint target $\pi ({\theta }_1,\dots ,{\theta }_d)$ in which each coordinate is updated by drawing from its full conditional $\pi ({\theta }_j\mid {\theta }_{-j})$, where ${\theta }_{-j}$ denotes all coordinates other than ${\theta }_j$. It applies whenever the full conditionals are available in closed form (or at least easy to sample from) even though the joint posterior is not. The canonical systematic-scan sweep is:

Given the current state ${\theta }^{(t)}=({\theta }_1^{(t)},\dots ,{\theta }_d^{(t)})$, produce ${\theta }^{(t+1)}$ by:

1. Draw ${\theta }_1^{(t+1)}\sim \pi \left({\theta }_1\right|{\theta }_2^{(t)},{\theta }_3^{(t)},\dots ,{\theta }_d^{(t)})$.
2. Draw ${\theta }_2^{(t+1)}\sim \pi \left({\theta }_2\right|{\theta }_1^{(t+1)},{\theta }_3^{(t)},\dots ,{\theta }_d^{(t)})$.
3. $⋮$
4. Draw ${\theta }_d^{(t+1)}\sim \pi \left({\theta }_d\right|{\theta }_1^{(t+1)},\dots ,{\theta }_{d-1}^{(t+1)})$.

Each coordinate update uses the most recently drawn values of every other coordinate.

Gibbs is the special case of Metropolis-Hastings in which the proposal for coordinate $j$ is its own full conditional and every proposal is accepted. The MH acceptance ratio reduces to $1$ identically:
$\alpha =\min \left(1,\frac{\pi ({\theta }_j^{\prime },{\theta }_{-j})\pi ({\theta }_j\mid {\theta }_{-j})}{\pi ({\theta }_j,{\theta }_{-j})\pi ({\theta }_j^{\prime }\mid {\theta }_{-j})}\right)=1,$
since $\pi ({\theta }_j,{\theta }_{-j})=\pi ({\theta }_{-j})\pi ({\theta }_j\mid {\theta }_{-j})$ and the $\pi ({\theta }_{-j})$ factors cancel. Each coordinate update therefore satisfies detailed balance with respect to $\pi$, and a full sweep preserves $\pi$ as the stationary distribution. Under irreducibility and aperiodicity, the ergodic theorem gives
$\frac{1}{N}{\sum }_{t=1}^{N}g({\theta }^{(t)})\stackrel{\text{a.s.}}{\to }\int g(\theta )\pi (\theta )d\theta .$

Random-scan Gibbs picks a coordinate $j\in \{1,\dots ,d\}$ uniformly at random at each step and updates only that coordinate; systematic-scan (above) cycles through all coordinates in a fixed order. Random-scan is reversible with respect to $\pi$; systematic-scan is not, but still admits $\pi$ as its stationary distribution.

Data augmentation (Tanner-Wong, 1987) extends Gibbs to models whose raw full conditionals are intractable by introducing latent variables $z$ such that $\pi (\theta ,z\mid y)$ has tractable conditionals $\pi (\theta \mid z,y)$ and $\pi (z\mid \theta ,y)$. Running Gibbs on $(\theta ,z)$ and discarding the $z$ component yields samples from the marginal $\pi (\theta \mid y)$. This is the standard route for probit regression (with Gaussian latent variables) and mixture models (with categorical class indicators).

Metropolis-within-Gibbs handles partial conjugacy: coordinates with closed-form conditionals are updated by Gibbs draws; coordinates without are updated by a Metropolis-Hastings step against the full conditional as target. The combined sweep still preserves $\pi$.

Intuition

The joint target $\pi$ is often intractable, but “condition on everything except ${\theta }_j$” frequently leaves a distribution whose kernel is a familiar family (normal, gamma, beta). Each Gibbs step exploits conditional tractability — the full conditional is a univariate problem, and univariate problems are tractable far more often than their joint counterpart. Cycling through the coordinates traces a path through the joint that keeps at least one axis moving at any given step; over many sweeps the path fills out the joint in proportion to $\pi$.

The central tension is coordinate correlation. If two coordinates ${\theta }_i,{\theta }_j$ are nearly collinear under $\pi$, conditioning on one pins the other in a narrow band and each update moves the state only along a thin ridge, producing strong autocorrelation and slow mixing. Reparameterizing to less-correlated coordinates, updating blocks of correlated coordinates jointly, or switching to a gradient-aware sampler (MALA, HMC) are the standard remedies. Gibbs excels at conditionally-conjugate hierarchies and at problems where latent-variable augmentation collapses a hard joint into easy conditionals.

Properties

(Always Accept) The Gibbs coordinate update is an MH step with acceptance probability identically $1$.

Full-conditional proposals match the conditional density of the target, so the Hastings ratio $\pi ({\theta }_j^{\prime },{\theta }_{-j})q({\theta }_j\mid {\theta }^{\prime })/\pi ({\theta }_j,{\theta }_{-j})q({\theta }_j^{\prime }\mid \theta )$ reduces via $\pi ({\theta }_j,{\theta }_{-j})=\pi ({\theta }_{-j})\pi ({\theta }_j\mid {\theta }_{-j})$ to $1$ exactly. This is the algebraic identity that makes Gibbs tune-free: there is no step size to pick, no acceptance rate to chase, and no proposal-scale trade-off. The cost is that the full conditional must be sampleable; when it is not, one falls back to MH.

(Coordinate Update Preserves $\pi$) Each single-coordinate update satisfies detailed balance; composing updates preserves the invariant distribution.

For a single-coordinate update at index $j$, the transition kernel is $K_j(\theta ,{\theta }^{\prime })=\pi ({\theta }_j^{\prime }\mid {\theta }_{-j})1\{{\theta }_{-j}^{\prime }={\theta }_{-j}\}$. Detailed balance follows because
$\pi (\theta )K_j(\theta ,{\theta }^{\prime })=\pi ({\theta }_{-j})\pi ({\theta }_j\mid {\theta }_{-j})\pi ({\theta }_j^{\prime }\mid {\theta }_{-j})=\pi ({\theta }_{-j}^{\prime })\pi ({\theta }_j^{\prime }\mid {\theta }_{-j}^{\prime })\pi ({\theta }_j\mid {\theta }_{-j}^{\prime })=\pi ({\theta }^{\prime })K_j({\theta }^{\prime },\theta ),$
using ${\theta }_{-j}^{\prime }={\theta }_{-j}$. A full systematic-scan sweep is the composition $K_d\circ K_{d-1}\circ \cdots \circ K_1$; composition of $\pi$-invariant kernels is $\pi$-invariant, so the sweep preserves $\pi$. Each individual $K_j$ is reversible; the composition is reversible if the scan order is randomized (random-scan Gibbs) but not in general for systematic scan.

(Data Augmentation) Introducing a latent variable $z$ can turn an intractable joint into a tractable pair of conditionals.

When the direct posterior $\pi (\theta \mid y)$ has no closed-form conditional, one can often introduce a latent $z$ such that the augmented posterior $\pi (\theta ,z\mid y)$ satisfies
$\pi (\theta \mid z,y)\text{ and }\pi (z\mid \theta ,y)\text{ both tractable}.$
Running Gibbs on $(\theta ,z)$ and discarding $z$ produces samples from $\pi (\theta \mid y)$ by marginalization. Canonical instances are probit regression ($z_i=x_i^{\top }\beta +{\epsilon }_i$ with ${\epsilon }_i\sim N(0,1)$ and $y_i=1\{z_i>0\}$, due to Albert-Chib, 1993) and mixture models ($z_i\in \{1,\dots ,K\}$ indicating the component of observation $i$, due to Diebolt-Robert, 1994). The same pattern covers missing data, censored data, and hierarchical models.

(Metropolis-within-Gibbs) A coordinate without a closed-form conditional can be updated by a Metropolis step against the conditional as target.

If $\pi ({\theta }_j\mid {\theta }_{-j})$ is known up to a normalizing constant but not sampleable, a single MH step with proposal $q({\theta }_j^{\prime }\mid {\theta }_j)$ and acceptance probability
$\alpha =\min \left(1,\frac{\pi ({\theta }_j^{\prime }\mid {\theta }_{-j})q({\theta }_j\mid {\theta }_j^{\prime })}{\pi ({\theta }_j\mid {\theta }_{-j})q({\theta }_j^{\prime }\mid {\theta }_j)}\right)$
preserves $\pi$ as the stationary distribution. Because only the coordinate-$j$ slice of $\pi$ matters, the $\pi ({\theta }_{-j})$ factor cancels and only the unnormalized conditional is needed. Hybrid sweeps — Gibbs for conjugate coordinates, MH (or slice sampling) for non-conjugate coordinates — are the standard way to handle partially conjugate models.

(Mixing Failure Under High Coordinate Correlation) Strong correlation between coordinates produces thin conditional ridges and slow exploration.

If $\theta$ is approximately bivariate normal with correlation $\rho$ close to $\pm 1$, the conditional $\pi ({\theta }_1\mid {\theta }_2)$ has variance ${\sigma }_1^2(1-{\rho }^2)$, narrow when $|\rho |$ is large. A Gibbs sweep moves ${\theta }_1$ within this narrow band, then ${\theta }_2$ within a symmetrically narrow band, producing a sequence of small orthogonal steps up a diagonal ridge. The chain eventually fills the joint but at an autocorrelation timescale that grows as $(1-{\rho }^2{)}^{-1}$. Standard remedies: block Gibbs (update correlated coordinates jointly from the multivariate conditional), reparameterize to uncorrelated coordinates (e.g., center-scale a hierarchical model to reduce prior-likelihood confounding), or switch to HMC which uses gradient information to step along ridge directions rather than axis-aligned.

(Systematic vs. Random Scan) Random-scan Gibbs is reversible; systematic-scan is generally not, but both have $\pi$ stationary.

Systematic scan cycles through coordinates in a fixed order $1,2,\dots ,d$; random scan picks a coordinate uniformly at each step and updates only that one. Systematic scan typically mixes faster because it guarantees every coordinate is updated each sweep, but its composite transition $K_d\circ \cdots \circ K_1$ is not symmetric under reversal. Random scan is reversible (the coordinate index and update are symmetric in time), which simplifies some theoretical arguments but wastes effort by revisiting the same coordinate before updating others. In practice, systematic scan is the default.

(Irreducibility Requirement) Gibbs can fail catastrophically when the support of $\pi$ is disconnected relative to coordinate-aligned moves.

The chain is irreducible only if coordinate-aligned moves can reach every region of positive $\pi$-mass. A target supported on two disjoint islands aligned diagonally (e.g., $\pi$ concentrated on $\{{\theta }_1={\theta }_2\}\cup \{{\theta }_1=-{\theta }_2\}$ with a gap between them) has no Gibbs path between the islands because changing one coordinate at a time cannot cross the gap. The chain gets stuck on the initial island forever, yielding samples that look like draws from a conditional restriction rather than $\pi$. Multi-modal posteriors (especially in mixture models with label-switching behavior) are the practical version of this trap; running multiple chains from overdispersed starting points is the standard diagnostic.

Examples

Bivariate normal

Let $({\theta }_1,{\theta }_2)$ have a mean-zero bivariate normal target with variances $1$ and correlation $\rho$:
$\pi ({\theta }_1,{\theta }_2)\propto \exp \left(-\frac{1}{2(1-{\rho }^2)}({\theta }_1^2-2\rho {\theta }_1{\theta }_2+{\theta }_2^2)\right).$
The full conditionals are
${\theta }_1\mid {\theta }_2\sim N(\rho {\theta }_2,1-{\rho }^2),{\theta }_2\mid {\theta }_1\sim N(\rho {\theta }_1,1-{\rho }^2),$
read off by completing the square in each coordinate. One Gibbs sweep: draw ${\theta }_1\sim N(\rho {\theta }_2,1-{\rho }^2)$, then draw ${\theta }_2\sim N(\rho {\theta }_1,1-{\rho }^2)$.

With $\rho =0$ the coordinates are independent and a single sweep produces an exact sample. With $\rho =0.99$ the conditional variance $1-{\rho }^2=0.0199$ is tiny and each update nudges the state along the main diagonal by $O(\sqrt{0.0199})\approx 0.14$; traversing the full support requires $O(1/(1-{\rho }^2))=100$ sweeps. A trace plot of ${\theta }_1$ in the high-$\rho$ case shows slow sinusoidal drift, the classic signature of coordinate-correlation-induced slow mixing.

Bayesian linear regression with a conditionally conjugate prior

Consider $y=X\beta +\epsilon$ with $\epsilon \sim N(0,{\sigma }^2{I}_n)$, and priors $\beta \sim N(0,{\tau }^2{I}_p)$ and ${\sigma }^2\sim \text{InvGamma}(a,b)$. The joint posterior is not a recognizable distribution, but each full conditional is conjugate:

For $\beta$, given ${\sigma }^2$:
$\beta \mid {\sigma }^2,y\sim N\left({\mu }_n,V_n\right),V_n={\left(\frac{1}{{\sigma }^2}{X}^{\top }X+\frac{1}{{\tau }^2}{I}_p\right)}^{-1},{\mu }_n=\frac{1}{{\sigma }^2}{V}_n{X}^{\top }y.$

For ${\sigma }^2$, given $\beta$:
${\sigma }^2\mid \beta ,y\sim \text{InvGamma}\left(a+\frac{n}{2},b+\frac{1}{2}(y-X\beta {)}^{\top }(y-X\beta )\right).$

Both conditionals are standard and easy to sample. A Gibbs sweep alternates the two updates. Initializing at ${\hat{\beta }}_{\text{OLS}}$ and ${\sigma }^2=\Vert y-X{\hat{\beta }}_{\text{OLS}}{\Vert }^2/(n-p)$ places the chain near the posterior mode and minimizes burn-in. This is the prototypical “conditionally conjugate hierarchy” Gibbs target — every hyperparameter in the vault’s hierarchical-models note (3.6 Hierarchical Models) is updated the same way.

Probit regression via data augmentation (Albert-Chib)

Let $y_i\in \{0,1\}$ and model $P(y_i=1\mid x_i)=\Phi (x_i^{\top }\beta )$ with a flat prior on $\beta \in {\mathbb{R}}^p$. The direct likelihood ${\prod }_i\Phi (x_i^{\top }\beta {)}^{y_i}(1-\Phi (x_i^{\top }\beta ){)}^{1-y_i}$ is log-concave in $\beta$ but has no conjugate conditional. Introduce latent $z_i=x_i^{\top }\beta +{\epsilon }_i$ with ${\epsilon }_i\sim N(0,1)$ and $y_i=1\{z_i>0\}$; marginalizing ${\epsilon }_i$ gives back the probit likelihood, so the augmented posterior has the same $\beta$-marginal as the original.

The augmented conditionals are tractable:

$\pi (\beta \mid z,y)\propto \pi (\beta ){\prod }_i\varphi (z_i-x_i^{\top }\beta )$, which under the flat prior gives
$\beta \mid z,y\sim N\left((X^{\top }X{)}^{-1}{X}^{\top }z,(X^{\top }X{)}^{-1}\right),$
a standard Gaussian conjugate update with $z$ playing the role of the “response.”

$\pi (z_i\mid \beta ,y_i)\propto \varphi (z_i-x_i^{\top }\beta )1\{z_i>0⟺y_i=1\}$, a truncated normal:
$z_i\mid \beta ,y_i=1\sim T{N}_{(0,\infty )}(x_i^{\top }\beta ,1),z_i\mid \beta ,y_i=0\sim T{N}_{(-\infty ,0)}(x_i^{\top }\beta ,1).$

A full Gibbs sweep draws every $z_i$ from its truncated normal and then updates $\beta$ from the resulting Gaussian. The $\beta$-marginal of the chain is the probit posterior $\pi (\beta \mid y)$. No tuning, no acceptance rate — the same data-augmentation trick generalizes to ordinal probit (multiple cutoffs) and multinomial probit (Gaussian latent vector).

Gaussian mixture via latent class indicators

Let $y_1,\dots ,y_n$ be drawn from a $K$-component Gaussian mixture
$y_i\sim {\sum }_{k=1}^K{\pi }_{k}N({\mu }_k,{\sigma }_k^2),$
with priors $\pi \sim \text{Dir}({\alpha }_1,\dots ,{\alpha }_K)$, ${\mu }_k\sim N({\mu }_0,{\tau }^2)$, ${\sigma }_k^2\sim \text{InvGamma}(a,b)$. The likelihood is a sum over $K^n$ assignments of data points to components — intractable directly. Introduce latent class indicators $z_i\in \{1,\dots ,K\}$ such that $z_i=k$ means $y_i$ came from component $k$.

The augmented full conditionals are all conjugate:

• $z_i\mid \text{rest}$: categorical with probabilities $P(z_i=k\mid \text{rest})\propto {\pi }_{k}N(y_i;{\mu }_k,{\sigma }_k^2)$.
• $\pi \mid z$: $\text{Dir}({\alpha }_1+n_1,\dots ,{\alpha }_K+n_K)$, where $n_k=\#\{i:z_i=k\}$.
• ${\mu }_k\mid \text{rest}$: Gaussian with mean and variance computed from the observations assigned to component $k$.
• ${\sigma }_k^2\mid \text{rest}$: inverse-gamma with shape and scale updated by the sum of squared residuals in component $k$.

A Gibbs sweep updates $z$, $\pi$, $\mu$, ${\sigma }^2$ in turn. This is the standard finite mixture Gibbs sampler; extending to infinitely many components replaces the Dirichlet prior with a Dirichlet process and adds a Chinese-restaurant-process step for $z$. The deterministic analog of this sampler — which maximizes the augmented posterior rather than sampling from it — is the EM algorithm, with the $z$-update replaced by its expectation (E-step) and the parameter updates replaced by their MAP (M-step).

Derivation

Goal. Show that each Gibbs coordinate update is an MH step with acceptance probability $1$, and that the update kernel preserves $\pi$ as its stationary distribution.

Fix coordinate $j$ and write $\pi (\theta )=\pi ({\theta }_{-j})\pi ({\theta }_j\mid {\theta }_{-j})$. Propose ${\theta }_j^{\prime }\sim \pi ({\theta }_j\mid {\theta }_{-j})$ with ${\theta }_{-j}^{\prime }={\theta }_{-j}$; the proposal kernel is $q({\theta }^{\prime }\mid \theta )=\pi ({\theta }_j^{\prime }\mid {\theta }_{-j})1\{{\theta }_{-j}^{\prime }={\theta }_{-j}\}$. The MH acceptance probability is
$\alpha (\theta ,{\theta }^{\prime })=\min \left(1,\frac{\pi ({\theta }^{\prime })q(\theta \mid {\theta }^{\prime })}{\pi (\theta )q({\theta }^{\prime }\mid \theta )}\right).$
Since ${\theta }_{-j}^{\prime }={\theta }_{-j}$, the indicator factors in $q$ are both $1$ when the proposal is made, and
$\frac{\pi ({\theta }^{\prime })q(\theta \mid {\theta }^{\prime })}{\pi (\theta )q({\theta }^{\prime }\mid \theta )}=\frac{\pi ({\theta }_{-j})\pi ({\theta }_j^{\prime }\mid {\theta }_{-j})\cdot \pi ({\theta }_j\mid {\theta }_{-j})}{\pi ({\theta }_{-j})\pi ({\theta }_j\mid {\theta }_{-j})\cdot \pi ({\theta }_j^{\prime }\mid {\theta }_{-j})}=1.$
So $\alpha \equiv 1$ and every proposal is accepted.

The resulting transition kernel $K_j(\theta ,{\theta }^{\prime })=\pi ({\theta }_j^{\prime }\mid {\theta }_{-j})1\{{\theta }_{-j}^{\prime }={\theta }_{-j}\}$ satisfies detailed balance with $\pi$:
$\pi (\theta )K_j(\theta ,{\theta }^{\prime })=\pi ({\theta }_{-j})\pi ({\theta }_j\mid {\theta }_{-j})\pi ({\theta }_j^{\prime }\mid {\theta }_{-j})=\pi ({\theta }_{-j})\pi ({\theta }_j^{\prime }\mid {\theta }_{-j})\pi ({\theta }_j\mid {\theta }_{-j})=\pi ({\theta }^{\prime })K_j({\theta }^{\prime },\theta ),$
because ${\theta }_{-j}={\theta }_{-j}^{\prime }$ and multiplication is commutative. Detailed balance implies $\pi$ is stationary for $K_j$. The composition $K_d\circ \cdots \circ K_1$ of $\pi$-stationary kernels is $\pi$-stationary, so a full systematic-scan sweep preserves $\pi$. Under irreducibility (every region of positive $\pi$-mass can be reached by coordinate-aligned moves) and aperiodicity, the ergodic theorem delivers a.s. convergence of time averages to posterior expectations.

See Also

• 3.11.1 Metropolis-Hastings — Gibbs is the special case with full-conditional proposals; MH handles cases where the conditional is unavailable.
• 3.11 MCMC — parent index: chain theory, ergodic theorem, diagnostics, tuning.
• 3.11.5 EM Algorithm — deterministic analog that maximizes the augmented posterior instead of sampling from it.
• 3.6 Hierarchical Models — conditionally conjugate hierarchies are the natural application.
• 3.8.3 Bayesian Logistic Regression — non-conjugate sibling that requires Metropolis-within-Gibbs or Laplace-adjusted proposals.
• 5.4 Reversibility — detailed-balance condition that single-coordinate updates satisfy.
• 5.3 Distribution of Markov Chains — stationary distributions and the ergodic theorem.