Nested sampling algorithm

Bayes' theorem can be applied to a pair of competing models $M1$ and $M2$ for data $D$, one of which must be true (though which one is not known) but which both cannot simultaneously be true, as follows:

${\begin{array}{lcl}P(M1|D)&=&P(D|M1)P(M1)/P(D)\\&=&P(D|M1)P(M1)/(P(D|M1)P(M1)+P(D|M2)P(M2)\\&=&1/(1+P(D|M2)P(M2)/P(D|M1)P(M1))\\&=&1/(1+(P(D|M2)/P(D|M1))(P(M2)/P(M1)))\\\end{array}}$

Given no a priori information in favor of $M1$ or $M2$, it is reasonable to assign prior probabilities $P(M1)=P(M2)=1/2$, so that $P(M2)/P(M1)=1$. The remaining ratio $P(D|M2)/P(D|M1)$ is not so easy to evaluate since in general it requires marginalization of nuisance parameters. Generally, $M1$ has a collection of parameters that can be lumped together and called $\theta$, and $M2$ has its own vector of parameters that may be of different dimensionality but is still referred to as $\theta$. The marginalization for $M1$ is

$P(D|M1)=\int d\theta P(D|\theta ,M1)P(\theta |M1)$

and likewise for $M2$. This integral is often analytically intractible, 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 |D,M1)$ [1].

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|M)$ where $M$ is $M1$ or $M2$:

  Start with $N$ points $\theta _{1},...,\theta _{N}$ sampled from prior.

  for $i=1$ to $j$ do        % The number of iterations j is chosen by guesswork.
      $Li:=\min($current likelihood values of the points$)$;
      $Xi:=\exp(-i/N);$
      $wi:=X_{i-1}-X_{i}$
      $Z:=Z+Li*wi;$

      Save the point with least likelihood as a sample point with weight wi.

      Update the point with least likelihood with some Markov Chain
      Monte Carlo steps according to the prior, accepting only steps that
      keep the likelihood above $Li$.
  end
  return $Z$;

At each iteration, $Xi$ 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 $wi$ is an estimate of the amount of prior mass that lies between two nested hypersurfaces $\{\theta |P(D|\theta ,M)=P(D|\theta _{i-1},M)\}$ and $\{\theta |P(D|\theta ,M)=P(D|\theta _{i},M)\}$. The update step $Z:=Z+Li*wi$ computes the sum over $i$ of $Li*wi$ to numerically approximate the integral

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

The idea is to chop up the range of $f(\theta )=P(D|\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.

References

[1] Nested Sampling for General Bayesian Computation. John Skilling. 2006. International Society for Bayesian Analysis.