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.
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)
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.
||Wrapper for the lower-level MCMC functions.|
||Create a set of MCMC walkers.|
||Run an MCMC chain using the Goodman & Weare (2010) algorithm.|
||Analyze an MCMC chain.|
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.
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
- 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
None, this random seed is used when generated the walker positions.
- walkers: array_like
A three-dimensional numpy array of dimensions [
2, nwalkers/2, nparams]. This array can be passed to the
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.
- 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
- 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
- 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
- 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
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.
- 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
- 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.
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.
- 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
- 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.
- x_mean: array_like
The mean of the chain for each parameter; has length
- x_median: array_like
The median of the chain for each parameter; has length
- x_stddev: array_like
The standard deviation of the chain for each parameter; has length
- 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.