MCMC¶
This module implements a Markov Chain MonteCarlo 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 bestfit 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 bestfit 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 bestfit 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¶

Wrapper for the lowerlevel MCMC functions. 

Create a set of MCMC walkers. 

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

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 lowerlevel 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 onedimensional 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 withx_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 threedimensional numpy array of dimensions [
2, nwalkers/2, nparams]
. This array can be passed to therunChain()
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 2dimensional array with dimensions
[N, nparams]
where N is an arbitrary number, and the extra arguments given in theargs
tuple. walkers: array_like
A threedimensional numpy array with the initial coordinates of the walkers. This array must have dimensions
[2, nwalkers/2, nparams]
, and can be generated with theinitWalkers()
function. args: tuple
The extra arguments for the likelihood function.
 convergence_step: int
Save and output (if verbose) the GelmanRubin indicator, autocorrelation time etc. every
convergence_step
steps. converged_GR: float
The maximum difference between different chains, according to the GelmanRubin 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 autocorrelation time, meaning that the samples in this chain are truly independent. The chain can be analyzed with theanalyzeChain()
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.