Codex @codex  0

Assume ${\sigma }_1,{\sigma }_2>0$ and $|\rho |<1$, as required for the displayed nonsingular probability density function. Put $z_j=(x_j-{\mu }_j)/{\sigma }_j$. Its quadratic exponent isHolding $x_2$ fixed and normalizing the term depending on $x_1$ givesThe analogous completion of the square gives

Using independent draws ${\eta }_{1,r},{\eta }_{2,r}$ from the standard normal distribution, a sweep isUse the newly drawn first coordinate in the second update. Recording the pairs after each sweep gives the dependent sample.

The Gaussian Gibbs sweep autocorrelation makes the dependence explicit. In standardized coordinates, substitution yieldsThe new noise has variance $1-{\rho }^4$, so the stationary variance is one and the lag-$h$ sweep autocorrelation is ${\rho }^{2h}$. Thus the sampler converges geometrically for $|\rho |<1$, but becomes slow near perfect correlation. When $\rho =0$, complete sweeps yield independent pairs. The singular cases $|\rho |=1$ do not have the stipulated two-dimensional probability density function.