The random-walk Metropolis sampler in the companion post has one fiddly part: a proposal width you have to tune, and a stream of rejected moves you have to accept as waste. When a model has a helpful structure, you can avoid both. If each parameter, holding the others fixed, has a posterior that is itself a standard distribution, you can draw from those one at a time and never reject anything. That is Gibbs sampling. This post codes it in base R for the simplest two-parameter problem an ecologist meets constantly: estimating the mean and the variance of a measured trait. Along the way it shows a mistake that is easy to make, treating the variance as if it were known once you have estimated it, which quietly shrinks the interval on the mean.
A dozen individuals of a hard-to-catch passerine are measured for wing length in millimetres. We want the population mean and also the population variance, because among-individual spread is itself an ecological quantity (it feeds selection gradients, size-based models, and detection).
set.seed(2718)
n <- 12L
y <- rnorm(n, mean = 76, sd = 3.2) # a small synthetic sample, wing length (mm)
ybar <- mean(y)
s2 <- var(y) # sample variance (denominator n - 1)
s <- sqrt(s2)
round(c(n = n, ybar = ybar, s = s, s2 = s2), 3) n ybar s s2
12.000 76.786 3.253 10.579 The sample mean is 76.79 mm and the sample standard deviation is 3.25 mm. With only 12 birds, neither quantity is pinned down, and the uncertainty in the variance should leak into the uncertainty in the mean.
Take the data as normal with unknown mean and variance, and use the standard reference prior \(p(\mu, \sigma^2) \propto 1/\sigma^2\), which is flat on \(\mu\) and on \(\log \sigma\). The joint posterior has no simple one-line form, but each parameter conditional on the other does. Completing the square in \(\mu\) and collecting powers of \(\sigma^2\) gives two familiar distributions:
\[\mu \mid \sigma^2, y \;\sim\; \text{Normal}\!\left(\bar y,\; \sigma^2/n\right), \qquad \sigma^2 \mid \mu, y \;\sim\; \text{Inv-Gamma}\!\left(\tfrac{n}{2},\; \tfrac{1}{2}\textstyle\sum_i (y_i - \mu)^2\right).\]
Both are things we can sample directly. A normal draw is rnorm; an inverse-gamma draw is the reciprocal of a gamma draw, since if \(X \sim \text{Gamma}(a, \text{rate}=b)\) then \(1/X \sim \text{Inv-Gamma}(a, b)\). Gibbs simply alternates: draw a new mean given the current variance, then a new variance given the new mean, and repeat.
gibbs <- function(n_iter, y, init_sigma2 = 1) {
n <- length(y); ybar <- mean(y)
mu <- numeric(n_iter); sig2 <- numeric(n_iter)
cur_sig2 <- init_sigma2
for (t in seq_len(n_iter)) {
cur_mu <- rnorm(1, ybar, sqrt(cur_sig2 / n)) # mu | sigma^2, y
ss <- sum((y - cur_mu)^2)
cur_sig2 <- 1 / rgamma(1, shape = n / 2, rate = ss / 2) # sigma^2 | mu, y
mu[t] <- cur_mu; sig2[t] <- cur_sig2
}
data.frame(mu = mu, sigma2 = sig2)
}There is no proposal width, no accept step, and no rejected draw. The price is that every conditional had to be a distribution we can sample from, which here it is.
n_iter <- 40000L; burn <- 2000L
set.seed(90210)
chain <- gibbs(n_iter, y, init_sigma2 = 1)
post <- chain[(burn + 1):n_iter, ]
round(c(acf1_mu = acf(post$mu, 1, plot = FALSE)$acf[2],
acf1_sigma2 = acf(post$sigma2, 1, plot = FALSE)$acf[2]), 3) acf1_mu acf1_sigma2
-0.004 0.096 The lag-one autocorrelation is -0.004 for the mean and 0.096 for the variance: successive draws are practically independent, so the chain mixes far better than a tuned random walk would on the same target. The figure below shows how Gibbs gets there. It updates one coordinate at a time, so its path through the joint space is a staircase of horizontal and vertical moves rather than a diagonal drift.
sig <- sqrt(post$sigma2)
K <- 30
mu <- post$mu[1:K]; sg <- sig[1:K]
vx <- mu[1]; vy <- sg[1]
for (t in 2:K) { vx <- c(vx, mu[t], mu[t]); vy <- c(vy, sg[t-1], sg[t]) }
stair <- data.frame(x = vx, y = vy)
cloud <- data.frame(mu = post$mu, sigma = sig)[seq(1, nrow(post), by = 4), ]
ggplot() +
geom_point(data = cloud, aes(mu, sigma), colour = sage, alpha = 0.20, size = 0.5) +
geom_path(data = stair, aes(x, y), colour = ink, linewidth = 0.35) +
geom_point(data = stair[1, ], aes(x, y), colour = brick, size = 2) +
labs(x = expression(mu ~ "(mean wing length, mm)"), y = expression(sigma ~ "(mm)"),
title = "Gibbs moves one coordinate at a time",
subtitle = "grey: posterior draws; line: first 30 alternating conditional updates")
This model has a closed-form marginal for the mean. Under the reference prior, integrating the variance out of the joint posterior leaves a Student-t: \(\mu \mid y \sim t_{n-1}(\bar y,\; s^2/n)\), the same t that underlies the one-sample confidence interval. That gives an exact target for the sampler and a clean contrast with the shortcut.
mu_ci_gibbs <- quantile(post$mu, c(0.025, 0.975))
mu_ci_exact <- ybar + c(-1, 1) * qt(0.975, df = n - 1) * s / sqrt(n) # exact t interval
mu_ci_naive <- ybar + c(-1, 1) * 1.96 * s / sqrt(n) # fix sigma at s, use normal
s2_ci_gibbs <- quantile(post$sigma2, c(0.025, 0.975))
s2_ci_exact <- (n - 1) * s2 / qchisq(c(0.975, 0.025), df = n - 1) # exact for the variance
rbind(mu_gibbs = mu_ci_gibbs, mu_exact_t = mu_ci_exact, mu_naive = mu_ci_naive) 2.5% 97.5%
mu_gibbs 74.70540 78.84898
mu_exact_t 74.71949 78.85258
mu_naive 74.94576 78.62631The Gibbs interval for the mean, 74.71 to 78.85, lands on the exact t interval, 74.72 to 78.85. The variance interval matches too: Gibbs gives 5.34 to 30.34 against the exact 5.31 to 30.5. The naive interval, which fixes the standard deviation at its sample value and then treats the mean as normal, runs 74.95 to 78.63. It is narrower, by 12.3 per cent here, and the gap grows as the sample shrinks. The next figure shows why: the mean’s true posterior has the heavier tails of a t, and pretending the variance is known throws those tails away.
xs <- seq(ybar - 4.5, ybar + 4.5, length.out = 400)
sc <- s / sqrt(n)
curv <- rbind(
data.frame(x = xs, dens = dt((xs - ybar) / sc, df = n - 1) / sc, src = "exact posterior t (df 11)"),
data.frame(x = xs, dens = dnorm(xs, ybar, sc), src = "naive fixed-sigma normal"))
curv$src <- factor(curv$src, levels = c("exact posterior t (df 11)", "naive fixed-sigma normal"))
ggplot() +
geom_histogram(data = data.frame(mu = post$mu), aes(mu, after_stat(density)),
binwidth = 0.15, fill = sage, colour = paper, alpha = 0.6) +
geom_line(data = curv, aes(x, dens, colour = src, linetype = src), linewidth = 0.85) +
scale_colour_manual(values = c(forest, brick), name = NULL) +
scale_linetype_manual(values = c("solid", "longdash"), name = NULL) +
labs(x = expression(mu ~ "(mean wing length, mm)"), y = "density",
title = "The mean's posterior is a t, not a normal",
subtitle = "accounting for unknown variance widens the interval (heavier tails)")
The joint draws also answer questions no single interval can. If a downstream model needs both a large mean and a wide spread, the posterior gives the joint probability straight from the stored columns: here about 19 per cent of draws have a mean above 77 mm and a standard deviation above 3.5 mm.
Gibbs and Metropolis are not rivals so much as tools for different situations. Gibbs needs conjugate conditionals but then runs with no tuning and little autocorrelation, which is why it dominates this normal model. When a conditional is not a standard distribution, you can still take a Metropolis step for that one parameter inside an otherwise-Gibbs loop, which is the pattern the hierarchical model in this series uses for its non-conjugate pieces. Either way, an independent check on convergence is still needed: a chain that looks settled in one run is not proof, and the diagnostics post shows how to test it with several chains and a numerical statistic.