version 1.8.0.0 (372 KB) by
Aslak Grinsted

An affine invariant ensemble Markov Chain Monte Carlo sampler

Cascaded affine invariant ensemble MCMC sampler. "The MCMC hammer"

gwmcmc is an implementation of the Goodman and Weare 2010 Affine

invariant ensemble Markov Chain Monte Carlo (MCMC) sampler. MCMC sampling

enables bayesian inference. The problem with many traditional MCMC samplers

is that they can have slow convergence for badly scaled problems, and that

it is difficult to optimize the random walk for high-dimensional problems.

This is where the GW-algorithm really excels as it is affine invariant. It

can achieve much better convergence on badly scaled problems. It is much

simpler to get to work straight out of the box, and for that reason it

truly deserves to be called the MCMC hammer.

(This code uses a cascaded variant of the Goodman and Weare algorithm).

USAGE:

[models,logP]=gwmcmc(minit,logPfuns,mccount,[Parameter,Value,Parameter,Value]);

INPUTS:

minit: an MxW matrix of initial values for each of the walkers in the

ensemble. (M:number of model params. W: number of walkers). W

should be atleast 2xM. (see e.g. mvnrnd).

logPfuns: a cell of function handles returning the log probality of a

proposed set of model parameters. Typically this cell will

contain two function handles: one to the logprior and another

to the loglikelihood. E.g. {@(m)logprior(m) @(m)loglike(m)}

mccount: What is the desired total number of monte carlo proposals.

This is the total number, -NOT the number per chain.

Named Parameter-Value pairs:

'StepSize': unit-less stepsize (default=2.5).

'ThinChain': Thin all the chains by only storing every N'th step (default=10)

'ProgressBar': Show a text progress bar (default=true)

'Parallel': Run in ensemble of walkers in parallel. (default=false)

OUTPUTS:

models: A MxWxT matrix with the thinned markov chains (with T samples

per walker). T=~mccount/p.ThinChain/W.

logP: A PxWxT matrix of log probabilities for each model in the

models. here P is the number of functions in logPfuns.

Note on cascaded evaluation of log probabilities:

The logPfuns-argument can be specifed as a cell-array to allow a cascaded

evaluation of the probabilities. The computationally cheapest function should be

placed first in the cell (this will typically the prior). This allows the

routine to avoid calculating the likelihood, if the proposed model can be

rejected based on the prior alone.

logPfuns={logprior loglike} is faster but equivalent to

logPfuns={@(m)logprior(m)+loglike(m)}

TIP: if you aim to analyze the entire set of ensemble members as a single

sample from the distribution then you may collapse output models-matrix

thus: models=models(:,:); This will reshape the MxWxT matrix into a

Mx(W*T)-matrix while preserving the order.

EXAMPLE: Here we sample a multivariate normal distribution.

%define problem:

mu = [5;-3;6];

C = [.5 -.4 0;-.4 .5 0; 0 0 1];

iC=pinv(C);

logPfuns={@(m)-0.5*sum((m-mu)'*iC*(m-mu))}

%make a set of starting points for the entire ensemble of walkers

minit=randn(length(mu),length(mu)*2);

%Apply the MCMC hammer

[models,logP]=gwmcmc(minit,logPfuns,100000);

models(:,:,1:floor(size(models,3)*.2))=[]; %remove 20% as burn-in

models=models(:,:)'; %reshape matrix to collapse the ensemble member dimension

scatter(models(:,1),models(:,2))

prctile(models,[5 50 95])

References:

Goodman & Weare (2010), Ensemble Samplers With Affine Invariance, Comm. App. Math. Comp. Sci., Vol. 5, No. 1, 65–80

Foreman-Mackey, Hogg, Lang, Goodman (2013), emcee: The MCMC Hammer, arXiv:1202.3665

WebPage: https://github.com/grinsted/gwmcmc

-Aslak Grinsted 2015

Aslak Grinsted (2021). Ensemble MCMC sampler (https://github.com/grinsted/gwmcmc), GitHub. Retrieved .

Created with
R2014a

Compatible with any release

**Inspired by:**
2D Bandwidth Estimator for KDE, Markov Chain Monte Carlo sampling of posterior distribution, Cornerplot

**Inspired:**
Markov Chain Monte Carlo sampling of posterior distribution

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!Create scripts with code, output, and formatted text in a single executable document.

Igor VaimanI have a serious suspicion that this very good-looking sampler is incorrect with its cascaded implementation. I believe the same consideration is referenced in the previous comment.

I have created issue on GitHub with more detailed description of what I believe is wrong here: https://github.com/grinsted/gwmcmc/issues/8

I hope to hear back from authors of this code (if it's still mainteined) and fix the problem. Meanwhile you can try using (and by all means compare with the original!) my own derivation from the sampler, with this issue improved: https://github.com/nj-vs-vh/gwmcmc2

Fahad EjazVery helpful. Just wondering, why it does not give the information of acceptance fraction at the end of simulations?

Darcy CordellThis is a great script. However, can someone elaborate on the line that reads: "lr(1)<(numel(proposedm(:,wix))-1)*log(zz(wix))". I find that the rejection rate is quite high (even for small step size) and it seems to be related to this. I looked in Goodman and Weare (2010) and couldn't find anything about the number of elements in the model and how that is relevant to the accept/reject criterion. Apologies if this is a simple question...

Liwei ChengVery helpful!

Jacob ThompsonWhenever I try to implement this function the paramater guesses do not change from the initial values. I suspect that it may have to do with how I declared logPfuns?

Sanjay ManoharSanjay ManoharNote- the parallel implementation generates one sample in each worker, then needs to coalesce before running the next sample. So only useful with very slow objective functions... Needs to be improved using spmd? if anyone knows how to do this that would be great!

name1 name2hi what is the difference between your algorithm and the rjmcmc of green?

Gilberto Gonzalez-ParraCan I get 95% confidence intervals for each parameter ?

Aslak GrinstedThe algorithm can only explore the dimensions of space spanned by the starting points of the walkers. If all the walkers start in the same point then they will never move.

Romain Fiévetworks great otherwise, nice good!

Romain FiévetI am having a problem describing a simple Gaussian prior with this code. The chain only has constant value! I tried the line_fit example and it works, but the moment I remove the randomness from the initial ensemble of walkers, it also contains only constant chains. Is this normal behavior ? Shouldnt the random walk move away from the initial position ?

Kusa AmuIs it possible to apply hierarchical bayes? If yes, could you provide an example to determine hyper-parameters of the slope from the linefit example?

That would be very helpful.

Liguang JiangNo boundary handing when new candidate is proposed!

AMRANE youcefcan help me

Nathan BowmanAslak Grinsted@Yalei: Good question. I guess this website shows uncompressed size. But it is out of my control.

Yalei DingFile Size: 1.35 MB

File ID: #49820

Version: 1.8

why the zip i downloaded was 391kb？

Aslak Grinsted@Arpad Rozsas: thank you. I have fixed the typo in the example

Arpad RozsasThanks for this great code!

I think there is a minor error in the example.

logP={@(m)-0.5*sum((m-mu)'*iC*(m-mu))}

should be:

logPfuns={@(m)-0.5*sum((m-mu)'*iC*(m-mu))}

Luis O'Shea