Solution (source code)

= Solution

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 is
$$
Q(x)=\frac{z_1^2-2\rho z_1z_2+z_2^2}{1-\rho^2}
=\frac{(z_1-\rho z_2)^2}{1-\rho^2}+z_2^2.
$$
Holding $x_2$ fixed and normalizing the term depending on $x_1$ gives
$$
\boxed{X_1\mid X_2=x_2\sim
N\!\left(\mu_1+\rho\frac{\sigma_1}{\sigma_2}(x_2-\mu_2),
\sigma_1^2(1-\rho^2)\right).}
$$
The analogous completion of the square gives
$$
\boxed{X_2\mid X_1=x_1\sim
N\!\left(\mu_2+\rho\frac{\sigma_2}{\sigma_1}(x_1-\mu_1),
\sigma_2^2(1-\rho^2)\right).}
$$

Using <independent> draws $\eta_{1,r},\eta_{2,r}$ from the <standard normal distribution>, a sweep is
$$
X_1^{(r)}=\mu_1+\rho\frac{\sigma_1}{\sigma_2}(X_2^{(r-1)}-\mu_2)
+\sigma_1\sqrt{1-\rho^2}\,\eta_{1,r},
$$
$$
X_2^{(r)}=\mu_2+\rho\frac{\sigma_2}{\sigma_1}(X_1^{(r)}-\mu_1)
+\sigma_2\sqrt{1-\rho^2}\,\eta_{2,r}.
$$
Use 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 yields
$$
Z_{2,r}=\rho^2Z_{2,r-1}
+\sqrt{1-\rho^2}(\rho\eta_{1,r}+\eta_{2,r}).
$$
The 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>.