Nested sampling algorithm

The nested sampling algorithm is a computational approach to the Bayesian statistics problems of comparing models and generating samples from posterior distributions. It was developed in 2004 by physicist John Skilling.

Background

Bayes' theorem can be used for model selection, where one has a pair of competing models $M_1$ and $M_2$ for data $D$, one of which may be true (though which one is unknown) but which both cannot be true simultaneously. Bayesian model selection provides a method for assessing the Bayes factor, which gives the relative merit of each model.

The posterior probability for $M_1$ may be calculated as:

$\begin{array}{cc}P(M_1\mid D)& =\frac{P(D\mid M_1)P(M_1)}{P(D)}\\ & =\frac{P(D\mid M_1)P(M_1)}{P(D\mid M_1)P(M_1)+P(D\mid M_2)P(M_2)}\\ & =\frac{1}{1+\frac{P(D\mid M_2)}{P(D\mid M_1)}\frac{P(M_2)}{P(M_1)}}\end{array}$

The prior probabilities $M_1$ and $M_2$ are already known, as they are chosen by the researcher ahead of time. However, the remaining Bayes factor $P(D\mid M_2)/P(D\mid M_1)$ is not so easy to evaluate, since in general it requires marginalizing nuisance parameters. Generally, $M_1$ has a set of parameters that can be grouped together and called $\theta$, and $M_2$ has its own vector of parameters that may be of different dimensionality, but is still termed $\theta$. The marginalization for $M_1$ is

$P(D\mid M_1)=\int d\theta P(D\mid \theta ,M_1)P(\theta \mid M_1)$

and likewise for $M_2$. This integral is often analytically intractable, and in these cases it is necessary to employ a numerical algorithm to find an approximation. The nested sampling algorithm was developed by John Skilling specifically to approximate these marginalization integrals, and it has the added benefit of generating samples from the posterior distribution $P(\theta \mid D,M_1)$. It is an alternative to methods from the Bayesian literature such as bridge sampling and defensive importance sampling.

Here is a simple version of the nested sampling algorithm, followed by a description of how it computes the marginal probability density $Z=P(D\mid M)$ where $M$ is $M_1$ or $M_2$:

Start with $N$ points ${\theta }_1,\dots ,{\theta }_N$ sampled from prior.
for $i=1$ to $j$ do        % The number of iterations j is chosen by guesswork.
    $L_i:=\min ($current likelihood values of the points$)$;
    $X_i:=\exp (-i/N);$
    $w_i:=X_{i-1}-X_i$
    $Z:=Z+L_i\cdot w_i;$
    Save the point with least likelihood as a sample point with weight $w_i$.
    Update the point with least likelihood by sampling from the prior restricted to likelihoods above $L_i$, for example with Markov chain Monte Carlo.
end
return $Z$;

At each iteration, $X_i$ is an estimate of the amount of prior mass covered by the hypervolume in parameter space of all points with likelihood greater than ${\theta }_i$. The weight factor $w_i$ is an estimate of the amount of prior mass that lies between two nested hypersurfaces $\{\theta \mid P(D\mid \theta ,M)=P(D\mid {\theta }_{i-1},M)\}$ and $\{\theta \mid P(D\mid \theta ,M)=P(D\mid {\theta }_i,M)\}$. The update step $Z:=Z+L_i{w}_i$ computes the sum over $i$ of $L_i{w}_i$ to numerically approximate the integral

$\begin{array}{cc}P(D\mid M)& =\int P(D\mid \theta ,M)P(\theta \mid M)d\theta \\ & =\int P(D\mid \theta ,M)dP(\theta \mid M)\end{array}$

In the limit $j\to \mathrm{\infty }$, this estimator has a positive bias of order $1/N$ which can be removed by using $(1-1/N)$ instead of the $\exp (-1/N)$ in the above algorithm.

The idea is to subdivide the range of $f(\theta )=P(D\mid \theta ,M)$ and estimate, for each interval $[f({\theta }_{i-1}),f({\theta }_i)]$, how likely it is a priori that a randomly chosen $\theta$ would map to this interval. This can be thought of as a Bayesian's way to numerically implement Lebesgue integration.

Likelihood-restricted prior sampling algorithms

The point with least likelihood can be updated with some Markov chain Monte Carlo steps according to the prior, accepting only steps that keep the likelihood above $L_i$. The original procedure outlined by Skilling (given above in pseudocode) does not specify what specific algorithm should be used to choose new points with better likelihood, but several algorithms have been developed.

