= Solution
At $\sigma=1$, complete the square:
$$
x^2+y^2-2\rho xy=(x-\rho y)^2+(1-\rho^2)y^2.
$$
For fixed $y$, all factors not involving $x$ enter its conditional <normalizing constant>. Interchanging the roles of the coordinates gives
$$
\boxed{X\mid Y=y\sim N(\rho y,1-\rho^2),\qquad Y\mid X=x\sim N(\rho x,1-\rho^2).}
$$
These are the <full conditional distributions> of the standard <bivariate normal distribution> with correlation $\rho$.
A systematic <Gibbs sampler> starts at a chosen $(x_0,y_0)$ and, at each sweep, draws <independent> standard normals $Z_{1,n},Z_{2,n}$ and sets
$$
x_{n+1}=\rho y_n+\sqrt{1-\rho^2}\,Z_{1,n},\qquad y_{n+1}=\rho x_{n+1}+\sqrt{1-\rho^2}\,Z_{2,n}.
$$
The second update uses the newly drawn $x_{n+1}$. Both coordinate updates preserve the target by part (a), hence so does the sweep. The <Gaussian two-coordinate Gibbs recursion> is explicitly
$$
y_{n+1}=\rho^2y_n+\sqrt{1-\rho^2}(\rho Z_{1,n}+Z_{2,n}),
$$
an <AR(1)> recursion with innovation <variance> $1-\rho^4$. Its stationary <variance> is one, and its correlation at sweep lag $k$ is $\rho^{2k}$, displaying the dependence of the sample.
The printed general-$\sigma$ expression has a separate normalization defect: for $\sigma>0$ its integral is $\sigma$, not one, since the exponent has <covariance> $\sigma^2\begin{pmatrix}1&\rho\\\rho&1\end{pmatrix}$, whose normalizing denominator is $2\pi\sigma^2\sqrt{1-\rho^2}$. Retaining that exponent requires $\sigma^2$ in the prefactor. The requested $\sigma=1$ target is already normalized, so this defect does not alter the <conditional distributions> above.
Back to article page
