The Uses Of MCMC and Nested Sampling Algorithms

There are various off-the-shelf samplers that make use of MCMC and nested sampling algorithms in Python, freely available for the public to use. The following webpage is a collection of demonstrations of how a handful of popular samplers can be used to analyse real-world, open-source data sets.

A GitHub repository containing all the examples on this webpage can be found here
A Docker container containing all the samplers used on this webpage can be found here

A Brief Introduction to Bayesian Inference

For any experiment, you might have some kind of belief on what the outcome should be. This belief can be tested by recording the results of the experiment, and finding the likelihood that your beliefs explain the observed data. Using this you can update your prior beliefs using Bayes' rule:

$P(M|D) = \frac{P(D|M) P(M)}{P(D)}$

Where $P(M)$ is the prior probability belief of the model $M$ , $P(D|M)$ is the likelihood of observing data $D$ given the model, $P(M|D)$ is the posterior probability of the model given the data, and $P(D)$ is the "marginalised likelihood" or the "evidence", which is the probability of observing the data set.

Bayesian inference can be used in physics to fit models to observed data, given only a prior belief on the parameters of the model and a model for the noise observed in the data (the likelihood). This has uses in astronomy, uses in particle physics, and everywhere in between.

Markov Chain Monte Carlo

Evaluating the posterior over the prior space for multi-parameter models can be computationally difficult, or in many cases impossible for all practical purposes, so instead a method called Markov chain Monte Carlo (MCMC) is employed. An MCMC algorithm has two components: Monte Carlo techniques imply that samples are randomly generated, and Markov Chain implies that the next step of the algorithm depends only on the current step.

The most popular MCMC algorithm is the MH algorithm, which will completely ignore the $P(D)$ term. Instead, on the $i$ 'th step of the algorithm, the sampler will randomly generate a proposal model $M^*$ which may (or may not) do a better job at explaining the data than the current model $M^{i-1}$ . The ratio of posteriors

$r = \frac{P(D|M^*) P(M^*)}{P(D|M^{i-1}) P(M^{i-1})} = \frac{P(M^*|D)}{P(M^{i-1}|D)}$

is then calculated, using the prior and likelihood distributions. This ratio describes which model is most likely to explain the data, and so the proposed model $M^*$ will be accepted if this ratio is greater than 1. Otherwise, the acceptance probability is $P = r$ . Although this allows less likely models to be accepted, it means a wider parameter space can be explored such that a global likelihood maximum can be found, rather than a local maximum.

Once a proposed model is accepted, it is added to a set of posterior samples which is returned by the MCMC algorithm when all iterations are complete. These samples can be used to construct marginalised posterior distributions for the individual model parameters without having to analytically or numerically integrate the posterior distribution.

Nested Sampling

Nested sampling algorithms are a variation on MCMC, which use numerical methods to calculate the marginal likelihood. The name "marginal likelihood" is literal, meaning that the evidence can be found as the integral of the likelihood with respect to the model (marginalising the model out). When models have tens or hundreds of parameters, this marginalisation becomes impossible.

A very simple nested sampling algorithm essentially uses the trapezium rule to estimate the value of the marginalisation integral. This is especially useful for comparing models, since the marginalised likelihood is a metric for how well a model describes a data set. For two models $M_1$ and $M_2$ , the Bayes' factor

$K = \frac{P(D|M_1)}{P(D|M_2)}$

describes whether or not $M_1$ is a more likely model than $M_2$ given the data $D$ .

The Samplers

The following samplers are used in demonstrations on this site:

The Data

All the data used on this webpage is open source, and publicly available on the following links: