Nested sampling algorithm

Background

Bayes' theorem can be used for model selection, where one has a pair of competing models ${\displaystyle M_1}$ and ${\displaystyle M_2}$ for data ${\displaystyle 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 ${\displaystyle M_1}$ may be calculated as:

${\displaystyle \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 ${\displaystyle M_1}$ and ${\displaystyle M_2}$ are already known, as they are chosen by the researcher ahead of time. However, the remaining Bayes factor ${\displaystyle 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, ${\displaystyle M_1}$ has a set of parameters that can be grouped together and called ${\displaystyle \theta }$, and ${\displaystyle M_2}$ has its own vector of parameters that may be of different dimensionality, but is still termed ${\displaystyle \theta }$. The marginalization for ${\displaystyle M_1}$ is

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

and likewise for ${\displaystyle 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 ${\displaystyle 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 ${\displaystyle Z=P(D\mid M)}$ where ${\displaystyle M}$ is ${\displaystyle M_1}$ or ${\displaystyle M_2}$:

Start with <span class="mwe-math-element mwe-math-element-inline" typeof="mw:Extension/math" about="#mwt37" data-mw='{"name":"math","attrs":{},"body":{"extsrc":"N"}}' id="mwOw"><span class="mwe-math-mathml-a11y" style="display: none;"><math xmlns="http://www.w3.org/1998/Math/MathML" class="mathjax_ignore" alttext="{\displaystyle N}"> <semantics> <mrow class="MJX-TeXAtom-ORD"> <mstyle displaystyle="true" scriptlevel="0"> <mi>N</mi> </mstyle> </mrow> <annotation encoding="application/x-tex">{\displaystyle N}</annotation> </semantics> </math></span><img src="https://wikimedia.org/api/rest_v1/media/math/render/svg/f5e3890c981ae85503089652feb48b191b57aae3" class="mwe-math-fallback-image-inline mw-invert skin-invert" aria-hidden="true" style="vertical-align: -0.338ex; width:2.064ex; height:2.176ex;" alt="{\displaystyle N}"/></span> points <span class="mwe-math-element mwe-math-element-inline" typeof="mw:Extension/math" about="#mwt38" data-mw='{"name":"math","attrs":{},"body":{"extsrc":"\\theta_1,\\ldots,\\theta_N"}}' id="mwPA"><span class="mwe-math-mathml-a11y" style="display: none;"><math xmlns="http://www.w3.org/1998/Math/MathML" class="mathjax_ignore" alttext="{\displaystyle \theta _{1},\ldots ,\theta _{N}}"> <semantics> <mrow class="MJX-TeXAtom-ORD"> <mstyle displaystyle="true" scriptlevel="0"> <msub> <mi>θ<!-- θ --></mi> <mrow class="MJX-TeXAtom-ORD"> <mn>1</mn> </mrow> </msub> <mo>,</mo> <mo>…<!-- … --></mo> <mo>,</mo> <msub> <mi>θ<!-- θ --></mi> <mrow class="MJX-TeXAtom-ORD"> <mi>N</mi> </mrow> </msub> </mstyle> </mrow> <annotation encoding="application/x-tex">{\displaystyle \theta _{1},\ldots ,\theta _{N}}</annotation> </semantics> </math></span><img src="https://wikimedia.org/api/rest_v1/media/math/render/svg/81cfb473ad54c2361a42e197d26edc1e866aab5d" class="mwe-math-fallback-image-inline mw-invert skin-invert" aria-hidden="true" style="vertical-align: -0.671ex; width:10.105ex; height:2.509ex;" alt="{\displaystyle \theta _{1},\ldots ,\theta _{N}}"/></span> sampled from prior. <b id="mwPQ">for</b> <span class="mwe-math-element mwe-math-element-inline" typeof="mw:Extension/math" about="#mwt39" data-mw='{"name":"math","attrs":{},"body":{"extsrc":"i=1"}}' id="mwPg"><span class="mwe-math-mathml-a11y" style="display: none;"><math xmlns="http://www.w3.org/1998/Math/MathML" class="mathjax_ignore" alttext="{\displaystyle i=1}"> <semantics> <mrow class="MJX-TeXAtom-ORD"> <mstyle displaystyle="true" scriptlevel="0"> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mstyle> </mrow> <annotation encoding="application/x-tex">{\displaystyle i=1}</annotation> </semantics> </math></span><img src="https://wikimedia.org/api/rest_v1/media/math/render/svg/9fb857369fbc28dab4b5d76bb835e6f196b1d1cf" class="mwe-math-fallback-image-inline mw-invert skin-invert" aria-hidden="true" style="vertical-align: -0.338ex; width:5.063ex; height:2.176ex;" alt="{\displaystyle i=1}"/></span> to <span class="mwe-math-element mwe-math-element-inline" typeof="mw:Extension/math" about="#mwt40" data-mw='{"name":"math","attrs":{},"body":{"extsrc":"j"}}' id="mwPw"><span class="mwe-math-mathml-a11y" style="display: none;"><math xmlns="http://www.w3.org/1998/Math/MathML" class="mathjax_ignore" alttext="{\displaystyle j}"> <semantics> <mrow class="MJX-TeXAtom-ORD"> <mstyle displaystyle="true" scriptlevel="0"> <mi>j</mi> </mstyle> </mrow> <annotation encoding="application/x-tex">{\displaystyle j}</annotation> </semantics> </math></span><img src="https://wikimedia.org/api/rest_v1/media/math/render/svg/2f461e54f5c093e92a55547b9764291390f0b5d0" class="mwe-math-fallback-image-inline mw-invert skin-invert" aria-hidden="true" style="vertical-align: -0.671ex; margin-left: -0.027ex; width:0.985ex; height:2.509ex;" alt="{\displaystyle j}"/></span> <b id="mwQA">do</b> % The number of iterations j is chosen by guesswork. <span class="mwe-math-element mwe-math-element-inline" typeof="mw:Extension/math" about="#mwt41" data-mw='{"name":"math","attrs":{},"body":{"extsrc":"L_i := \\min("}}' id="mwQQ"><span class="mwe-math-mathml-a11y" style="display: none;"><math xmlns="http://www.w3.org/1998/Math/MathML" class="mathjax_ignore" alttext="{\displaystyle L_{i}:=\min(}"> <semantics> <mrow class="MJX-TeXAtom-ORD"> <mstyle displaystyle="true" scriptlevel="0"> <msub> <mi>L</mi> <mrow class="MJX-TeXAtom-ORD"> <mi>i</mi> </mrow> </msub> <mo>:=</mo> <mo movablelimits="true" form="prefix">min</mo> <mo stretchy="false">(</mo> </mstyle> </mrow> <annotation encoding="application/x-tex">{\displaystyle L_{i}:=\min(}</annotation> </semantics> </math></span><img src="https://wikimedia.org/api/rest_v1/media/math/render/svg/625fb2d72ad0be182a7e2fc3c74aecf178c1ca88" class="mwe-math-fallback-image-inline mw-invert skin-invert" aria-hidden="true" style="vertical-align: -0.838ex; width:10.908ex; height:2.843ex;" alt="{\displaystyle L_{i}:=\min(}"/></span>current likelihood values of the points<span class="mwe-math-element mwe-math-element-inline" typeof="mw:Extension/math" about="#mwt42" data-mw='{"name":"math","attrs":{},"body":{"extsrc":")"}}' id="mwQg"><span class="mwe-math-mathml-a11y" style="display: none;"><math xmlns="http://www.w3.org/1998/Math/MathML" class="mathjax_ignore" alttext="{\displaystyle )}"> <semantics> <mrow class="MJX-TeXAtom-ORD"> <mstyle displaystyle="true" scriptlevel="0"> <mo stretchy="false">)</mo> </mstyle> </mrow> <annotation encoding="application/x-tex">{\displaystyle )}</annotation> </semantics> </math></span><img src="https://wikimedia.org/api/rest_v1/media/math/render/svg/775c8f99fbc2db4ef20dd618a468f110bae7bd76" class="mwe-math-fallback-image-inline mw-invert skin-invert" aria-hidden="true" style="vertical-align: -0.838ex; width:0.905ex; height:2.843ex;" alt="{\displaystyle )}"/></span>; <span class="mwe-math-element mwe-math-element-inline" typeof="mw:Extension/math" about="#mwt43" data-mw='{"name":"math","attrs":{},"body":{"extsrc":"X_i := \\exp(-i/N);"}}' id="mwQw"><span class="mwe-math-mathml-a11y" style="display: none;"><math xmlns="http://www.w3.org/1998/Math/MathML" class="mathjax_ignore" alttext="{\displaystyle X_{i}:=\exp(-i/N);}"> <semantics> <mrow class="MJX-TeXAtom-ORD"> <mstyle displaystyle="true" scriptlevel="0"> <msub> <mi>X</mi> <mrow class="MJX-TeXAtom-ORD"> <mi>i</mi> </mrow> </msub> <mo>:=</mo> <mi>exp</mi> <mo>⁡<!-- ⁡ --></mo> <mo stretchy="false">(</mo> <mo>−<!-- − --></mo> <mi>i</mi> <mrow class="MJX-TeXAtom-ORD"> <mo>/</mo> </mrow> <mi>N</mi> <mo stretchy="false">)</mo> <mo>;</mo> </mstyle> </mrow> <annotation encoding="application/x-tex">{\displaystyle X_{i}:=\exp(-i/N);}</annotation> </semantics> </math></span><img src="https://wikimedia.org/api/rest_v1/media/math/render/svg/acd186327738278ba0b89813944a7aaa12f2a5a1" class="mwe-math-fallback-image-inline mw-invert skin-invert" aria-hidden="true" style="vertical-align: -0.838ex; width:18.314ex; height:2.843ex;" alt="{\displaystyle X_{i}:=\exp(-i/N);}"/></span> <span class="mwe-math-element mwe-math-element-inline" typeof="mw:Extension/math" about="#mwt44" data-mw='{"name":"math","attrs":{},"body":{"extsrc":"w_i := X_{i-1} - X_i"}}' id="mwRA"><span class="mwe-math-mathml-a11y" style="display: none;"><math xmlns="http://www.w3.org/1998/Math/MathML" class="mathjax_ignore" alttext="{\displaystyle w_{i}:=X_{i-1}-X_{i}}"> <semantics> <mrow class="MJX-TeXAtom-ORD"> <mstyle displaystyle="true" scriptlevel="0"> <msub> <mi>w</mi> <mrow class="MJX-TeXAtom-ORD"> <mi>i</mi> </mrow> </msub> <mo>:=</mo> <msub> <mi>X</mi> <mrow class="MJX-TeXAtom-ORD"> <mi>i</mi> <mo>−<!-- − --></mo> <mn>1</mn> </mrow> </msub> <mo>−<!-- − --></mo> <msub> <mi>X</mi> <mrow class="MJX-TeXAtom-ORD"> <mi>i</mi> </mrow> </msub> </mstyle> </mrow> <annotation encoding="application/x-tex">{\displaystyle w_{i}:=X_{i-1}-X_{i}}</annotation> </semantics> </math></span><img src="https://wikimedia.org/api/rest_v1/media/math/render/svg/89745af0b7fb350531aff0aedc09a9ba7eebf3d0" class="mwe-math-fallback-image-inline mw-invert skin-invert" aria-hidden="true" style="vertical-align: -0.671ex; width:16.598ex; height:2.509ex;" alt="{\displaystyle w_{i}:=X_{i-1}-X_{i}}"/></span> <span class="mwe-math-element mwe-math-element-inline" typeof="mw:Extension/math" about="#mwt45" data-mw='{"name":"math","attrs":{},"body":{"extsrc":"Z := Z + L_i\\cdot w_i;"}}' id="mwRQ"><span class="mwe-math-mathml-a11y" style="display: none;"><math xmlns="http://www.w3.org/1998/Math/MathML" class="mathjax_ignore" alttext="{\displaystyle Z:=Z+L_{i}\cdot w_{i};}"> <semantics> <mrow class="MJX-TeXAtom-ORD"> <mstyle displaystyle="true" scriptlevel="0"> <mi>Z</mi> <mo>:=</mo> <mi>Z</mi> <mo>+</mo> <msub> <mi>L</mi> <mrow class="MJX-TeXAtom-ORD"> <mi>i</mi> </mrow> </msub> <mo>⋅<!-- ⋅ --></mo> <msub> <mi>w</mi> <mrow class="MJX-TeXAtom-ORD"> <mi>i</mi> </mrow> </msub> <mo>;</mo> </mstyle> </mrow> <annotation encoding="application/x-tex">{\displaystyle Z:=Z+L_{i}\cdot w_{i};}</annotation> </semantics> </math></span><img src="https://wikimedia.org/api/rest_v1/media/math/render/svg/f3a0c74d2884e1967ad44277cb3e7ce704bb9a7a" class="mwe-math-fallback-image-inline mw-invert skin-invert" aria-hidden="true" style="vertical-align: -0.671ex; width:17.119ex; height:2.509ex;" alt="{\displaystyle Z:=Z+L_{i}\cdot w_{i};}"/></span> Save the point with least likelihood as a sample point with weight <span class="mwe-math-element mwe-math-element-inline" typeof="mw:Extension/math" about="#mwt46" data-mw='{"name":"math","attrs":{},"body":{"extsrc":"w_i"}}' id="mwRg"><span class="mwe-math-mathml-a11y" style="display: none;"><math xmlns="http://www.w3.org/1998/Math/MathML" class="mathjax_ignore" alttext="{\displaystyle w_{i}}"> <semantics> <mrow class="MJX-TeXAtom-ORD"> <mstyle displaystyle="true" scriptlevel="0"> <msub> <mi>w</mi> <mrow class="MJX-TeXAtom-ORD"> <mi>i</mi> </mrow> </msub> </mstyle> </mrow> <annotation encoding="application/x-tex">{\displaystyle w_{i}}</annotation> </semantics> </math></span><img src="https://wikimedia.org/api/rest_v1/media/math/render/svg/fe22f0329d3ecb2e1880d44d191aba0e5475db68" class="mwe-math-fallback-image-inline mw-invert skin-invert" aria-hidden="true" style="vertical-align: -0.671ex; width:2.464ex; height:2.009ex;" alt="{\displaystyle w_{i}}"/></span>. Update the point with least likelihood by sampling from the prior restricted to likelihoods above <span class="mwe-math-element mwe-math-element-inline" typeof="mw:Extension/math" about="#mwt47" data-mw='{"name":"math","attrs":{},"body":{"extsrc":"L_i"}}' id="mwRw"><span class="mwe-math-mathml-a11y" style="display: none;"><math xmlns="http://www.w3.org/1998/Math/MathML" class="mathjax_ignore" alttext="{\displaystyle L_{i}}"> <semantics> <mrow class="MJX-TeXAtom-ORD"> <mstyle displaystyle="true" scriptlevel="0"> <msub> <mi>L</mi> <mrow class="MJX-TeXAtom-ORD"> <mi>i</mi> </mrow> </msub> </mstyle> </mrow> <annotation encoding="application/x-tex">{\displaystyle L_{i}}</annotation> </semantics> </math></span><img src="https://wikimedia.org/api/rest_v1/media/math/render/svg/fbbedf9f5f71ca52f8f24392c3cc18fca6c420dc" class="mwe-math-fallback-image-inline mw-invert skin-invert" aria-hidden="true" style="vertical-align: -0.671ex; width:2.383ex; height:2.509ex;" alt="{\displaystyle L_{i}}"/></span>, for example with <a rel="mw:WikiLink" href="./Markov_chain_Monte_Carlo" title="Markov chain Monte Carlo" id="mwSA">Markov chain Monte Carlo</a>. <b id="mwSQ">end</b> <b id="mwSg">return</b> <span class="mwe-math-element mwe-math-element-inline" typeof="mw:Extension/math" about="#mwt48" data-mw='{"name":"math","attrs":{},"body":{"extsrc":"Z"}}' id="mwSw"><span class="mwe-math-mathml-a11y" style="display: none;"><math xmlns="http://www.w3.org/1998/Math/MathML" class="mathjax_ignore" alttext="{\displaystyle Z}"> <semantics> <mrow class="MJX-TeXAtom-ORD"> <mstyle displaystyle="true" scriptlevel="0"> <mi>Z</mi> </mstyle> </mrow> <annotation encoding="application/x-tex">{\displaystyle Z}</annotation> </semantics> </math></span><img src="https://wikimedia.org/api/rest_v1/media/math/render/svg/1cc6b75e09a8aa3f04d8584b11db534f88fb56bd" class="mwe-math-fallback-image-inline mw-invert skin-invert" aria-hidden="true" style="vertical-align: -0.338ex; width:1.68ex; height:2.176ex;" alt="{\displaystyle Z}"/></span>;

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

${\displaystyle \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 ${\displaystyle j\to \mathrm{\infty }}$, this estimator has a positive bias of order ${\displaystyle 1/N}$ which can be removed by using ${\displaystyle (1-1/N)}$ instead of the ${\displaystyle \exp (-1/N)}$ in the above algorithm.

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