I am very new to statistics and trying to implement a Gibbs sampler. However, according to wikipedia https://en.wikipedia.org/wiki/Gibbs_sampling and this discussion thread http://metaoptimize.com/qa/questions/7831/disadvatage-of-gibbs-sampler, it is possible for Gibbs sampling to get stuck in a local mode. When this is the case, how do people tend to resolve this issue?

I am attempting to sample a distribution in the context of bioinformatics, where the joint distribution is difficult to compute analytically. The reference paper is here http://bioinformatics.oxfordjournals.org/content/31/12/i97.long, though it might be too wordy to go through.

Any general directions will be appreciated.enter image description here

  • $\begingroup$ I got a "failed to load" error for that PDF. $\endgroup$ – eric_kernfeld Jan 16 '16 at 21:18
  • $\begingroup$ @eric_kernfeld Sorry for not being able to attach the full article, I cannot find attach button anywhere. But that's the Gibbs sampling part, the exponentials and the probability on RHS of the very last equations could be easily computed, though some exponentials quickly outweighs everything else making the sampler stuck in local mode. I'm only after general directions of how people tend to do in this senario. $\endgroup$ – TuanDT Jan 16 '16 at 21:30

This is certainly possible. It often happens when variables are strongly correlated.

For simplicity, consider a two-parameter model. Because Gibbs sampling alters only one variable at a time, it can only move vertically or horizontally on the Cartesian plane. It will be unable to reach regions of high posterior probability that lie diagonally from one another.

You may be able to resolve this problem by reparameterizing your model so the posterior is roughly spherical or by jointly updated blocks of correlated parameters. For an example of how to reparameterize, you could take a regression model and standardize the covariates. This makes is so that the slope won't have much effect on the fit of the intercept and vice versa. Without standardizing, data may sit far from the origin, and a small change in slope could demand a huge change in intercept for the regression line to fall across the data.

MCMC implementation is difficult and error prone, so you may want to consider tools such as STAN. STAN provides already-implemented state-of-the-art MCMC methods via a language similar in spirit to R's formula syntax. It has interfaces with many common scientific computing languages. (I am not affiliated with STAN -- merely a believer.)

  • $\begingroup$ Also, would the initial condition mater? $\endgroup$ – TuanDT Jan 16 '16 at 21:51
  • $\begingroup$ Would another sampler, I have metropolis hasting in mind, avoid this issue in general? $\endgroup$ – TuanDT Jan 17 '16 at 13:13
  • $\begingroup$ I prefer to initialize somewhere safe, such as a method-of-moments or maximum-likelihood estimate. Other people arrive at a suitable starting position via "burn-in" (running the sampler itself). The initial condition only matters a lot if the sampler sucks or if it's deliberately outrageous (theta[1] = MAX_FLOAT). Metropolis-hastings is more flexible than Gibbs; it only has trouble if you use a bad proposal. There are adaptive versions that help find a good proposal automatically. STAN uses Hamiltonian Monte Carlo, a form of metropolis-hastings with a fancy physics-inspired proposal. $\endgroup$ – eric_kernfeld Jan 17 '16 at 22:10
  • $\begingroup$ STAN seems like a good software package. Thank you for the recommendation. $\endgroup$ – TuanDT Jan 19 '16 at 21:06

Here's a simple case where Gibbs sampling gets stuck:

Imagine $(\theta_1,\theta_2)$ is a 50-50 mixture of two bivariate normals which each have independent components with variance 1.

The first mixture component is centered at $(10,-10)$. The second is centered at $(-10,10)$.

![enter image description here

If you're in the top-left (green) mode, you'll stay there. If you're in the red mode, you'll stay there instead.

Each component is a local mode in the bivariate distribution. When you start, you'll right away get stuck in one mode or the other -- you can easily move around within either mode, but not between them.

In a situation like this one, there are a variety of strategies that might be used. For example, sometimes it is feasible to sample affected variables as a block. Sometimes you can integrate some variables out, which can reduce or even fix the problem. Sometimes you might be able to reparameterize (rewriting as $(\psi_1,\psi_2)=(\theta_1+\theta_2,\theta_1-\theta_2)$ would work for the abobe example). You may be able to add a variable (or more than one) that allow movement between modes that could not previously communicate. There are numerous other strategies; some more suited to particular kinds of multiple modes.

You might also move to a completely different approach, such as Hamiltonian Monte Carlo.

  • $\begingroup$ I don't see how Hamiltonian Monte Carlo would be helpful in such case. It would also get stuck in the local mode, no? It would even get stuck if the modes were separated vertically or horizontally, while Gibbs sampler only suffers if the modes are separated diagonally (but not vertically or horizontally). $\endgroup$ – Davor Josipovic Nov 4 '18 at 20:17

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.