I'm currently using a bayesian approach to estimate parameters for a model consisting of several ODEs. As I have 15 parameters to estimate, my sampling space is 15-dimensional and my searched for posterior distribution seems to have many local maxima which are very isolated by large regions of very low probability.
This leads to mixing problems of my Monte Carlo chains as it is very unlikely that one chain "jumps" out of one local maximum and accidentially hits one of the other maxima.
There seems to be a lot of research in this area as its easy to find papers dealing with this problem (see below), but finding an actual implementation is hard. I only found packages relating to molecular dynamics, but not bayesian inference. Are there implementations of (MC)MC samplers which are able to deal with isolated local maxima out there?
I am forced to work with Matlab as that's what my ODE model is written in, so proposals regarding Matlab are most welcome ;-). However if there's a "killer app" in some other language, maybe I can convince my PI to switch ;-).
I'm currently working with a Delayed-Rejection/Adaptive Monte Carlo sampler written by Haario, Laine et al., and that's also the only sampler I could find so far which is more sophisticated than the standard Metropolis-Hastings algorithm
Notable approaches seem to be:
EDIT Updated on 2017-Mar-07 with what I have learnt in the meanwhile
Multiple similar chains with different starting points
Inter-Chain adaption. Use the empirical covariance matrix of the pooled samples generated by multiple independent chains to update the covariance matrices of the chain's proposal distributions.(1)
Multiple Chains with different tempering
Tempering: Some sort of "temperature" seems to modify the posterior landscape, making mixing of chains more probable. (I have not dived into this very much yet)(1) The purpose of tempering is to flatten the (high-dimensional) probability landscape formed by the posterior probability distribution. It is commonly accomplished by taking the posterior probability to the power of $1/T$, where the posterior landscape is flattened for $T>1$ (3, p.298). This means, instead of calculating the posterior probability $p(\theta\mid D)$ of a state $\theta$, given data $D$ the tempered posterior probability is calculated
$$p(\theta\mid D)^{1/T} \propto \left( p(D\mid\theta)\cdot p(\theta)\right)^{1/T} $$
The higher $T$ is chosen, the flatter and broader peaks in the probability landscape become. Therefore, higher values of $T$ lead to a higher probability of the sampler to switch from one local maximum to another. However, $p(\theta\mid D)^{1/T}$ is not the posterior distribution searched for if $T\neq1$. Therefore, the chain of samples of that distribution must be used to enable sampling from $p(\theta\mid D)$ afterwards.
Samples from the original, untempered posterior distribution, given samples from a tempered version of that distribution can be obtained by several methods:
Metropolis coupled MCMC Run multiple chains concurrently, each one having a different but constant value for $T$. Switch the states of two chains probabilistically. Only use the samples from the chain with $T=1$ for downstream estimations; the other chains just make sure that all peaks are sampled. Ref. (4) has a parallel algorithm and cites a conference article and a textbook for the idea (5,6)
Small-World MCMC. Sampler switches between two proposals. Most often a proposal distribution with small variance is used, seldomly a proposal with a big variance is used. Choice between these two proposals is stochastic. Proposals with big variance may also be drawn from another chain which only makes very big jumps, sampling as much as possible of the sample space in a rough fashion.(2,7)
Hamiltonian Monte Carlo (HMC)
I don't know much about that, but the No-U-Turn sampler (NUTS) from JAGS seems to use it. See ref. (8). Alex Rogozhnikov has created a visual tutorial on the topic.
References:
(1) Craiu et al., 2009: Learn From Thy Neighbor: Parallel-Chain and Regional Adaptive MCMC. J Am Stat Assoc 104:488, pp. 1454-1466. http://www.jstor.org/stable/40592353
(2) Guam et al., 2012: Small World MCMC with tempering: Ergocity and spectral gap. https://arxiv.org/abs/1211.4675 (only on arXiv)
(3): Brooks et al. (2011). Handbook of Markov Chain Monte Carlo. CRC press.
(4): Altekar et al. (2004): Parallel Metropolis coupled Markov chain Monte Carlo for Bayesian phylogenetic inference. Bioinformatics 20 (3) 2004, pp. 407–415, http://dx.doi.org/10.1093/bioinformatics/btg427
(5): Geyer CJ (1991) Markov chain Monte Carlo maximum likelihood. In: Keramidas (ed.), Computing Science and Statistics: Proceedings of the 23rd Symposium on the Interface. Interface Foundation, Fairfax Station, pp. 156–163.
(6): Gilks WR and Roberts GO (1996). Strategies for improving MCMC. In: Gilks WR, Richardson S and Spiegelhalter (eds) Markov chain Monte Carlo in Practice. Chapman&Hall, p. 89–114.
(7):Guan Y, et al. Markov Chain Monte Carlo in small worlds. Statistics and Computing (2006) 16 (2), pp. 193-202. http://dx.doi.org/10.1007/s11222-006-6966-6
(8): Hoffmann M and Gelman A (2014):The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo. Journal of Machine Learning Research, 15, pp. 1351-1381. https://arxiv.org/abs/1111.4246
