MCMC

This module implements a Markov Chain Monte-Carlo sampler based on the Goodman & Weare 2010 algorithm. The code was originally written by Andrey Kravtsov and adapted for Colossus by Benedikt Diemer.

Basics

Say we want to find the best-fit parameters x of some model. The likelihood for this model should take on a form like this:

def likelihood(x, data, some_other_variables):

Then all we need to do to obtain the best-fit parameters is to make some initial guess for x that is stored in the vector x_initial, and run:

from colossus.utils import mcmc

args = data, some_other_variables
mcmc.run(x_initial, likelihood, args = args)

The run() function is a wrapper around the main stages of the MCMC sampling process. If, for example, we wish to obtain the mean best-fit parameters, we can execute the following code:

args = data, some_other_variables
walkers = mcmc.initWalkers(x_initial)
chain_thin, chain_full, _ = mcmc.runChain(likelihood, walkers, args = args)
mean, _, _, _ = mcmc.analyzeChain(chain_thin, param_names = param_names)

There are numerous more advanced parameters as listed below. Please see the Tutorials for more code examples, including a function to plot MCMC chains generated with this module.

Module contents

run(x_initial, L_func[, args, verbose, ...])	Wrapper for the lower-level MCMC functions.
initWalkers(x_initial[, initial_step, ...])	Create a set of MCMC walkers.
runChain(L_func, walkers[, args, ...])	Run an MCMC chain using the Goodman & Weare (2010) algorithm.
analyzeChain(chain[, param_names, ...])	Analyze an MCMC chain.

Module reference

utils.mcmc.run(x_initial, L_func, args=(), verbose=True, initial_step=0.1, nwalkers=100, random_seed=None, convergence_step=100, converged_GR=0.01, param_names=None, percentiles=[68.27, 95.45, 99.73])

Wrapper for the lower-level MCMC functions.

This function combines the creation of walkers, running of the chain and analysis functions into one. Please see the documentation of the individual functions for the meaning of the input parameters.

utils.mcmc.initWalkers(x_initial, initial_step=0.1, nwalkers=100, random_seed=None)

Create a set of MCMC walkers.

This function distributes the initial positions of the walkers in an isotropic Gaussian around the initial guess provided by the user. The output of this function serves as an input for the runChain() function.

Parameters

x_initial: array_like

A one-dimensional numpy array. The length of this array is the number of parameters varied in the MCMC run.

initial_step: array_like

The width of the Gaussian in parameter i. Can either be a numpy array of the same length as x_initial, or a number. In the latter case, the number is multiplied with x_initial, i.e. the same fractional width of the Gaussian is applied in all parameter dimensions.

nwalkers: int

The number of walkers created. In this implementation, the walkers are split into two subgroups, meaning their number must be divisible by two.

random_seed: int

If not None, this random seed is used when generated the walker positions.

Returns

walkers: array_like

A three-dimensional numpy array of dimensions [2, nwalkers/2, nparams]. This array can be passed to the runChain() function.

utils.mcmc.runChain(L_func, walkers, args=(), convergence_step=100, converged_GR=0.01, verbose=True, output_every_n=100)

Run an MCMC chain using the Goodman & Weare (2010) algorithm.

Parameters

L_func: function

The likelihood function which is maximized. This function needs to accept two parameters: a 2-dimensional array with dimensions [N, nparams] where N is an arbitrary number, and the extra arguments given in the args tuple.

walkers: array_like

A three-dimensional numpy array with the initial coordinates of the walkers. This array must have dimensions [2, nwalkers/2, nparams], and can be generated with the initWalkers() function.

args: tuple

The extra arguments for the likelihood function.

convergence_step: int

Save and output (if verbose) the Gelman-Rubin indicator, autocorrelation time etc. every convergence_step steps.

converged_GR: float

The maximum difference between different chains, according to the Gelman-Rubin criterion. Once the GR indicator is lower than this number in all parameters, the chain is ended.

verbose: bool

If False, this function outputs no information.

output_every_n: int

Output information about the progress of the chain every n steps. This parameter only has an effect if verbose == True.

Returns

chain_thin: array_like

A numpy array of dimensions [n_independent_samples, nparams] with the parameters at each step in the chain. In this thin chain, only every nth step is output, where n is the auto-correlation time, meaning that the samples in this chain are truly independent. The chain can be analyzed with the analyzeChain() function.

chain_full: array_like

Like the thin chain, but including all steps. Thus, the samples in this chain are not indepedent from each other. However, the full chain often gives better plotting results.

R: array_like

A numpy array containing the GR indicator at each step when it was saved.

utils.mcmc.analyzeChain(chain, param_names=None, percentiles=[68.27, 95.45, 99.73], verbose=True)

Analyze an MCMC chain.

An MCMC chain represents a statistical sample of the likelihood in question. This function computes more convenient parameters such as the mean, median, and various percentiles for each of the parameters.

Parameters

chain: array_like

A numpy array of dimensions [nsteps, nparams] with the parameters at each step in the chain. The chain is created by the runChain() function.

param_names: array_like

Optional; a list of strings which are used when outputting the parameters.

percentiles: array_like

A list with percentages which are output. By default, the classical 1, 2, and 3 sigma probabilities are used.

verbose: bool

Print the results.

Returns

x_mean: array_like

The mean of the chain for each parameter; has length nparams.

x_median: array_like

The median of the chain for each parameter; has length nparams.

x_stddev: array_like

The standard deviation of the chain for each parameter; has length nparams.

x_percentiles: array_like

The lower and upper values of each parameter that contain a certain percentile of the probability; has dimensions [n_percentages, 2, nparams] where the second dimension contains the lower/upper values.