Minimal example for using the RedBlueMove

[mmenen]

Hello everyone,

so I wanted to try out the RedBlueMove to try and speed up my sampling a bit.

However, when using the command

sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob, moves=emcee.moves.RedBlueMove())

i'm getting a "Not Implemented" error message.

I guess this is because I haven't defined a propose (as listed under the RedBlueMove), but I'm also not sure how to do this.

Is there someone who already uses this successfully?

Kind regards

Marco

[mmenen]

Another question:

I have a code with more than 10 free parameters and I plotted the free parameters against the corresponding chi^2 value.

The model I'm using is sort of a black box model where I don't have control over, but the parameter space appears to be very complex.

As one can see in the attached picture of one of the plots, Emcee samples the regions of high likelyhood very well, but the points near the end of the 1- or 2-sigma region (green and yellow) are sparse.

Is there a possibility to modify the update rule of the StretchMove algorithm to include some sort of entropy or random walk behaviour, so that Emcee is more likely to also explore the less "border regions"?

Thanks in advance

Marco

[Dan Foreman-Mackey]

Hi Marco,

1. The "RedBlue" move is an abstract base class that implements the logic that is used by many other moves. You should never need to instantiate it directly. What made you want to use it? I do see that the API docs are not as clear as perhaps they should be.

2. emcee is designed to sample the target distribution that you hand it. If you want to sample a different target distribution, it's up to you to specify that, but you'll need to think about *why* you want to sample differently and how to interpret the results that you get.

Best,

Dan

--
You received this message because you are subscribed to the Google Groups "emcee users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to emcee-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/emcee-users/825208ef-63cd-4173-b834-995d57af739fn%40googlegroups.com.

--

Dan Foreman-Mackey

Associate Research Scientist

Flatiron Institute

http://dfm.io

[William Heymann]

The first thing I would do is stop using the StretchMove at least on all my systems the DEmove is VASTLY superior. I also have some complex spaces and typically DEmove is 10x more effective or more at sampling. I would use the following combination of moves based on the DEsnooker paper.

Also look at your integrated autocorrelation time. Convergence is very slow on the tails of the distribution anyways and in high dimensions it gets much harder.

Also look at the arviz library. It has tools to read emcee data and you can see what your confidence is on your posterior distribution. You might just not have enough samples yet. I have some systems where I needed to run for 59,000 steps.

moves=[

                (DESnookerMove(), 0.1),

                (DEMove(), 0.9 * 0.9),

                (DEMove(gamma0=1.0), 0.9 * 0.1),

            ]

--

[mmenen]

Hello,

first of all a (late) thank you for your answers.

Looking at the autocorrelation time definitely helped, although the sampling didn't really get any more efficient with the combined DE moves.

There is however another strange thing going on in my program. The stored samples in the chain which are accessed by sampler.get_chains are not equal to the one printed out while actually doing the sampling.

Here is a minimal example of what I'm doing:

def lnlike(theta):

    a,b = theta

    print(a,b)

    lnlike = my_program.initialize(a,b) #I call an external program to compute my log likelihood
    return lnlike

def lnprior(theta):
    a,b = theta
    if a > 1. or b > 1. or a < 0. or b < 0.:
        return -np.inf
    else:
        return 0.0

def lnprob(theta):
    lp = lnprior(theta)
    if lp < -1e+9 :
        return -np.inf
    else:
        return lp + lnlike(theta)

initial = np.array([0.5,0.5])

ndim = len(initial)

p0 = [np.array(initial) + 1e-2 * np.random.randn(ndim) for i in range(nwalkers)]

def main(p0,nwalkers,Npoints,Ntune,ndim,lnprob):

    sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob)

    p0, prob, state = sampler.run_mcmc(p0, Ntune)

    sampler.reset()

    pos, prob, state = sampler.run_mcmc(p0, Npoints, progress=True)

    return sampler, pos, prob, state

sampler, pos, prob, state = main(p0,nwalkers,Npoints,Ntune,ndim,lnprob)

samples = sampler.get_chain(flat=True)

a = samples[:,0]

b = samples[:,1]

Now I have the arrays for a and b from calling the samples as well as the values for a and b that get printed out when the lnlike function is called.

The results (for a) look like this:

a from samples / a from print statement:

0.086 / 0.086

0.114 / 0.114

0.088 / 0.093 (?)

0.132 / 0.112 (?)

0.089 / 0.089 (suddenly they are equal again)

Another thing is that the array from samples has the correct length of nwalkers*Npoints, while the lnlike function is only called some less number of times. This is due to the lnprior restrictions on the parameters I guess.

Help with this would be much appreciated, I can't figure out why the values are different.

I need them because I also want to write out additional values from my program which are wrong when the values from a and b differ.

Kind regards,

Marco

[Dan Foreman-Mackey]

Hi Marco,

This is not surprising because MCMC only "accepts" a (often small) fraction of the proposed steps. If you want to track some quantities alongside you samples, use the blobs feature: https://emcee.readthedocs.io/en/stable/user/blobs/

Dan

To view this discussion on the web visit https://groups.google.com/d/msgid/emcee-users/8540f984-9063-420d-b5a0-e441f220d7e3n%40googlegroups.com.

--

Dan Foreman-Mackey

Associate Research Scientist

Flatiron Institute

https://dfm.io

[mmenen]

Hi all,

this might be a "dumb" question but I have some problems understanding the acceptance ratio q of the Stretch Move.

When proposing a new point Y = X_j + Z * (X_k - X_j), so Z > 1 means that X_k will run away from X_j and vice versa.
Now the acceptance ratio is q = Z^(n-1) * p(Y) / p(X_k), so it heavily depends on Z. For high dimensions n, the term Z^(n-1) dominates q.
To me this looks like X_k will reject almost all moves towards X_j (as Z < 1) and accept almost all moves away from it, as p(Y) / p(X_k) barely has an impact then.
I know that q has to look the way it is because of detailed balance, but it still feels weird as the sampler clearly collects all walkers at regions of high likelihood.
Am I missing a log() term somewhere? I took all information directly from your paper https://arxiv.org/pdf/1202.3665.pdf

Kind regards

Marco