reddemcee

Adaptive Parallel Tempering MCMC Ensemble Sampler, made for the exoplanet finder algorithm EMPEROR. This sampler works as a stand-alone program, so the community might find it useful.

Overview

reddemcee is an Adaptive Parallel Tempering MCMC implementation based on the excellent emcee code, and a modified version of the Vousden et al. implementation.

It's coded in such a way that minimal differences in input are required compared to emcee (v. 3.1.3). Make sure to check reddemcee's documentation, as this readme might not be updated.

Dependencies

This code makes use of:

• numpy
• tqdm
• emcee
• scipy

tqdm is used to display progress bars on the terminal. emcee is used for calculating the autocorrelation times.

Installation

Pip

In the console type

pip install reddemcee

From Source

In the console type

git clone https://github.com/ReddTea/reddemcee
cd reddemcee
python -m pip install -e .

Usage

For a complete tutorial please refer to the documentation's tutorial or the test file in the tests folder.

import numpy as np
import reddemcee

def log_like(x, ivar):
    return -0.5 * np.sum(ivar * x ** 2)

def log_prior(x):
    return 0.0

ndim = 2
ntemps, nwalkers, nsweeps, nsteps = 5, 50, 100, 2

ivar = 1. / np.random.rand(ndim)
p0 = np.random.randn(10, nwalkers, ndim)
sampler = reddemcee.PTSampler(nwalkers,
                             ndim,
                             log_like,
                             log_prior,
                             ntemps=ntemps,
                             loglargs=[ivar],
                             )
                             
sampler.run_mcmc(p0, nsweeps, nsteps)  # starting coords, nsweeps, nsteps

Configuring the Sampler

When setting up PTSampler, you can use the arguments:

The adaptation modes try to equalise the following quantity:

Additional Functions

Additional functions on the sampler include:

All these functions accept as arguments:

For example, the previous chain would have shape (5, 200, 50, 2), for 5 temperatures, 200 steps (nsweeps*nsteps), 50 walkers, and 2 dimensions.

sampler.get_chain(discard=100)

Would return the samples with shape (5, 100, 50, 2). Dropping the first 100 for every walker.

sampler.get_chain(discard=100, flat=True)

Would return the samples with shape (5, 5000, 2), 'linearising' the walkers.