Skilling's own code examples (such as one in Sivia and Skilling (2006), available on Skilling's website) chooses a random existing point and selects a nearby point chosen by a random distance from the existing point; if the likelihood is better, then the point is accepted, else it is rejected and the process repeated. Subsequently, a variety of MCMC algorithms tailored for nested sampling have been developed, including slice sampling, popularized by PolyChord, and constrained Hamiltonian Monte Carlo.

An alternative line of algorithms is based on rejection sampling. Mukherjee et al. (2006) found higher acceptance rates by selecting points randomly within an ellipsoid drawn around the existing points; this idea was refined by the MultiNest algorithm which handles multimodal posteriors better using multiple ellipsoids built from clustering of the live points. Rejection methods can be efficient up to 20-30 dimensions.

Implementations

Example implementations demonstrating the nested sampling algorithm are publicly available for download, written in several programming languages.

• Simple examples in C, R, or Python are on John Skilling's website.
• A Haskell port of the above simple codes is on Hackage.
• An example in R originally designed for fitting spectra is described on Bojan Nikolic's website and is available on GitHub.
• A NestedSampler is part of the Python toolbox BayesicFitting for generic model fitting and evidence calculation. It is available on GitHub.
• An implementation in C++, named DIAMONDS, is on GitHub.
• A highly modular Python parallel example for statistical physics and condensed matter physics uses is on GitHub.
• pymatnest is a package designed for exploring the energy landscape of different materials, calculating thermodynamic variables at arbitrary temperatures and locating phase transitions is on GitHub
• The MultiNest software package is capable of performing nested sampling on multi-modal posterior distributions. It has interfaces for C++, Fortran and Python inputs, and is available on GitHub.
• PolyChord is another nested sampling software package available on GitHub. PolyChord's computational efficiency scales better with an increase in the number of parameters than MultiNest, meaning PolyChord can be more efficient for high dimensional problems. It has interfaces to likelihood functions written in Python, Fortran, C, or C++. PolyChord can be used jointly with Cobaya, a Python-based code for Bayesian analysis of hierarchical physical models. Cobaya facilitates exploration of posteriors using various Monte Carlo samplers, allows for maximization and importance-reweighting of samples, and includes interfaces to cosmological theory codes and likelihoods.
• NestedSamplers.jl, a Julia package for implementing single- and multi-ellipsoidal nested sampling algorithms is on GitHub.
• Korali is a high-performance framework for uncertainty quantification, optimization, and deep reinforcement learning, which also implements nested sampling.
• The UltraNest software package implements a fast MPI-capable, generalized multi-ellipsoidal dynamic nested sampling algorithm. The user can also choose slice sampling algorithms. Written in Python, it has interfaces for Python, C, Fortran, Python, R and Julia and is available on GitHub.

Applications

Since nested sampling was proposed in 2004, it has been used in many science areas, in particular in astronomy. One paper suggested using nested sampling for cosmological model selection and object detection, as it "uniquely combines accuracy, general applicability and computational feasibility." A refinement of the algorithm to handle multimodal posteriors has been suggested as a means to detect astronomical objects in extant datasets. Other applications of nested sampling are in the field of finite element updating where the algorithm is used to choose an optimal finite element model, and this was applied to structural dynamics. This sampling method has also been used in the field of materials modeling. It can be used to learn the partition function from statistical mechanics and derive thermodynamic properties.

Diagnostics

Dedicated diagnostics for nested sampling have been developed for verifying that a nested sampling run is performing well. This includes a U test that the rank of the likelihood of the replacement point is uniformly distributed among the live points, Markov Chain Monte Carlo jump distance, comparison of the consistency of several independent nested sampling runs, including reruns with increased number of MCMC steps. The computation can also be checked with generically applicable techniques such as simulation-based calibration.

Dynamic nested sampling

Dynamic nested sampling is a generalisation of the nested sampling algorithm in which the number of samples taken in different regions of the parameter space is dynamically adjusted to maximise calculation accuracy. This can lead to improvements in accuracy and computational efficiency when compared to the original nested sampling algorithm, in which the allocation of samples cannot be changed and often many samples are taken in regions which have little effect on calculation accuracy.

Publicly available dynamic nested sampling software packages include:

• dynesty – a Python implementation of dynamic nested sampling which can be downloaded from GitHub.
• dyPolyChord: a software package which can be used with Python, C++ and Fortran likelihood and prior distributions. dyPolyChord is available on GitHub.
• UltraNest (see above).

Dynamic nested sampling has been applied to a variety of scientific problems, including analysis of gravitational waves, mapping distances in space and exoplanet detection.

See also

• Bayesian model comparison
• List of algorithms