A Worked Gibbs Sampler for a Bivariate Gaussian

From a Gaussian Target to Conditional Updates

A Gibbs sampler replaces a difficult joint draw with a sequence of draws from conditional distributions. For a bivariate Gaussian target, both conditional distributions can be found explicitly, so the algorithm is simple to implement. Yet the draws within its output sequence are not independent: the target correlation determines how much information one coordinate carries into the next sweep.

Let \((\Omega,\mathcal{F},\mathbb{P})\) be a probability space, and consider a target random vector \((X,Y)\) with a nondegenerate bivariate normal distribution. Write its mean vector as \((\mu_x,\mu_y)\), its standard deviations as \(\sigma_x>0,\sigma_y>0\), and its correlation as \(\rho\), where \(|\rho|<1\). The strict inequality ensures that the covariance matrix is positive definite. Define the standardized coordinates \(z_x=(x-\mu_x)/\sigma_x\) and \(z_y=(y-\mu_y)/\sigma_y\). The joint density is proportional to

To condition on \(Y=y\), regard \(z_y\) as fixed and complete the square in \(z_x\). The identity \(z_x^2-2\rho z_xz_y+z_y^2=(z_x-\rho z_y)^2+(1-\rho^2)z_y^2\) shows that the terms depending on \(x\) have the kernel of a normal density with mean \(\rho z_y\) and variance \(1-\rho^2\) on the standardized scale. Conditioning on \(X=x\) and completing the square in \(z_y\) gives the other conditional.

Both conditional variances are positive under the stated assumption. In a Gibbs sweep, first draw \(X_{n+1}\) from the conditional distribution given the current \(Y_n\), then draw \(Y_{n+1}\) given the newly updated \(X_{n+1}\). The order matters when writing the updates: the second draw conditions on \(X_{n+1}\), not \(X_n\).

The Systematic Gibbs Sampler

For a direct implementation, use two sequences \((Z_n)_{n\ge1}\) and \((W_n)_{n\ge1}\) of independent standard normal random variables, with all these variables independent of the initial state \((X_0,Y_0)\). The process \((X_n,Y_n)_{n\ge0}\) has state space \(S=\mathbb{R}^2\) and discrete time \(n\). For each \(n\ge0\), make one complete sweep using

These equations implement the two conditional normal draws because a normal variable with mean \(m\) and variance \(v\) can be generated as \(m+\sqrt{v}\,Z\) for a standard normal \(Z\). Fresh independent noise is used for each update. The state after the second update is \((X_{n+1},Y_{n+1})\).

Proof. Let \(f(x,y)\) be the joint target density, \(f_Y(y)\) its \(Y\)-marginal density, and \(f_{X\mid Y}(x\mid y)\) its conditional density. The factorization \(f(x,y)=f_Y(y)f_{X\mid Y}(x\mid y)\) follows from the definition of conditional density. If the initial pair has density \(f\), then after replacing \(X\) by a draw from its conditional distribution given the retained value \(Y=y\), the density of the pair is \(f_Y(y)f_{X\mid Y}(x\mid y)=f(x,y)\). Thus the \(X\)-update preserves the joint target. The same argument with \(X\) retained and \(Y\) replaced uses \(f(x,y)=f_X(x)f_{Y\mid X}(y\mid x)\), so the \(Y\)-update also preserves it. Applying these two preserving updates in sequence leaves the density equal to \(f\). \(\square\)

This result says that the target is invariant: if the sampler starts with the target distribution, each complete-sweep state has that distribution. It does not say that successive states are independent. The dependence can be calculated directly from the update equations.

Dependence Across Sweeps

For clarity, now consider the standard bivariate normal target with zero means and unit variances, retaining \(|\rho|<1\). In this case the sweep equations are \(X_{n+1}=\rho Y_n+\sqrt{1-\rho^2}Z_{n+1}\) and \(Y_{n+1}=\rho X_{n+1}+\sqrt{1-\rho^2}W_{n+1}\). Let \(\mathcal{F}_n\) be the sigma-field generated by \((X_0,Y_0)\) and all the noises \(Z_1,W_1,\ldots,Z_n,W_n\). The current state is measurable with respect to \(\mathcal{F}_n\), and \(Z_{n+1},W_{n+1}\) are independent of \(\mathcal{F}_n\) and of each other.

