= Solution
Start from any $(a^{(0)},b^{(0)})$. For independent draws $G_{1,k},G_{2,k}$ from the <standard normal distribution> at each sweep, the <Gibbs sampler> can be implemented as
$$
\boxed{a^{(k+1)}=\frac{r-Cb^{(k)}}A+\frac{G_{1,k}}{\sqrt A},\qquad
b^{(k+1)}=\frac{s-Ca^{(k+1)}}D+\frac{G_{2,k}}{\sqrt D}.}
$$
The second update must use the newly drawn first coordinate. The transform in part (a) supplies the needed independent Gaussian draws. Each conditional update preserves the joint <posterior>, so their composition does too.
Convergence is particularly transparent here. Centering at the <posterior> mean gives
$$
b^{(k+1)}-\mu_b=\frac{C^2}{AD}(b^{(k)}-\mu_b)-\frac{C}{D\sqrt A}G_{1,k}+\frac{G_{2,k}}{\sqrt D}.
$$
Positive definiteness gives $C^2/(AD)<1$, so this is a stable Gaussian autoregression. The full sampler has the desired <posterior> as its limiting invariant law. This is the <linear contraction of a two-coordinate Gaussian Gibbs sweep>.
For a <posterior>-integrable function $\phi$, its <posterior> expectation is estimated after a burn-in $K$ by
$$
\boxed{\widehat{\mathbb E}_\pi\phi=\frac1L\sum_{k=K+1}^{K+L}\phi(a^{(k)},b^{(k)}).}
$$
The ergodic theorem justifies this Monte Carlo average. The retained values of $\phi$ also approximate its <posterior> distribution and quantiles. If a Bayesian point estimate is requested under squared-error loss, the <posterior> mean is the appropriate estimate, when its needed moments exist. An arbitrary function need not have a finite <posterior> mean; integrability must be assumed for the displayed target. Successive Gibbs draws are correlated, so uncertainty in the Monte Carlo average should use chain-aware error estimates rather than treating them as iid.
Back to article page
