A bivariate-normal Gibbs sampler

Topic: Markov chains, Gibbs sampling, and Metropolis–Hastings · more questions on this topic
‹ Adapt, discard, fix Metropolis–Hastings for the logit-normal posterior ›

Let {\boldsymbol{X}}=(X_1,X_2)^{\mathsf T} have a bivariate normal distribution with mean zero, unit marginal variances, and correlation \rho, so that

{\boldsymbol{X}}\sim\operatorname{N}_2\!\left( \begin{bmatrix}0\\0\end{bmatrix}, \begin{bmatrix}1&\rho\\ \rho&1\end{bmatrix} \right),\qquad |\rho|<1.

  1. Derive the conditional distributions of X_1\mid X_2 and X_2\mid X_1.

  2. Write pseudocode for a systematic-scan Gibbs sampler.

  3. Implement it in R or Python.

  4. From a chosen initial vector, generate 10{,}000 dependent draws for \rho=0, 0.7, and 0.97. Plot the draws and explain how \rho affects serial movement even though the stationary marginal variance remains one.

Question. Let {\boldsymbol{X}}=(X_1,X_2)^{\mathsf T} have a bivariate normal distribution with mean zero, unit marginal variances, and correlation \rho, so that

{\boldsymbol{X}}\sim\operatorname{N}_2\!\left( \begin{bmatrix}0\\0\end{bmatrix}, \begin{bmatrix}1&\rho\\ \rho&1\end{bmatrix} \right),\qquad |\rho|<1.

  1. Derive the conditional distributions of X_1\mid X_2 and X_2\mid X_1.

  2. Write pseudocode for a systematic-scan Gibbs sampler.

  3. Implement it in R or Python.

  4. From a chosen initial vector, generate 10{,}000 dependent draws for \rho=0, 0.7, and 0.97. Plot the draws and explain how \rho affects serial movement even though the stationary marginal variance remains one.

Solution. The joint density has kernel

p(x_1,x_2)\propto \exp\!\left[-\frac{x_1^2-2\rho x_1x_2+x_2^2} {2(1-\rho^2)}\right].

When x_2 is held fixed, \begin{align*} x_1^2-2\rho x_1x_2+x_2^2 &=(x_1-\rho x_2)^2-\rho^2x_2^2+x_2^2\\ &=(x_1-\rho x_2)^2+(1-\rho^2)x_2^2. \end{align*} The final term is constant as a function of x_1, so \begin{align*} p(x_1\mid x_2) &\propto \exp\!\left[-\frac{(x_1-\rho x_2)^2}{2(1-\rho^2)}\right],\\ &\boxed{X_1\mid X_2=x_2\sim\operatorname{N}(\rho x_2,1-\rho^2)}. \end{align*} Symmetry gives

\boxed{X_2\mid X_1=x_1\sim\operatorname{N}(\rho x_1,1-\rho^2)}.

Starting from (x_1^{(0)},x_2^{(0)}), one systematic Gibbs sweep is \begin{align*} x_1^{(t)}&\sim\operatorname{N}\!\left(\rho x_2^{(t-1)},1-\rho^2\right),\\ x_2^{(t)}&\sim\operatorname{N}\!\left(\rho x_1^{(t)},1-\rho^2\right). \end{align*} The newest available coordinate is used in the second update.

gibbs_bvn <- function(rho, iterations = 10000, initial = c(1, 1)) {
  stopifnot(abs(rho) < 1)
  out <- matrix(NA_real_, iterations, 2,
                dimnames = list(NULL, c("x1", "x2")))
  out[1, ] <- initial
  conditional_sd <- sqrt(1 - rho^2)

  for (t in 2:iterations) {
    out[t, 1] <- rnorm(1, rho * out[t - 1, 2], conditional_sd)
    out[t, 2] <- rnorm(1, rho * out[t, 1], conditional_sd)
  }
  out
}

set.seed(3500)
par(mfrow = c(1, 3))
for (rho in c(0, 0.7, 0.97)) {
  z <- gibbs_bvn(rho)
  plot(z, pch = 20, cex = 0.35,
       xlab = expression(x[1]), ylab = expression(x[2]),
       main = bquote(rho == .(rho)))
}
par(mfrow = c(1, 1))

When \rho=0, each conditional is \operatorname{N}(0,1) and a sweep produces a new independent pair. When \rho=0.97, the conditional variance is

1-0.97^2=1-0.9409=0.0591,

far smaller than the marginal variance one. The chain therefore moves in small, strongly correlated steps along the narrow diagonal ellipse. The stationary marginal distributions remain \operatorname{N}(0,1); what deteriorates is the effective information per iteration.