covariance matrix for estimating $(\mu_1, \mu_2)^T$ with the Monte Carlo average. The proof is provided in Appendix \[subsec:Gibbs\_var\].

\[thm:bivariate\_gibbs\] For the two variable deterministic scan Gibbs sampler described by the full conditionals in and , the asymptotic covariance matrix in the CLT for the sample mean is $$\Sigma = \begin{bmatrix}
\omega_1\left(\dfrac{\omega_1\omega_2 + \rho^{2}}{\omega_1 \omega_2 - \rho^{2}}\right) & \dfrac{2 \omega_1 \omega_2\rho}{\omega_1 \omega_2 - \rho^{2}}\\
\dfrac{2\omega_1 \omega_2\rho}{\omega_1\omega_2 - \rho^{2}} & \omega_2\left(\dfrac{\omega_1\omega_2 + \rho^{2}}{\omega_1\omega_2 - \rho^{2}}\right)
\end{bmatrix} \; .$$

Theorem \[thm:bivariate\_gibbs\] is a rare result in MCMC output analysis. To the best of our knowledge, this is the first time a closed-form expression of the asymptotic covariance matrix has been obtained for a non-trivial MCMC algorithm and it allows us to compare our estimates of $\Sigma$ against the truth.

We compare the performance of RBM, ABM, and the naive estimator for two settings, $\rho = 0.5, 0.999$ with $\omega_1 = \omega_2 = 1$. For Setting 1 ($\rho = .5$), the Markov chains travel freely across the state space so that in only a few steps, independent runs of the Markov chains from over-dispersed starting values look similar to each other. For Setting 2 $(\rho = .999)$, each Markov chain moves slowly across the 