# Implementing MCMC

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})$?

• Not exactly related, but I wouldn't use MCMC for a two parameter problem. Quadrature would be much better. Aug 18, 2013 at 0:30
• @probabilityislogic when first learning MCMC it's often worth starting on problems you can check the answers to. Aug 18, 2013 at 2:48
• @Glen_b - I take your point, but there's nothing explicitly in the question that says this is a pure learning exercise. If it is a "real" problem, then mcmc would be crazy for a two parameter problem. MCMC should always be the method of last resort imo. Aug 18, 2013 at 3:46
• @probabilityislogic MCMC works fine for arbitrary number of parameters. Otherwise hierarchical models wouldn't work very well, since you generally need to track the posteriors for multiple parameters at each level of the model. Hell, bayesian linear regression wouldn't work either sinec you'd need to track the posterior of each coefficient. Aug 18, 2013 at 5:02
• @DavidMarx I think probabilityislogic is suggesting that it is better to use MCMC on large problems rather than very small problems that can be done to high accuracy in other ways. Aug 18, 2013 at 5:31

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)$$

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.