Implementing MCMC

Question

I'm writing an MCMC algorithm in R and I'm wondering about the following: say we have two parameters, $\theta_1$ and $\theta_2$. I want to update each one at a time from the corresponding posterior conditional distributions. Say $\theta_1^{(0)}$ and $\theta_2^{(0)}$ are the initial values. Then at iteration 1 I update first $\theta_1$ from

$\theta_1^{(1)} \sim f(\theta_1 | \theta_2^{(0)}, \text{Data})$

Now when updating $\theta_2$, I should use

$\theta_2^{(1)} \sim f(\theta_2 | \theta_1^{(1)}, \text{Data})$

Is there a theoretical problem if when updating $\theta_2$ I use $\theta_1^{(0)}$ I instead of $\theta_1^{(1)}$? That is

$\theta_2^{(1)} \sim f(\theta_2 | \theta_1^{(0)}, \text{Data})$?

Answer 1

Yes, there is a big problem with that. You don't have the correct stationary distribution anymore when you do that. The relevant question to ask is whether the transition kernel leaves the posterior invariant i.e. $\int K(\theta \mid \theta^{old}) p(\theta^{old}) \ d\theta^{old} = p(\theta)$. The transition kernel you are proposing is $$ K(\theta_1, \theta_2 \mid \theta_1^{old}, \theta_2^{old}) = p(\theta_1 \mid \theta_2^{old}, y) p(\theta_2 | \theta_1^{old}, y). $$ This does not leave the posterior $p(\theta_1, \theta_2 \mid y)$ invariant: $$ \int \int p(\theta_1\mid\theta_2^{old}, y) p(\theta_2\mid\theta_1^{old}, y) p(\theta_1^{old}, \theta_2^{old} \mid y) \ d\theta_1^{old} d\theta_2^{old} \ne p(\theta_1, \theta_2 \mid y). $$ Actually, it leaves a different distribution invariant: $$ \int\int p(\theta_1\mid\theta_2^{old}, y) p(\theta_2\mid\theta_1^{old}, y) p(\theta_2^{old}\mid y)p(\theta_1^{old}\mid y) \ d\theta_1^{old}d\theta_2^{old} = p(\theta_1\mid y)p(\theta_2 \mid y) $$ which is not the target. Hence none of the usual theory applies - actually, it seems to me like you will effectively be running two separate, correct, Gibbs samplers so in some sense you could still use the output if you ran this chain but it would be with a different justification and it might not extend to $3$ or more parameters. Compare with the correct transition kernel, $K(\theta_1, \theta_2 \mid \theta_1^{old}, \theta_2^{old}) = p(\theta_1 \mid \theta_2)p(\theta_2 \mid \theta_1^{old})$: $$ \int \int p(\theta_1\mid\theta_2, y) p(\theta_2\mid\theta_1^{old}, y) p(\theta_1^{old}, \theta_2^{old} \mid y) \ d\theta_2^{old} d\theta_1^{old} \\ = p(\theta_1\mid\theta_2, y) \int p(\theta_1^{old}, \theta_2 \mid y) \ d\theta_1^{old} \\ = p(\theta_1\mid \theta_2, y) p(\theta_2 \mid y) = p(\theta_1, \theta_2 \mid y) $$

Answer 2

For that first iterate? That shouldn't matter, since in effect it just affects your starting point - as long as subsequently you sample correctly, everything should still work as intended - the convergence still works and so on. (It's useless but should be harmless.)

If after that first time through, your scheme is doing one of these:

(1) $\theta_1^{(i+1)} \sim f(\theta_1 | \theta_2^{(i)}, \text{Data})$

$\quad\,\,\theta_2^{(i+1)} \sim f(\theta_2 | \theta_1^{(i+1)}, \text{Data})$

OR

(2) $\theta_2^{(i+1)} \sim f(\theta_2 | \theta_1^{(i)}, \text{Data})$

$\quad\,\,\theta_1^{(i+1)} \sim f(\theta_1 | \theta_2^{(i+1)}, \text{Data})$

then you have a Gibbs sampler and from the convergence of the Markov Chain to its stationary distribution, you'll in the long run be sampling from $f(\theta_1 , \theta_2| \text{Data})$. It shouldn't matter that you did something nonstandard on the first step.