# MCMC chain getting stuck

I am trying to use a Metropolis-within-Gibbs type algorithm to sample $\theta$ and $x$ from the following model. Starting with Bayes theorem I can write: $$P(\theta, x | y) = \frac{P(y | x, \theta) P(\theta)}{P(y)}$$ I can easily draw samples from $P(\theta)$ and $P(x | \theta)$, so I factor $P(y | x, \theta)$ into $P(y | x) P(x | \theta)$.

Bayes theorem again is now, $$P(\theta, x | y) = \frac{P(y | x) P(x | \theta) P(\theta)} {P(y)}$$

• $P(y | x)$ is multivariate Gaussian.
• $P(x | \theta)$ is multivariate t.
• $P(\theta)$ is uniform.

My sampling steps:

1. Start with some value of $\theta$, say, $\theta_1$.

2. Gibbs step: Draw a sample $x_1$ from $P(x | \theta_1)$.

3. Plug my $x_1$ into $P(y | x)$ and evaluate (y was observed data).

4. Metropolis step: Accept $\theta$, $x_1$, if $P(y | x_1) > P(y | x_0)$. If not, accept $\theta$, $x_1$ with probability $a = \frac{P(y | x_1)}{P(y | x_0)}$.

5. Go back to step 2, using $\theta_1$ if it was accepted, otherwise, draw a new $\theta$ from the proposal, then go back to step 2.

The problem is that $\theta$ gets stuck. Sometimes I will get a draw of $x$ from $P(x | \theta)$, where that $x$ will be nearly identical to y, making result from $P(y | x)$ exceptionally high. After that happens, essentially all my steps are rejected. $a = \frac{P(y | x_1)}{P(y | x_0)}$ for pretty much all subsequent $\theta$ steps is a very small number.

Here is why I think this happens: The likelihood $P(y | x)$ is very sharply peaked. On the other hand the distribution $P(x | \theta)$ is very broad in comparison. When $\theta$ changes a small amount, $P(x | \theta)$ is essentially unchanged.

Would a different sampling strategy work better in such a situation? Let me know if further info is required.

• Have you tried different starting values for $\theta$? – Jens Kouros Mar 8 '15 at 9:09
• The denominator in the Bayes theorem is $p(y)$, not $p(\theta)$ – alberto Mar 8 '15 at 9:19
• This can happen if you have large variation in your x-sampling; you get an unlikely value so that you never accept another proposal. If you do some kind of random walk proposals then you could try to adjust the variance. – Hunaphu Mar 8 '15 at 17:24
• Yes, thanks for the catch alberto and @Jens-Kouros. Edited question regarding starting values of $\theta$ – bill_e Mar 8 '15 at 17:27
• @Hunaphu, could you elaborate more? Yes, I do have large variation in my x sampling. What do you mean by 'do some random walk proposals' to try to adjust the variance? It is what it is right? – bill_e Mar 8 '15 at 17:47

In step 4, you don't have to reject the proposal $x,\theta$ every time its new likelihood is lower; if you do so, you are doing a sort of optimization instead of sampling from the posterior distribution.

Instead, if the proposal is worse then you still accept it with an acceptance probability $a$.

With pure Gibbs sampling, the general strategy to sample this would be:

Gibbs

Iteratively sample: \begin{align} p(x | \theta, y) &\propto p(y | x) p(x |\theta)\\ p(\theta | x, y) &\propto p(x | \theta) p(\theta) \end{align}

Gibbs with Metropolis steps for non-conjugate cases:

If you some of the conditionals above is not a familar distribution (because you are multiplying non-conjugates; this is your case) you can sample with Metropolis Hastings:

• From the current $x$, generate some proposal, e.g.: $$x^* \sim \mathcal{N(x, \sigma)}$$

• Accept $x^*$ with probability [1]: $$a = min \left(1, \frac{p(x^*)}{p(x)}\right) = min \left(1, \frac{p(x^* | \theta, y)}{p(x | \theta, y)}\right)$$

[1] If the proposal distribution wasn't symmetric then there is another multiplying factor.

Appendix: $$p(\theta | x, y) = \frac{p(y|x)p(x| \theta)p(\theta)} {\int p(y|x)p(x| \theta)p(\theta) \text{d}\theta}= \frac{p(x| \theta)p(\theta)} {\int p(x| \theta)p(\theta) \text{d}\theta} \propto p(x| \theta)p(\theta)$$

• Sorry, I'll edit the question. I do what you suggest, in step 4, I refer to this as the "metropolis step ... etc". I will edit the question to be more clear – bill_e Mar 8 '15 at 16:35
• How did you factor $P(\theta | x, y)$ into $P(x | \theta)P(\theta)$? – bill_e Mar 8 '15 at 17:31
• See the Appendix :) – alberto Mar 8 '15 at 17:41
• Btw after your last edit I would say you are missing a factor in $a$. You consider that your proposal is symmetric, but it is not. See en.wikipedia.org/wiki/Metropolis%E2%80%93Hastings_algorithm – alberto Mar 8 '15 at 17:49
• My proposal distribution for proposing values of $\theta$ is symmetric. Does drawing from $P(x | \theta)$ change that? One other small point. In this, y is my "data", so I don't think I can sample $P(y | x)$ in $P(y | x)P(x | \theta)$. – bill_e Mar 8 '15 at 17:51

Here is an R code in the univariate case for the above Metropolis-within-Gibbs approach drafted by @alberto. No indication of the chain getting stuck: the acceptance rate for the $x$ component is close to 50%.

First, I picked some pseudo-values to run the algorithm:

#observation from N(x,1)
y=3.081927
#latent x from t(nu,theta,1)
nu=3


Second, I simulated the location $\theta$ from the full condition distribution, namely a Student's $t$-distribution with location parameter $x$ and the latent parameter $x$ by a Metropolis-within-Gibbs step, making a proposal from a Student's $t$-distribution with location parameter $\theta$ and accepting this proposal based on the second part of the full conditional, namely the normal pdf centred in $y$.

#Metropolis-within-Gibbs
T=10^4
mcmc=matrix(NA,T,2)
#initialisation
mcmc[1,1]=rnorm(1,mean=y)
mcmc[1,2]=rt(1,df=nu)+mcmc[1,1]
#Gibbs iterations
for (t in 2:T){
mcmc[t,1]=rt(1,df=nu)+mcmc[t-1,2] #theta
mcmc[t,2]=proposal=rt(1,df=nu)+mcmc[t,1] #x
#acceptance probability:
accept=dnorm(proposal,mean=y)/dnorm(mcmc[t-1,2],mean=y)
if (runif(1)>accept) mcmc[t,2]=mcmc[t-1,2]
}


As seen from the contour plot below, the resulting chain $(\theta_t,x_t)$ is correctly located on the highest contours of the target density.

• Cool plot! But alas I feel more confused than when I started. Column 1 of your mcmc matrix is ... $\theta$, and column 2 is $x$? The the sample from row 1 is added to row 2 and visa versa..? Lost! Thank you for your answer, I think it is my fault I am struggling to understand. – bill_e Mar 8 '15 at 19:27