Gibbs Sampling

A Bayesian posterior over several parameters at once — say, a population mean, a population variance, and a dozen group-level effects — often has no clean formula you can sample from directly. You know the shape of the distribution up to a constant, but you can't just plug it into a random-number generator the way you can for a normal or a beta. Gibbs sampling sidesteps the problem entirely: instead of sampling all the parameters jointly, which is hard, it samples them one at a time, each conditional on the current values of all the others — which, for a huge and useful class of models, turns out to be easy, because each individual conditional distribution is something simple like a normal or a beta. Cycle through all the variables enough times, and the sequence of samples you get behaves exactly as if you'd drawn from the full, complicated joint posterior all along.

An analogy: guessing a room's temperature and humidity by adjusting one dial at a time

Imagine you're trying to find the combination of thermostat and humidifier settings that a room's occupants find comfortable, but you can only test one dial at a time — turn the thermostat while keeping humidity fixed, see how comfort responds, then turn the humidifier while keeping the new temperature fixed, and so on. You never directly search the two-dial space at once. And yet, alternating between the two one-dial adjustments long enough, your sequence of settings wanders around and eventually spends time in each combination roughly in proportion to how comfortable it is — you've explored the full two-dimensional comfort landscape using only one-dimensional moves. Gibbs sampling is this exact strategy applied to probability: sample each parameter's conditional distribution — "given everything else is fixed, what's this one's distribution?" — one at a time, in a loop, and the accumulated samples trace out the full joint posterior.

The algorithm, one symbol at a time

Suppose the posterior is over two parameters, $\theta_1$ and $\theta_2$ ("theta sub 1," "theta sub 2") — the idea extends to any number. Gibbs sampling requires that you can sample from the two full conditional distributions: $p(\theta_1 \mid \theta_2, \text{data})$ and $p(\theta_2 \mid \theta_1, \text{data})$ — "the distribution of $\theta_1$ if you already knew $\theta_2$," and vice versa. Starting from any initial guess $(\theta_1^{(0)}, \theta_2^{(0)})$, at each iteration $t$:

$$\theta_1^{(t+1)} \sim p\big(\theta_1 \mid \theta_2^{(t)},\ \text{data}\big), \qquad \theta_2^{(t+1)} \sim p\big(\theta_2 \mid \theta_1^{(t+1)},\ \text{data}\big) .$$

In words: draw a new value for $\theta_1$ from its conditional distribution, using the current value of $\theta_2$; then draw a new value for $\theta_2$ from its conditional distribution, using the value of $\theta_1$ you just drew — and repeat. Note the second step uses the freshly updated $\theta_1$, not the old one; each new draw is used as soon as it's available. Repeating this for many iterations produces a sequence $(\theta_1^{(1)}, \theta_2^{(1)}), (\theta_1^{(2)}, \theta_2^{(2)}), \dots$ that is a Markov chain: each state depends only on the immediately preceding one. Under mild conditions, this chain's long-run distribution — the distribution of states it visits, once it's run long enough to forget its starting point — converges exactly to the true joint posterior $p(\theta_1, \theta_2 \mid \text{data})$, even though no step ever directly sampled from that joint distribution.

The early samples, taken before the chain has "forgotten" its arbitrary starting point, are discarded as burn-in; the remaining samples are treated as (correlated, but asymptotically valid) draws from the posterior, and any posterior summary — a mean, a credible interval — is computed by simply averaging or ranking across the retained samples.

Worked example 1: a two-step Gibbs sampler for a normal mean and variance

Data: 5 observed daily returns (in percent) $y = [1.2, -0.5, 0.8, 2.1, 0.4]$, sample mean $\bar y = 0.80$. Model: $y_i \sim N(\mu, \sigma^2)$, with a simplified setup where, conditional on $\sigma^2$, the posterior for $\mu$ is $N\!\left(\bar y,\ \sigma^2/n\right)$, and conditional on $\mu$, the posterior for $\sigma^2$ (using an inverse-gamma prior, omitted here for brevity) works out to something proportional to an inverse-gamma with rate depending on $\sum (y_i-\mu)^2$.

Iteration 1. Start with $\sigma^{2(0)} = 1.0$. Draw $\mu^{(1)} \sim N(0.80,\ 1.0/5) = N(0.80, 0.20)$ — suppose this draw comes out $\mu^{(1)} = 0.65$. Now compute $\sum(y_i - 0.65)^2 = (0.55)^2+(-1.15)^2+(0.15)^2+(1.45)^2+(-0.25)^2 = 0.303+1.323+0.023+2.103+0.063=3.814$. Draw $\sigma^{2(1)}$ from its conditional using this sum — suppose it comes out $\sigma^{2(1)} = 0.85$.

Iteration 2. Draw $\mu^{(2)} \sim N(0.80,\ 0.85/5) = N(0.80, 0.17)$ — suppose $\mu^{(2)} = 0.92$. Recompute $\sum(y_i-0.92)^2$ and draw $\sigma^{2(2)}$ from the updated conditional. Notice each new $\mu$ draw uses whatever $\sigma^2$ was most recently sampled, and each new $\sigma^2$ draw uses whatever $\mu$ was most recently sampled — the chain never needs the full joint distribution of $(\mu, \sigma^2)$ together, only these two conditionals, each of which is a standard, easy-to-sample distribution.

Worked example 2: why the chain still converges even from a bad start

Start the same sampler from a deliberately bad initial guess, $\mu^{(0)} = 10.0$ (wildly off — no return in the data is anywhere near 10%). Iteration 1 draws $\sigma^{2(1)}$ conditional on $\mu=10.0$: $\sum(y_i-10)^2$ is enormous (roughly $402$), so the conditional posterior for $\sigma^2$ comes out very large, reflecting "if $\mu$ really were 10, the data would look wildly overdispersed." But the next step draws $\mu^{(1)} \sim N(\bar y, \sigma^{2(1)}/n) = N(0.80,\ \text{large}/5)$ — centered correctly at $0.80$ regardless of how bad the starting $\mu$ was, just with a temporarily wide spread. Within a handful of iterations, $\mu$ has moved to sit near $0.80$ and $\sigma^2$ has shrunk back toward a value consistent with the data's actual spread — the chain "forgets" its bad start exponentially fast. This is exactly why practitioners discard the first several hundred iterations (burn-in): early samples are still contaminated by the arbitrary starting point, and only the samples after convergence are treated as genuine posterior draws.

The explorer's running average settling down after an initial noisy stretch is the same phenomenon as burn-in: early iterations reflect the starting point more than the truth, and only the later, settled portion of the chain should be trusted.

What this means in practice

• Hierarchical Bayesian models are Gibbs sampling's home turf — the population mean, population variance, and every group-level effect each have simple conditional distributions even though their joint posterior is intractable, which is exactly the setup Gibbs sampling was built for.
• Convergence diagnostics matter. Trace plots (the sampled value over iterations), the Gelman-Rubin $\hat R$ statistic comparing multiple independent chains, and effective sample size all check whether the chain has actually settled — treating early, unconverged samples as valid posterior draws silently corrupts every downstream summary.
• Samples are autocorrelated, not independent. Consecutive Gibbs draws are correlated with each other (each depends on the last), so the effective number of independent samples is smaller than the raw count — thinning (keeping every $k$-th sample) or simply running longer chains compensates.
• When conditionals aren't simple, Gibbs doesn't apply directly — that's the regime for Metropolis-Hastings or, for speed at scale, variational inference, which trades exactness for a much faster approximation.