Substitution of the first update into the second gives

The combined noise is independent of \(\mathcal{F}_n\), has mean zero, and has variance \((1-\rho^2)(\rho^2+1)=1-\rho^4\). Consequently this coordinate itself follows a one-dimensional autoregression: its next value depends on its current value and a fresh independent normal innovation. In particular, the conditional-mean recursion can be iterated legitimately; the independence of the fresh noise from the full history supplies the needed conditioning step.

Proof. Write \(\varepsilon_{n+1}=\sqrt{1-\rho^2}(\rho Z_{n+1}+W_{n+1})\). The independence and zero mean just established give \(\mathbb{E}[\varepsilon_{n+1}\mid\mathcal{F}_n]=0\), so the recurrence implies \(\mathbb{E}[Y_{n+1}\mid\mathcal{F}_n]=\rho^2Y_n\). To iterate, suppose \(\mathbb{E}[Y_{n+k}\mid\mathcal{F}_n]=\rho^{2k}Y_n\). The recurrence at time \(n+k\), together with the independence and zero conditional mean of \(\varepsilon_{n+k+1}\) given \(\mathcal{F}_{n+k}\), gives \(\mathbb{E}[Y_{n+k+1}\mid\mathcal{F}_{n+k}]=\rho^2Y_{n+k}\). Taking conditional expectation given \(\mathcal{F}_n\) and using the tower property yields \(\mathbb{E}[Y_{n+k+1}\mid\mathcal{F}_n]=\rho^2\mathbb{E}[Y_{n+k}\mid\mathcal{F}_n]=\rho^{2(k+1)}Y_n\). The identity holds at \(k=0\), so induction proves the conditional-mean formula for every \(k\). If the chain is initialized in the target, invariance gives \(\mathbb{E}Y_n=0\) and \(\operatorname{Var}(Y_n)=1\) at every sweep. Multiplying the conditional-mean identity by \(Y_n\) and taking expectations gives \(\operatorname{Cov}(Y_{n+k},Y_n)=\rho^{2k}\mathbb{E}[Y_n^2]=\rho^{2k}\). Division by the two standard deviations, both equal to \(1\), gives the stated correlation. \(\square\)

Interpreting the Dependence and Avoiding a Common Error

The parameter \(\rho\) describes the correlation between \(X\) and \(Y\) in a single target draw. By contrast, \(\rho^2\) is the coefficient governing the dependence from one completed sweep to the next for the \(Y\)-coordinate in this update order. Its lag-\(k\) correlation is \(\rho^{2k}\) when the chain is stationary. Thus a target with \(\rho\) close to \(1\) or \(-1\) produces a coordinate that changes slowly across sweeps: in either case \(\rho^2\) is close to \(1\). When \(\rho=0\), the recurrence has no dependence on \(Y_n\), and successive \(Y\)-values are independent under stationary initialization.

Common mistake. It is incorrect to assert that the lag-one correlation of \(Y\) is \(\rho\) just because each conditional mean contains the coefficient \(\rho\). For example, if \(\rho=\tfrac12\), the two updates give \(X_{n+1}=\tfrac12Y_n+\sqrt{\tfrac34}Z_{n+1}\) and \(Y_{n+1}=\tfrac12X_{n+1}+\sqrt{\tfrac34}W_{n+1}\). Substitution makes the coefficient of \(Y_n\) equal to \(\tfrac12\cdot\tfrac12=\tfrac14\), not \(\tfrac12\). The independent noise has zero covariance with \(Y_n\), so the stationary lag-one correlation is \(\tfrac14\). The first \(\rho\) transmits \(Y_n\)'s information to \(X_{n+1}\); the second transmits that information back to \(Y_{n+1}\).

The two conditional formulas, the order of updates, and the fresh-noise independence together determine the sampler's dependence. In particular, the one-step conditional mean alone is not enough to justify iterating a conditional-mean formula: here the explicit autoregressive representation and independence from \(\mathcal{F}_n\) provide that justification. The formulas require a nondegenerate target with \(|\rho|<1\); at \(|\rho|=1\), the covariance matrix is singular and the stated positive-variance conditional updates do not apply.