(interacting) MCMC for multimodal posterior

I am trying to sample from a posterior having many modes particularly far from each others using MCMC. It appears that in most cases, only one of these modes contains the 95% hpd I am looking for. I tried to implement solutions based on tempered simulating but this does not provide satisfactory results as in practice going from one "capture range" to another is to costly.

As a consequence, it appears to me that a more efficient solution would be to run many simple MCMCs from different starting points and to dive into the dominant solution by making the MCMCs interact each other. Do you know if there is some proper way to implement such an idea ?

Note: I found that paper http://lccc.eecs.berkeley.edu/Papers/dmcmc_short.pdf (Distributed Markov chain Monte Carlo, Lawrence Murray) that looks close to what I am looking for but I really do not understand the design of the function $R_i$.

[EDIT]: the lack of answers seems to indicate that there is no obvious solution to my initial problem (making several MCMCs sampling from the same target distribution from different starting point interact each other). Is that true ? why is it so complicated ? Thanks

• Sounds like what I've previously used "Population MCMC" for. See "Estimating Bayes factors via thermodynamic integration and population MCMC" by Calderhead and Girolami. Sorry if this comes a bit late, but maybe other people will appreciate the reference! – Sam Mason Feb 23 '15 at 13:01
• Possibly the same as Population MCMC is Sequential Monte Carlo. – Robert P. Goldman Jun 6 at 3:27

Thirst of all I would recommend to look for a better method, or at least a method with more in-depth description, since the "Distributed Markov chain Monte Carlo" from the paper that you are refering doesn't seem to be clearly stated. The advantages and disadvantages are not well explored. There is a method, that showed up in arxiv quite recently called "Wormhole Hamiltonian Monte Carlo", I would recommend to check it.

Going back to the paper that you gave reference, the remote proposal $R_{i}(\theta_{i})$ is very vaguely discribed. In the application part it is described as "maximum likelihood Gaussian over the preceeding t/2 samples". Maybe this means that you average the last t/2 values of the $i^{th}$ chain? A bit hard to guess with the poor description given in the reference.

[UPDATE:] Interaction between several chains and the application of this idea to sample from posterior distribution can be found in parallel MCMC methods, for example here. However, running several chains and forcing them to interact may not fit for the multimodal posterior: for example, if there is a very pronounced region where most of the posterior distribution is concentrated the interaction of the chains may even worsten things by sticking to that specific region and not exploring other, less pronounced, regions/modes. So, I would strongly recommend to look for MCMC specifically designed for multimodal problems. And if you want to create another/new method, then after you know what is available in the "market", you can create more efficient method.

• Thanks for your help. This is not totally what I expected but maybe I am looking for unreachable things. I am currently investigating (I am not sure to be able to derive everything needed by an Hamiltonian MCMC), the so-called darting MC that I found through your references. I give you the bounty, thks again for your advices. – peuhp Jun 18 '13 at 7:37
• @peuhp, This is closely related to the Hamiltonian methods referred to above, but the Equi-Energy Sampler is an MCMC method that is specifically billed as a method for multi-modal posterior distributions. I know it operates by starting several parallel chains and allowing jumps between chains but I do not have enough knowledge of/experience with the method to post an answer. Have a look at the link. I Hope this helps. – Macro Jun 18 '13 at 14:20

You should try multinest: https://arxiv.org/pdf/0809.3437.pdf https://github.com/JohannesBuchner/MultiNest It's a bayesian inference engine that will give you parameter samples for a multimodal distribution.

The github link contains multinest source code that you compile and install as per instructions. it also has a python wrapper that's easier to use. The example codes have a prior section which serves to constrain your parameters, and a likelihood section which contains your likelihood. the settings file contains all your settings, and chains folder multinest output after fitting. it will give you samples of your parameters

• Welcome to CV and thank you for your great answer. Maybe you can improve your answer by summarizing the Githublink. – Ferdi Jul 4 '18 at 10:35
• The github link contains multinest source code that you compile and install as per instructions. it also has a python wrapper that's easier to use. The example codes have a prior section which serves to constrain your parameters, and a likelihood section which contains your likelihood. the settings file contains all your settings, and chains folder multinest output after fitting. it will give you samples of your parameters. – Eric Kamau Jul 9 '18 at 16:23
• wow. That is awesome. I upvoted your answer. In the future you can include such a description everytime you share a link. – Ferdi Jul 10 '18 at 8:07

This appears to be a difficult and ongoing problem in computational stats. However, there are a couple of less state-of-the-art methods that should work ok.

Say you have already found several distinct modes of the posterior and you are happy that these are the most important modes, and if the posterior around these modes is reasonably normal. Then you can calculate the hessian at these modes (say, using optim in R with hessian=T) and you can approximate the posterior as a mixture of normals (or t distributions). See p318-319 in Gelman et al. (2003) "Bayesian Data Analysis" for details. Then you can use the normal/t-mixture approximation as the proposal distribution in an independence sampler to obtain samples from the full posterior.

Another idea, that I haven't tried, is Annealed Importance Sampling (Radford Neal, 1998, link here).

• Another thing, if you already have a set of chains that each appear to have "converged" locally, would be to do some sort of importance-resampling from the combination of all of your chains. I'm not sure about that though. – Mat Jun 26 '15 at 14:50