Using different moves

One of the most important new features included in the version 3 release of emcee is the interface for using different “moves” (see Moves for the API docs). To demonstrate this interface, we’ll set up a slightly contrived example where we’re sampling from a mixture of two Gaussians in 1D:

Show code cell content Hide code cell content

%config InlineBackend.figure_format = "retina"

from matplotlib import rcParams
rcParams["savefig.dpi"] = 100
rcParams["figure.dpi"] = 100
rcParams["font.size"] = 20

import numpy as np
import matplotlib.pyplot as plt


def logprob(x):
    return np.sum(
        np.logaddexp(
            -0.5 * (x - 2) ** 2,
            -0.5 * (x + 2) ** 2,
        )
        - 0.5 * np.log(2 * np.pi)
        - np.log(2)
    )


x = np.linspace(-5.5, 5.5, 5000)
plt.plot(x, np.exp(list(map(logprob, x))), "k")
plt.yticks([])
plt.xlim(-5.5, 5.5)
plt.ylabel("p(x)")
plt.xlabel("x");

Now we can sample this using emcee and the default moves.StretchMove:

import emcee

np.random.seed(589403)

init = np.random.randn(32, 1)
nwalkers, ndim = init.shape

sampler0 = emcee.EnsembleSampler(nwalkers, ndim, logprob)
sampler0.run_mcmc(init, 5000)

print("Autocorrelation time: {0:.2f} steps".format(sampler0.get_autocorr_time()[0]))

Autocorrelation time: 40.03 steps

This autocorrelation time seems long for a 1D problem! We can also see this effect qualitatively by looking at the trace for one of the walkers:

plt.plot(sampler0.get_chain()[:, 0, 0], "k", lw=0.5)
plt.xlim(0, 5000)
plt.ylim(-5.5, 5.5)
plt.title("move: StretchMove", fontsize=14)
plt.xlabel("step number")
plt.ylabel("x");

For “lightly” multimodal problems like these, some combination of the moves.DEMove and moves.DESnookerMove can often perform better than the default. In this case, let’s use a weighted mixture of the two moves. In deatil, this means that, at each step, we’ll randomly select either a :class:moves.DEMove (with 80% probability) or a moves.DESnookerMove (with 20% probability).

np.random.seed(93284)

sampler = emcee.EnsembleSampler(
    nwalkers,
    ndim,
    logprob,
    moves=[
        (emcee.moves.DEMove(), 0.8),
        (emcee.moves.DESnookerMove(), 0.2),
    ],
)
sampler.run_mcmc(init, 5000)

print("Autocorrelation time: {0:.2f} steps".format(sampler.get_autocorr_time()[0]))

plt.plot(sampler.get_chain()[:, 0, 0], "k", lw=0.5)
plt.xlim(0, 5000)
plt.ylim(-5.5, 5.5)
plt.title("move: [(DEMove, 0.8), (DESnookerMove, 0.2)]", fontsize=14)
plt.xlabel("step number")
plt.ylabel("x");

Autocorrelation time: 6.49 steps

That looks a lot better!

The idea with the Moves interface is that it should be easy for users to try out several different moves to find the combination that works best for their problem so you should head over to Moves to see all the details!