4
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Have you tried different starting values for $\theta$? $\endgroup$ – Jens Kouros Mar 8 '15 at 9:09
  • $\begingroup$ The denominator in the Bayes theorem is $p(y)$, not $p(\theta)$ $\endgroup$ – alberto Mar 8 '15 at 9:19
  • $\begingroup$ 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. $\endgroup$ – Hunaphu Mar 8 '15 at 17:24
  • $\begingroup$ Yes, thanks for the catch alberto and @Jens-Kouros. Edited question regarding starting values of $\theta$ $\endgroup$ – bill_e Mar 8 '15 at 17:27
  • $\begingroup$ @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? $\endgroup$ – bill_e Mar 8 '15 at 17:47
7
$\begingroup$

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

$\endgroup$
  • $\begingroup$ 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 $\endgroup$ – bill_e Mar 8 '15 at 16:35
  • $\begingroup$ How did you factor $P(\theta | x, y)$ into $P(x | \theta)P(\theta)$? $\endgroup$ – bill_e Mar 8 '15 at 17:31
  • $\begingroup$ See the Appendix :) $\endgroup$ – alberto Mar 8 '15 at 17:41
  • $\begingroup$ 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 $\endgroup$ – alberto Mar 8 '15 at 17:49
  • $\begingroup$ 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)$. $\endgroup$ – bill_e Mar 8 '15 at 17:51
4
$\begingroup$

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.

representation of the log-posterior against a 10⁴ MCMC sample produced by the above code

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – bill_e Mar 8 '15 at 19:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.