Nested Sampling for physical scientists

Since NS is primarily an algorithm for integration, let us write a general multi-dimensional integral of a function $L$ over parameters $\Theta$ as

In many scientific problems we need to be able to integrate in high dimensions and for challenging integrands. We assume that the integrand is positive, $L(\Theta)\geq 0$.

Often, $Z$ may be a physical quantity such as the total mass of an object distributed with density $L$ across volumes $\text{d}\mu(\Theta)$. Whilst NS is a general method for integration, for concreteness, we view all such applications through the lens of Bayesian inference (see ), with $\text{d}\mu(\Theta)\equiv\pi(\Theta)\text{d}\Theta$ seen as an element of prior probability with $\pi$ the prior , normalized by its nature to $\int\pi(\Theta)\text{d}\Theta=1$. The integrand, $L$, is the modulating likelihood function (hence the symbol) and $Z$ is the evidence . In scientific inference problems, the integral could be over tens if not hundreds of parameters, required to model fundamental effects as well as the calibration and systematics of complicated experimental measurements [13].

We may rewrite the integrand through the elementary factorization known as Bayes’ theorem,

where

is the posterior , normalized to $\int P(\Theta)\,\text{d}\Theta=1$. Notation apart, however, all we are doing here is decomposing the integrand into a magnitude $Z$ and a shape $P(\Theta)$.

Historically, Bayesian computation focused on only the shape $P(\Theta)$, partly owing to controversies around Bayesian model comparison such as its sensitivity to the choices of prior [12], and partly due to computational difficulties [14]. However, shapes and magnitudes both matter, especially in the general setting of multi-dimensional integration beyond Bayesian model comparison. NS [1, 2] surmounts the challenge by computing shapes and magnitudes simultaneously.

Before introducing NS, let us attempt to simplify the general integral in eq. 1. Consider traditional Riemann-style integration. This decomposes the space into volume elements $\Delta\Theta$, typically small cubes, and performs a sum over them. Small cubes, however, rapidly become infeasible in multi-dimensional integration because their cost grows exponentially with dimension — this is the “curse of dimensionality.”

We don’t, however, need to decompose our space into little cubes; our cells can be any shape we want. The integrals needed for quantification,

are defined as limiting sums over volume elements which should be small enough to keep $P(\Theta)$ almost constant regardless of shape. We may therefore combine the cells in which the integrand is almost constant. Schematically, we may write

where $\Delta X$ is the volume of cells that share likelihood $L(X)$ weighted by the prior $\pi(\Theta)$. This is illustrated schematically in fig. 1b and works whether the integrand is uni- or multi-modal.

[box:math]Mathematical details The idea of NS is to transform eq. 1 into eq. 8 which can be approximated more efficiently using the described Monte Carlo approach with active samples:

where $\Omega$ is the parameter domain and $\tilde{L}$ is an overloaded form of $L$, as described next. To achieve this transformation, we define the survival function $X:\mathbb{R}\to[0,1]$, $X(\lambda)=\mu(\{z\in\Omega:L(z)>\lambda\})$, that is the $\mu$-measure of the $\lambda$-super-level-sets. Then, we introduce a mapping $\Phi:\Omega\to[0,1]$ with

This mapping $\Phi$ is the transformation which allows us to shift the integration from $\Omega$ to $[0,1]$, by virtue of the push-forward measure $\mu\circ\Phi^{-1}$ on $[0,1]$. Now, we can define

Thus the integral transformation above which provides an alternative characterization of $Z$ is true because

Here, $\tilde{L}$ is a generalized inverse of $X$ under suitable assumptions on $L$, see ref. [15, 16]. This integral transformation holds at least in the case that $L$ has no plateaus of positive prior measure, owing to the fact that in this case $\mu\circ\Phi^{-1}$ is indeed the uniform measure on $[0,1]$, that is $\text{d}\mu\circ\Phi^{-1}(x)=\text{d}x$.

If $L$ has a plateau of non-negligible mass, i.e., there exists a level $\lambda^{\star}$ such that $\mu(\{z\in\Omega:L(z)=\lambda^{\star}\})>0$, then the derivation is more challenging, see also 6 Limitations and optimizations.

We can reach eq. 5 more concretely by noting that the evidence is the expectation of a non-negative random variable, such that it may be written as

where the volume variable $X$,

is the volume enclosed by contour $L^{\star}$. This result can be readily proven by integration by parts (see also sec. 21 in ref. [17]). Applying integration by parts again to eq. 6, we obtain the familiar NS evidence identity,

providing that $L(X)$, the inverse of $X(\lambda)$, indeed exists and that the evidence is finite. This formalizes the schematic eq. 5. We discuss this result more formally in .

As the multi-dimensional integral in eq. 7 is impractical in high dimension, some sort of statistical estimation is inevitable. NS starts with an ensemble of $n_{\text{live}}$ random locations $\Theta$ drawn from the prior, $\pi(\Theta)$, each of which has its likelihood $L(\Theta)$ which we can place in ascending order. Crudely, if we discarded the lowest half of the values, the survivors would be random samples taken within the restricted volume $L>\text{Median}[L]$, which would statistically be roughly half the original volume. This allows us to make a statistical estimate of the volume variable in eq. 7. Repeating that $n_{\text{iter}}$ times would yield compression by a factor of about $2^{n_{\text{iter}}}$. This is the exponential behaviour required to overcome the curse of dimensionality.

Thus NS works by statistical estimates of the compression, which is a general and fundamental operation that can be used in various ways not limited to those in table 1. The evidence identity in eq. 6 isn’t required in every application; many applications only use the compression in eq. 7. We present scientific applications in 4 Applications.

We now present the NS algorithm in more detail. We assume that there are no regions of constant likelihood resulting in likelihood plateaus (see 6 Limitations and optimizations for further discussion). The NS algorithm begins by drawing an ensemble of $n_{\text{live}}$ samples from the prior. We compute the likelihood for each sample. We denote the smallest likelihood by $L^{\star}$ and we discard that point. The remaining live points are now distributed over a compressed volume; we denote the factor by which the volume compressed by $t$. Finally, a replacement point is drawn from the prior subject to $L>L^{\star}$, that is, from the constrained prior,

This leaves a new ensemble with $n_{\text{live}}$ samples obeying a likelihood constraint $L>L^{\star}$.

As they are drawn from the constrained prior, the volumes $X$ associated with the live points are uniformly distributed. Thus the compression associated with the discarded outermost sample, $t$, corresponds to the smallest of $n_{\text{live}}$ uniform random variables. This follows a $\text{Beta}(n_{\text{live}},1)$ distribution,

The first factor accounts for the fact that any live point could be the outermost and the second factor for the fact that $n_{\text{live}}-1$ uniformly distributed samples lie above the outermost sample at $t$ (the remaining sample lies at $t$).

We started from unconstrained samples drawn from the full original volume $X_{0}=1$. As we repeat this process of replacing the outermost points, the successive compressions by factors $t_{1},t_{2},t_{3},\dots$ lead to exponentially decreasing inferred volumes $X_{1},X_{2},X_{3},\dots$

We illustrate a single iteration of NS in fig. 1c. If we are most interested in the magnitude of the evidence, $\log Z$, we should consider,

as under repeated multiplication $X_{k}=t_{1}t_{2}\dots t_{k}$ it’s the logarithms that add. For a more complete inference, the compression factors $t$ may be sampled directly from $\text{Beta}(n_{\text{live}},1)$. See appendix A for further discussion. NS therefore estimates volumes through probability not geometry, topology, or even dimension. We do not get a definite compression value, only a distribution of what it might have been.

Having obtained estimates of the volume $X(L^{\star})$ at each of $n_{\text{iter}}$ iterations, we may accumulate the evidence via eq. 6 or eq. 8, for example by the trapezium rule,

where the weights are

The sum in eq. 13 converges to the desired integral [18, 19, 20]. This is the magnitude; we obtain shape by assigning weights to each sample (see ref. [18] for discussion),

normalized such that $\sum P_{i}=1$. These are the posterior weights of the dead points; the shape may be recovered by for example a weighted histogram or other density estimation methods. We show the whole algorithm schematically in algorithm 1 and the summation for a two-dimensional problem in fig. 1b.

Compared to the sketch of NS in section 1.3, we replace a single live point per iteration, rather than half of the live points. Whilst the number of replacements can be varied, one replacement is optimal (though see considerations in 2.4 Parallelization). This is because for $r\ll n_{\text{live}}$, replacing $r$ points per iteration would reach the posterior bulk in about $r$ times fewer iterations. While the computational expense wouldn’t change as $r$ replacements are required per iteration, the error estimates would scale as $\sqrt{r}$ because reducing the number of iterations increases the relative Poisson noise in the number of iterations.

In the last decade or so, analogies between NS and statistical mechanics (see ) and other statistical methods, including sequential Monte Carlo (SMC; [21]) and subset simulation in rare event sampling [24, 25, 22, 23], among others [26, 27, 28, 29, 30], have been recognized. The connections to SMC and an SMC variant of NS are discussed in . There are, of course, other strategies for computing the evidence; see ref. [31, 32, 33, 34] for reviews. Notable examples include approximating the integrand by a tractable function[35], Chib’s method using density estimation and Gibbs’ sampling[36], importance sampling[37], and techniques that re-use MCMC draws[38]. Broadly speaking, NS lies in a class of algorithms that form a path of bridging distributions, and evolves samples along that path [39, 40]. NS stands out because the path is automatic and smooth — we compress in $\log X$ by on average $1/n_{\text{live}}$ at each iteration — and because along the path we compress through constrained priors, rather than from the prior to the posterior. This was in fact a motivation for NS as it avoids phase transitions — abrupt changes in the bridging distributions — that cause problems for other methods including path samplers such as annealing. We further discuss NS’s historical background and contrast it with annealing in appendix B.

Our estimates of the magnitude and shape in eq. 1 must be accompanied by a discussion of the bias and statistical uncertainty. The latter originates from our noisy estimates of the compression factors. We may estimate the resulting statistical uncertainty in the evidence by considering the compression required to reach the bulk of the posterior. This may be quantified by the information content [41, 42]

known in statistics as the Kullback-Leibler (KL) divergence. We can write it using the volume variable as

where $P(X)\equiv L(X)/Z$ is the posterior density of the volume. In eq. 18 we see that the KL divergence equals minus the posterior expectation of $\log X$ minus the differential entropy associated with the posterior of $\log X$. As the first term typically dominates, the KL divergence provides a measure of compression.

Thus from eq. 12 it’s likely to take about

iterations to compress to the bulk of the posterior at $\log X=-H$, and this count is likely to be subject to $\sqrt{n_{H}}$ Poisson variability. Neglecting contributions from outside the bulk, we may write the evidence sum in eq. 13 as

We see that the final estimate of $\log Z$ will be plausibly subject to an approximately Gaussian uncertainty from the first factor,

Thus NS statistical uncertainties scale as $1/\sqrt{n_{\text{live}}}$ as usual for statistical uncertainties (see ref. [18] for an alternative proof and discussion).

Although NS was first introduced with this estimate [1], it can be unreliable. Authority rests with repeated simulation through eq. 10 of what the compressions might actually have been. See ref. [43] for further discussion of the statistical uncertainty in the NS estimates and appendix D for discussion of uncertainties in NS estimates of the posterior. As well as this statistical uncertainty, there are four potential sources of bias: bias originating from failure to faithfully sample from the constrained prior, bias originating from the choice of estimator for the compression factor, the generally negligible quadrature error, and the potentially important truncation error in eq. 13. The latter occurs as we stop after a finite number of iterations and is further discussed in 2.5 Stopping conditions. Provided that NS is appropriately configured, the statistical uncertainty usually dominates.

We see that difficulty in NS does not in fact lie in dimension but in compression from the prior to the posterior: the compression and the resolution $n_{\text{live}}$ alone determine the uncertainty and the run-time. To maintain a given uncertainty, by eq. 21 we require $n_{\text{live}}\propto H$ live points and by eq. 19 $n_{\text{iter}}\propto H^{2}$ iterations. If the prior and posterior are factorizable into a term for each dimension, by eq. 16 the KL divergence is additive, and so scales linearly with dimension, $D$, and so run-time goes like $\mathcal{O}\mathopen{}\mathclose{{\left(D^{2}}}\right)$. NS beats the exponential scaling with dimension expected from the curse of dimensionality.

[box:stat_mech]Statistical mechanics analogy There is a strong analogy between Bayesian inference and statistical mechanics. This suggests that NS might be useful in exploring problems that are typically the subject of statistical mechanical analysis. Consider $\Theta$ a microstate and $E(\Theta)=-\log L(\Theta)$ its energy. Then the microcanonical ensemble includes all states with $E(\Theta)=\epsilon$ where $\epsilon$ is a constant energy level. The volume of state space corresponding to a given energy $\epsilon$ is given by the density of states,

and taking the prior to be uniform corresponds to the ergodic hypothesis, namely that each microstate that is allowed by applicable conservation laws is equally likely to be observed. The Laplace transform of the density of states is the canonical partition function $Z(\beta)=\int e^{-\beta E}g(E)\text{d}E$, which corresponds to the generalized evidence (eq. 38). The canonical ensemble is an alternative description of thermodynamic states that is based on the inverse temperature $\beta$ rather than the energy level $\epsilon$. In the canonical ensemble, states follow the Boltzmann distribution $p(\Theta\mid\beta)=\exp\{-\beta E(\Theta)\}/Z(\beta)$. NS in essence tracks the cumulative density of states via:

The practical approximation corresponding to eqs. 13 and 14 is

During NS, states are generated from the prior constrained by an upper energy limit $\epsilon=-\log L^{\star}$, which can be achieved with any number of techniques, such as simple rejection sampling, Galilean Monte Carlo and Hamiltonian dynamics[44] or Creutz’ microcanonical demon algorithm [45, 46]. The information entropy of the constrained prior is the volume entropy $\log X(\epsilon)$ (also known as Gibbs entropy). As NS progresses from one energy limit $\epsilon$ to the next $\epsilon^{\prime}<\epsilon$, the volume entropy changes by $\Delta H=\log[X(\epsilon)/X(\epsilon^{\prime})]$ at a rate that is constant on average: $\langle\Delta H\rangle=1/n_{\text{live}}$.

Widely-used tempering methods such as simulated annealing work in the canonical ensemble and use the inverse temperature $\beta$ as a control parameter to weight each microstates. Alternatives include the Wang-Landau method[47] and NS, both of which use energy as a control parameter. In contrast to the ensemble property $\beta$, a key advantage of $E(\Theta)$ is that it can be evaluated for a single microstate. The Wang-Landau method uses a fixed, predefined set of energy bins, while NS constructs a sequence of energy levels at runtime. The sequence is optimal in that it achieves constant thermodynamic speed because changes in volume entropy are constant on average. Therefore, NS elegantly avoids a major problem of canonical annealing methods: designing a good temperature schedule.

As discussed in 1.4 Formulation, the number of live points controls the rate of exponential compression during an NS run; we compress inwards by about $\Delta\log X\approx 1/n_{\text{live}}$ per iteration. This means that run-time scales as $\mathcal{O}\mathopen{}\mathclose{{\left(n_{\text{live}}}}\right)$ (eq. 19) and that the dominant uncertainties on the evidence integral (eq. 21) and posterior scale as $\mathcal{O}\mathopen{}\mathclose{{\left(1/\sqrt{n_{\text{live}}}}}\right)$. The above considerations are represented graphically in fig. 2a.

There is thus a trade-off between run-time and uncertainty. Despite that, there isn’t a straight-forward method for choosing the number of live points if the required compression isn’t known ahead of time. The number should be chosen, furthermore, bearing in mind the alternative role that it plays as a resolution parameter for NS[48, 49], especially in multi-modal problems. In particular, $n_{\text{live}}$ should be large enough that at any time the constrained prior splits into disjoint modes, at least one live point lies inside the footprint of each mode. As a rough rule of thumb, if the constrained prior occupies a total volume $X$, only modes with a footprint greater than about $X/n_{\text{live}}$ may be reliably found[8, 49]. This defines a resolution down to which the posterior is reliably sampled. Modes with smaller footprints are typically not located and correctly sampled, and hence will also not contribute to the evidence estimate. Moreover, to sample reliably and efficiently from the constrained prior, it is usually advisable that $n_{\text{live}}$ exceeds the dimensionality of the parameter space.

[box:smc]Sequential Monte Carlo Connections between NS and a form of rare-event SMC sampler [50] can be exploited to develop an NS-SMC algorithm with desirable theoretical properties [21].

SMC samplers evolve $n$ samples (referred to as “particles”) through a series of distributions $\pi_{t}$ for $t=0,\ldots,T$ using reweighting, resampling and mutation steps. Of particular interest is a form of rare-event SMC [50] that involves sampling from the sequentially constrained prior distribution. This SMC sampler has the sequence of distributions

where $\mathbb{1}(\cdot)$ denotes an indicator function and $0=L^{\star}_{0}<\cdots<L^{\star}_{t}<\cdots<L^{\star}_{T+1}=\infty$. Each iteration $t$ of the sampler involves sampling $n$ particles with replacement from the current set of points satisfying $L(\Theta)\geq L^{\star}_{t}$ and then diversifying those particles through several iterations of an $\eta_{t}$-invariant MCMC kernel.

Noticing the similarities between the two algorithms, ref. [21] built upon these rare-event SMC samplers to create the NS-SMC algorithm for evidence and posterior estimation. Approximation of $Z$ and $P$ in NS-SMC is done by weighting the particles from iteration $t$ to target a shell of the posterior,

While original NS and NS-SMC both sample from the sequentially constrained prior distribution, there are several key differences. Crucially, NS-SMC uses weights based on importance sampling rather than numerical quadrature, it naturally handles the region beyond the largest threshold with no truncation error and it uses a different sampling mechanism. As a result, NS-SMC avoids the issues with bias in NS. Under mild conditions, including the situation where a finite number of MCMC steps are used, NS-SMC produces unbiased and consistent estimates of $Z$ and consistent estimates of $P$ as $n\rightarrow\infty$ [21]. Practically, NS can result in noticeably biased estimates of the evidence compared to NS-SMC for the same computational cost when the MCMC kernel is inefficient.

We have so far only considered NS with a fixed number of live points, and noted that the uncertainties in both posterior distributions and evidence estimates are reduced by increasing this number. However, the evidence depends on an accurate estimate of the total compression when we reach the posterior bulk. The posterior, on the other hand, depends only on an accurate estimate of the relative compression once inside the posterior bulk. The former uncertainty cancels in the posterior weights in eq. 15, as they are invariant under rescaling the estimates of the volume variable. This reflects the fact that whereas the evidence may depend strongly on the size of the prior, the posterior usually depends only weakly on its shape.

We are thus motivated to consider a dynamic number of live points to efficiently reduce uncertainties in parameter inference. As shown in fig. 2a, we may quickly compress to the posterior bulk using few live points. Upon reaching it, we may increase the number of live points, reducing uncertainty in the bulk of the posterior mass. We denote schemes that vary the number of live points as dynamic NS [51]. Open-source dynamic NS software packages include dynesty [52] and dyPolyChord [53]; see table 2. The gains are greatest in problems with substantial compression to the posterior bulk.

In dynamic NS schemes, the number of live points can be automatically adjusted to maximize a user-specified objective for a fixed computational budget. Usually, the run starts with an exploratory NS run using a constant number of live points. We spend the remaining computational budget by repeatedly increasing the number of live points in the most important regions of volume judged according to the objective (for example the shaded region in fig. 2a). Because dynamic NS re-winds a run and adds extra samples anywhere, running for longer reduces uncertainties and increases the effective sample size, unlike in ordinary NS.

In the original dynamic NS algorithm [51], a user specifies their objective by assigning a relative importance to reducing uncertainties in the posterior and the evidence. When focusing on posterior inferences, dynamic NS can achieve orders of magnitude reductions in computational cost for a fixed uncertainty. The approach can also improve evidence calculation accuracy for a fixed number of samples, and can improve posterior inferences and evidence calculations simultaneously. This objective was generalised by the reactive NS [54] variant of dynamic NS. This considers the computation as a graph, with the nodes being the live and dead points, and edges indicating the replacement of a point by another. Multiple agents can then add live points (edges) where needed, and optimize towards additional goals, such as the effective sample size, or the number of samples per cluster. Lastly, we note a significant variant of NS, diffusive NS [55], in which the number of live points at a given likelihood threshold can change. In this variant, random walks starting from the existing live points are permitted to step down as well as up in likelihood, refining the typical likelihood in a volume range $X$.

NS progresses by replacing live points by independent samples drawn from the constrained prior in eq. 9, that is the prior restricted to regions in which the likelihood exceeds a threshold. This is the major difficulty in efficiently and reliably implementing NS, especially in multi-modal problems. Provided one successfully samples from the constrained prior, however, NS works identically in uni-modal and multi-modal settings.

Whilst we could simply sample from the entire prior until we find a sample for which the likelihood exceeds the threshold, this rapidly becomes incredibly inefficient due to the exponential reduction in the volume contained within the constrained prior at every iteration. Fortunately, the current set of live points and the estimate of the volume enclosed by the contour may guide our search for new live points. There are two main classes of methods for sampling from the constrained prior [56]: region samplers and step samplers. They are illustrated in fig. 2b. Analogous to the choices of transition kernels in MCMC, the choices here lead to various flavors of NS with different performance characteristics and different behaviour as dimension grows. Although a priori they require just as much tuning as MCMC, the live points allow them to build proposal structures and apply clustering algorithms (analogously to ensemble samplers[57]), such that NS is naturally amenable to being robustly self-tuning. As they are guided by them, their reliability and efficiency usually improve when the number of live points is increased; see ref. [58, 54] for numerical investigations. For both region and step samplers, if no current live points lie inside a mode in a multi-modal problem, that region of the constrained prior almost certainly won’t be sampled, i.e., they miss that mode (see section 2.1).

Both samplers usually operate in the hypercube, a parameterization in which the prior is uniform over a unit hypercube (though this is not a requirement in NS). This slightly simplifies the problem to that of sampling uniformly from within a contour defined by the threshold. This is usually achieved by the inverse transformation method. In this case, users specify their priors by an inverse-cumulative density function rather than a density function. For example, suppose we desired a Gaussian prior for a parameter, $\theta\sim\mathcal{N}(\mu,\sigma^{2})$. We could transform a unit hypercube parameter, $u\sim\mathcal{U}(0,1)$, using the standard normal distribution’s inverse-cumulative density function, $\Phi^{-1}$,

This in principle allows all manner of priors, though see 6 Limitations and optimizations for further discussion.

To utilize computing resources, we could parallelize the computation of the likelihood function. If that is impossible or impractical, we may wish to design an efficient scheme for parallelizing an NS run itself. As discussed in ref. [51, 52], statistically independent NS runs may be combined into an equivalent NS run with $\sum n_{\text{live}}$ live points. To achieve this, pool the initial live points from the independent runs $a$, $b$, $\ldots$ together, giving $\sum n_{\text{live}}$ live points drawn from the prior. Remove the worst, supposing it lay at $L^{\star}$ and originated from run $a$. To draw a replacement from the prior subject to $L>L^{\star}$, simply take the replacement used in run $a$, since it is already a draw from the prior subject to $L>L^{\star}$. We may continue in this fashion, weaving together independent NS runs to build a new NS run with $\sum n_{\text{live}}$ live points. This allows parallelization of an NS run with $n_{\text{live}}$ live points into $n_{\text{CPU}}$ independent runs with about $n_{\text{live}}/n_{\text{CPU}}$ live points each. The independent runs themselves proceed linearly and may be ultimately combined, resulting in a speedup of about $n_{\text{CPU}}$. Even simpler, estimates of the evidence integrals may be combined by weighted averaging.

However, the reduction in the number of live points per run, $n_{\text{live}}/n_{\text{CPU}}$, impacts the exploration schemes discussed in section 2.3, especially for ellipsoidal rejection sampling where it could lead to inefficient or faulty bounding ellipsoids. It may therefore be desirable to utilize parallelization within individual NS runs. For example, by drawing $n_{\text{CPU}}$ candidate replacement points and evaluating their likelihoods in parallel. We could subsequently replace the worst $n_{\text{CPU}}$ live points at each iteration [75, 76], replace a single live point and consider the other evaluated points at subsequent iterations [65], or replace a single live point and discard as many as $n_{\text{CPU}}-1$ acceptable candidate points [48]. The latter is wasteful if more than one viable point are likely to be found among the $n_{\text{CPU}}$ candidates,

If points are considered at subsequent iterations,

and so a speedup of about $n_{\text{CPU}}$ if $n_{\text{CPU}}\ll n_{\text{live}}$. The expression originates from the fact that the threshold increases as the run progresses, meaning that points drawn from the constrained prior might not be valid at a subsequent iteration. Lastly, replacing $n_{\text{CPU}}$ points in parallel per iteration results in a speedup of about $n_{\text{CPU}}$ but increases the variance in the evidence integral by about $\sqrt{n_{\text{CPU}}}$. Increasing $n_{\text{live}}$ by a factor $\sqrt{n_{\text{CPU}}}$ to maintain the same uncertainty decreases the speedup to about $\sqrt{n_{\text{CPU}}}$.

We must decide when to stop an NS run. The fact that we can only perform a finite number of iterations introduces a truncation error in eq. 13 that we wish to be negligible. Skilling originally proposed to stop NS once we reached the posterior bulk at $X\simeq e^{-H}$ at iteration $n_{\text{iter}}\simeq n_{\text{live}}H$, or using an estimate of the remaining evidence. Popular NS software later adopted the latter. In MultiNest, this was based on the maximum likelihood found so far, $\max LX$, whereas PolyChord chose the mean likelihood, $\bar{}LX$. They stopped once $\Delta Z/Z\leq\textsf{tol}$, where tol is a user-specified parameter. Upon deciding to stop, the truncation error in the evidence was corrected by either adding an estimate of the remaining evidence or by killing the live points one by one without replacement and incrementing the evidence in the usual manner until no live points remain. The latter is in keeping with the NS approach. For the former the remaining evidence may be estimated by $\bar{}LX$. See ref. [43] for further discussion of the statistical properties of this estimate of the remainder. When NS is used to calculate the partition function of a material system (), physically motivated stopping conditions can be based on the expected minimum energy (negative log likelihood) or the sampled temperature, which is proportional to the derivative of the limiting energy with respect to NS iteration [77]

None of these approaches guarantees that summation hasn’t been terminated too early; there could well be a spike of enormous likelihood lurking inward. The computational budget, furthermore, cannot be easily anticipated ahead of time. However, runs that are terminated prematurely may still be used to illustrate what was learned so far about the posterior. Proposals to construct termination criteria for a fixed computational budget that result in unbiased estimates of the evidence have been suggested [22].

The rapid spread of Bayesian methods in cosmology [9] in the early 2000s was generated by the growth of data, specifically the new cosmic microwave background (CMB) temperature power spectrum measured by the Wilkinson Microwave Anisotropy Probe (WMAP) [115], and the type-Ia supernovae distance measurements [116, 117]. However, these new and powerful cosmological datasets had raised difficult questions regarding the cosmological model. The favoured model of the Universe included a mysterious accelerating force, the dark energy, which accounted for 70% of the energy density today. And the initial spectrum of density fluctuations in the early Universe, which were shown to be Gaussian and adiabatic by WMAP, indicated a period of accelerated expansion at very early times (about $10^{-23}$ seconds after the Big Bang), known as cosmic inflation [118]. Finally, there was still the question of the missing mass of the Universe, which generates the gravitational fields required for cosmic structure, known as cold dark matter (CDM). All three of these phenomena (cold dark matter, dark energy, and cosmic inflation) had proposed explanations from the field of high-energy theoretical physics, and these model predictions could be combined with the new wealth of cosmological data to be evaluated with respect to each other, model by model. Thus model selection in general, and the computation of the Bayesian evidence using NS specifically, became a tool of choice.

The simplest model of cosmic inflation is one driven by a single scalar field, a particle physics object similar (but not identical) to the Higgs boson. The behaviour of the scalar field driving inflation is determined by its potential $V(\phi)$ ($\phi$ being the value of the scalar field)[118, 119]. The formulation of $V(\phi)$ (as a function of some set of parameters $\Theta$) determines the duration of inflation along with the distribution of anisotropies found in the CMB. Ref. [120] performed a detailed Bayesian model comparison between 193 inflationary models using the region sampler MultiNest. They used the cosmological observations (including CMB data from WMAP mission) to discriminate between alternative models for $V(\phi)$ and found slight preference for so-called Small Field Inflation (SFI) models over Large Field Inflation (LFI) models.

Another area where NS has been used extensively in cosmology is the modelling of galaxy clusters [121, 122]. Clusters of galaxies are the most massive gravitationally bound objects in the Universe and therefore, can be used to trace the formation of large scale structure in the Universe. Galaxy clusters can be observed through several methods including X-ray observations, weak gravitational lensing and by exploiting the Sunyaev–Zel’dovich effect (SZ). Weak gravitational lensing involves the distortion of the images of the background galaxies by the presence of a large mass lying along the line-of-sight. Weak lensing allows to probe the total mass distribution, including the dark matter, of the galaxy cluster. Ref. [123] presented a Bayesian approach using NS for the detection of galaxy clusters in SZ data, including estimation of parameters associated with the physical model assumed for the cluster and quantification of the detection using Bayesian model selection, making use of the statistics of the CMB anisotropies in the likelihood function. Ref. [124] presented a joint Bayesian analysis of weak lensing and SZ observations of several galaxy clusters, allowing for the estimation of gas fraction of individual clusters. In modern weak lensing surveys [125] such as DES [126] and KiDS [127], NS forms a critical part of their parameter estimation, model comparison and tension quantification pipelines [128] (fig. 3).

The expansion history of the Universe can also be measured through the use of distance measurements, such as observations of type Ia supernovae (SNIa), which can be used as ‘standard candles’. The important results from the Supernova Cosmology Project and the High-Z Supernova Search Team [116, 117] presented evidence for the accelerated expansion of the Universe, requiring the existence of a mysterious dark energy acting against gravity. Estimation of cosmological parameters from the observations of SNIa light curves has historically been done using a $\chi^{2}$ approach (see e.g. [129]) which lacks a rigorous approach for the determination of systematic uncertainties. Ref. [130] introduced a Bayesian hierarchical model for the determination of cosmological parameters. By using NS, they showed that their principled Bayesian approach delivered tighter statistical constraints on the cosmological parameters over 90% of the time, reduced statistical bias by a factor $\sim$ 2-3 times and that it had better coverage properties than the standard $\chi^{2}$ approach.

The measurement of the CMB anisotropies by Planck [132] was an increase in the statistical power over the previous experiment (WMAP), but required more sophisticated modeling of galactic foregrounds and instrumental calibration, introducing $\mathcal{O}\mathopen{}\mathclose{{\left(20}}\right)$ “nuisance parameters”. Metropolis-Hastings techniques for parameter estimation were able to accommodate these by exploiting the fact that these parameters were “fast” in comparison to the cosmological ones; namely, by caching results from previous calculations, these parameters can be changed with negligible computational cost providing the cosmological parameters remain fixed.

Region-based samplers such as MultiNest are not able to exploit this hierarchy of parameter speeds, and could not navigate the now high-dimensional cosmological + nuisance parameter space reliably. To address this, cosmologists turned to step sampling based strategies as instantiated in the PolyChord algorithm [133, 65] which uses slice sampling.¹¹ 1 Axial slice sampling had been applied in the then near past to systems biology [134] although unsurprisingly the cosmological authors were not aware of this work. PolyChord improved upon the existing slice-sampling method by implementing covariance-based non-axial steps and mode clustering in addition to the ability to exploit a hierarchy of parameter speeds. These were successfully applied throughout the 2015 Planck inflation paper [135] to non-parametric reconstructions [136] (fig. 3) and general inflationary model comparison, in particular to the challenging example of axion monodromy models. The non-parametric reconstruction approach, which only became possible with the ability to reliably (and fully) explore the parameter space of a large number of parameters, was an important model-independent demonstration of the simple power-law behavior of the primordial power spectrum of density perturbations.

Although step sampling was originally introduced in cosmology to exploit a fast-slow hierarchy of parameters, the better dimensional scaling opened up a new range of cosmological analysis possibilities that were inaccessible to region samplers. It has been applied to constraining kinetically dominated inflation models [137], model comparison for the quantum mechanical initial conditions for the Universe [138], additional reconstructions of the dark energy equation of state [139, 140], astronomical sparse reconstruction [141], and played a critical role in the GAMBIT combined cosmology and particle analyses [142, 143]. Step sampling takes a leading role in the REACH 21cm global cosmology analysis [144], and at the other end of the astrophysical scale it has also been applied to high-dimensional exoplanet analyses [111, 113].

Particle physics is a field related to cosmology that has also seen various applications of NS. Around 2010, when NS tools such as MultiNest were reaching maturity, the particle physics community was focused on the first results from the Large Hadron Collider (LHC). A particularly favoured theoretical framework was supersymmetry (SUSY; see e.g., ref. [145]). SUSY introduces an array of new particles and unknown parameters that can be fitted to collider data as well as observations from other experiments in a so-called global fit. The physics goals of a global fit are to understand what the model predicts in future experiments and how best to discover the new particles.

These global fits of $\mathcal{O}\mathopen{}\mathclose{{\left(10}}\right)$ free parameters in SUSY models presented a problem for which NS tools were naturally well-suited, as they were typically multi-modal. A package named SuperBayeS [147, 148, 146] utilized the MultiNest implementation of NS and was used to make a number of early LHC predictions and fits[149, 150, 151, 152, 153, 154]. As the LHC results have poured in throughout the 2010s, the theoretical landscape shifted and a wider set of models have been considered using NS (see e.g., ref. [155, 156, 157, 158]). More recently, the GAMBIT collaboration has driven many such global fits making use of MultiNest and PolyChord to sample parameter spaces and compute Bayesian evidence s, and benchmarked NS against a number of other MC based and gradient-free methods [159, 160]. Related problems such as tuning phenomenological parameters in event generators have seen preliminary work, and could be rich avenues of NS application in the field. Lastly, we note that NS was recently applied to the sampling space, rather than the parameter space, of a statistical model. This enables efficient computation of small $p$-values that are used in the discovery of new particles at the LHC[161].

Gravitational-wave (GW) astronomy is a field that has blossomed since the first observation in 2015 of two colliding black holes [162] by the LIGO [163] and Virgo [164] interferometers. The signals are produced by non-axisymmetric changes in the gravitational field, typically sourced by the rapid motion of neutron stars and black holes in binary systems [165]. The binary orbits decay through GW emission, increasing the orbital frequency and the rate of energy loss until the objects merge. Such events, known as compact binary coalescences (CBCs), produce signals at frequencies of 10–1000 Hz which are recorded in the detectors as a time-series.

Early development of Bayesian methods for GW analysis was contemporaneous with the first publication of NS, and it was found that NS provided an efficient means to robustly sample the posterior distributions of GW signals [166]. This is important as the posteriors provide a rich new astrophysical view: from measuring the masses of neutron stars to the expansion rate of the Universe itself. The small signal-to-noise ratios of observed signals and degeneracies in the model parameter space produce posterior distributions that are often multi-modal and highly correlated. To date, only NS and MCMC approaches have been able to robustly sample the posteriors of all observed events. While each has its own advantages, agreement between the NS and MCMC approaches has been critical in building confidence in results. However, the NS approach is better-suited to providing robust evidence estimates for model comparisons [167]; the efficiency of the NS approach does not depend so much on the use of problem-specific proposals; and massively parallel approaches using the dynesty code have made analyses with more advanced signal models [168] computationally tractable. The success of NS in analysing merging binary signals has inspired many other efforts such as the analysis of continuous signals from individual rapidly-rotating non-axisymmetric neutron stars [169, 170]; detection of unmodelled sources [171]; model selection between different physical mechanisms of core-collapse supernova encoded in the GW signal [172]; and detection of a stochastic superposition of weak merger sources [173].

CBCs consisting of two black holes have been the dominant sources of GW signals seen by the LIGO and Virgo detectors [165]. The signals are modelled through a combination of post-Newtonian approximations to General Relativity and relativistic numerical modelling [174, 175, 176, 177]. The signal model, as observed in a detector, is parameterized by 8 parameters that are intrinsic to the binary system (the individual masses, $m_{1}$ and $m_{2}$, and their 3D spin vectors) and 7 parameters related to the relative orientation and position of the system with respect to the detector (including the source’s luminosity distance $D_{\rm L}$ and location on the sky). For systems that include at least one neutron star, the signal model’s phase evolution requires additional parameters related to the neutron star equation-of-state [178].

Signals from black hole mergers last only a few seconds within the sensitive frequency regime of the detectors. On these timescales the detector data can be approximated as having noise drawn from a stationary Gaussian process described by a known power spectral density, which can either be estimated from data surrounding the signal [180], or inferred directly using a parameterized model of its shape [181]. For inference, a Gaussian likelihood function of the form given by Whittle [182] can be used [180]. For multiple detectors, which have independent noise, the likelihoods can be coherently combined using the product rule of probability. The prior distributions used are discussed in ref. [180], but are generally set to be uninformative (e.g., uniform over a sphere for sky coordinates), or constrained to be within a physically reasonable range. In addition to the parameters related to the source, the likelihood can also contain $\sim$ 60 unknown parameters that relate the frequency-dependent uncertainties in the phase and amplitude calibration of the detectors [183].

For an observed signal, using the likelihood and priors discussed above, the joint posterior of all these 15 (or more) parameters can be extracted after application of NS. The posterior distribution of the typical events observed so far are not uni-modal or Gaussian. Due to the relatively small signal-to-noise ratio, they exhibit significant degeneracies and correlations. For example, the posterior of the source sky location is largely determined by differences in arrival time of the signal at different detectors. This produces ring-shaped degeneracies as demonstrated in fig. 4b.

The many degeneracies encountered for typical CBC inference have led to significant work to identify optimal parameterization and, where possible, utilized a marginalized likelihood (see ref. [184] for a review). An optimal parameterization involves identifying a mapping between the physical model parameters and combinations of these which reduce the complexity of the target density. Usually, a good reparameterization involves identifying the combinations of physical parameters which are ‘best measured’. As an example, the physical mass of the two component stellar objects are $m_{1}$ and $m_{2}$. However, there is a strong banana-like correlation between the two masses (see, e.g. fig. 4b) and an exact degeneracy under exchange of $m_{1}$ and $m_{2}$. To enable efficient sampling, Veitch et al [180] proposed sampling in the chirp-mass, $\mathcal{M}$, and mass ratio $q$: two algebraic combinations of the component masses. In fig. 4b, we show that the posterior, as viewed in $m_{1}-m_{2}$ space, follows contours of the chirp mass. Physically, this is because the chirp-mass is the most well measured mass parameter (i.e. has the smallest posterior width) followed by the mass ratio.

The chirp-mass and mass ratio parameterization, along with numerous others have greatly improved sampling efficiencies. We note that the choice of sampling parameters, which we choose for computational efficiency, is separate from the choice of prior distributions. If required, a Jacobian transformation [185] can be made to enable sampling in the optimal sampling parameters while setting priors on the physical parameters.

As hinted at in above, NS can be used to study the thermodynamic properties of molecules and materials, which ultimately derive from the partition function

where $q\in\mathbb{R}^{3N}$ specifies the spatial coordinates of the $N$ particles in the system, $p\in\mathbb{R}^{3N}$ is the corresponding momentum vector, $H(q,p)$ is the Hamiltonian, and $\beta=1/k_{B}T$ is the inverse temperature.

The classical Hamiltonian can be separated into the configuration dependent potential energy, $U(q)$, and momentum dependent kinetic energy, $K(p)=\sum p^{2}/2m$ (where $m$ is the mass of each particle), giving

(For the fully quantum treatment using NS, see ref. [186].) The first factor is a Gaussian that may be computed analytically,

and we tackle the second, configuration dependent factor by NS. From the partition function all thermodynamic quantities of relevance can be calculated, for example, the average energy,

and the heat capacity

and these may be interpreted as posterior expectations to be computed after the NS run. In this application,

• •

  the factor $e^{-\beta U(q)}$ plays the role of the likelihood function, parameterized by the inverse temperature, $\beta$;

• •

  the prior, according to the ergodic hypothesis, is uniform;

• •

  the dimensionality is commonly on the order of $10^{2}$ to $10^{3}$ or more.

The striking difference from standard Bayesian inference is the inverse temperature parameter $\beta$ in the likelihood. Rather than having a single inference problem for some fixed value of $\beta$, almost always we want to know the behaviour of observables as a function of temperature, so effectively we have a continuous family of inference problems. And for this, NS (and other density-of-states methods, see ) have a remarkable feature: the likelihood is a monotonic function of $\beta$ and hence the entire NS algorithm is invariant to changes in $\beta$. In practice this means that a single NS run can be used to calculate observables at all temperatures. For large $\beta$ values (corresponding to low temperatures), the partition function is dominated by the lowest energy minima (highest likelihood modes), and for molecules and materials these configurations correspond to the globally stable structures. On the other hand, when $\beta$ is small (corresponding to high temperatures), the partition function is dominated by the large volume associated with high energy states — or in the parlance of materials science, “entropic effects”[187, 188].

Some of the most interesting phenomena in molecular and materials science (and more generally in statistical mechanics) are associated with the above change of regime. Collectively known as phase transitions, they are characterized by a dramatic change of where the bulk of the posterior mass lies, as the temperature (or other system parameters such as pressure) is varied — and this is what makes it very challenging to study them using numerical sampling schemes. Experimentally, phase transitions are often observed indirectly as changes in the expectation value of observables (e.g. discontinuously, in the case of first order phase transitions such as melting and evaporation), or more directly as sharp peaks in response functions, such as the heat capacity or magnetic susceptibility. The promise of NS is to enable the calculation of such response functions in general with high reliability and the minimum of fuss.

There are some aspects of the materials application of NS that turn out to be favourable in comparison with the general inference problem. The first is a rather easy stopping criterion for the NS iterations. For any given model of the potential energy, it is generally not hard to come up with a good global lower bound on the energy, which translates into an upper bound in the likelihood. As the likelihood values sampled by NS appear to converge, if this is close to the known bound, the iterations can be stopped without the risk of missing the highest likelihood mode.

Secondly, convergence of NS with the number of live points and other sampling parameters is desired and observed in terms of convergence of the heat capacity peak locations, and this typically occurs far earlier than the decorrelation of the sampler chains that are used to explore the constrained prior.

Although NS is broadly applicable, there are several potential limitations. The first limitation relates to the prior: we sample from the constrained prior and thus require a proper prior. For many NS implementations, this proper prior must be transformed from the unit hypercube. Normalizing flows were recently proposed [213] for cases where this is inconvenient or impossible analytically, including using the posterior from an NS run as a prior.

There are, furthermore, limitations related to the likelihood. Integration using NS requires a non-negative integrand (though the compression itself makes no such restriction). Whilst this condition is always fulfilled in statistical applications, it could be violated when NS is used as a general purpose integrator. We wrote eq. 6 assuming that $L\geq 0$ whereas eq. 8 makes no restriction. NS, however, compresses upwards in likelihood, such that positive and negative likelihood regions of the integral would be treated differently, with the latter explored at inadequate resolution. NS, furthermore, requires a tractable likelihood. For cases in which the likelihood is intractable, ref. [99] proposed a likelihood-free NS in the context of systems biology, assuming that an unbiased estimator of the likelihood was available.

Lastly, plateaus in the likelihood function spoil the estimates of the compression [2, 214, 215, 216]. That is, sets $A$ with non-negligible prior mass, $\mu(A)>0$, and constant likelihood, which means $L(\Theta)\equiv c$ for all $\Theta\in A$. We may modify the likelihood function to remove plateaus by adding a negligible unique iid tie-breaking random draw to the likelihood or promoting that draw to a parameter and increasing the dimension of the problem. We may alternatively modify the algorithm itself such that it sums plateaus correctly but reduces to the ordinary NS algorithm in their absence. In the presence of plateaus, the minimum likelihood among the live points, $\min L$, may be shared by several live points.

For example, we could modify NS by removing all $q$ points in the plateau and then replacing them all [217]. In this case, the compression factors follow $\text{Beta}(n_{\text{live}}+1-q,q)$. This may be applied retrospectively to any NS runs, though suffers from increased uncertainty as the number of live points is temporarily reduced to as few as $n_{\text{live}}-q$. Alternatively, we could first increase the live points by sampling subject to $L>L^{\star}$ until $n_{\text{live}}-1$ lie at $L>\min L$. Note that $\min L$ may decrease as we add live points at $\min L>L>L^{\star}$. We then remove $q$ points in the outermost contour. The compression factor follows $\text{Beta}(n_{\text{live}},q)$. This cannot be applied retrospectively but the uncertainties are reduced (at the cost of greater run-time) as the number of live points is temporarily increased beyond $n_{\text{live}}$. This may in fact be seen as a dynamic version of the first modification and we could increase the number of live points according to a different criterion. Lastly, we could split the integral into an integral over plateaus and an integral without plateaus, and perform only the latter with NS [216].

In addition, there are computational limitations. In ordinary NS running for longer doesn’t reduce uncertainties or increase the number of posterior samples. If one is only interested in using NS for optimization, this isn’t a drawback. In any case, as discussed in 2.2 Dynamic nested sampling, dynamic NS overcomes this problem by rewinding and resuming an NS run. Besides that, NS requires $\mathcal{O}\mathopen{}\mathclose{{\left(D^{2}}}\right)$ iterations (eq. 19), the memory required to store the co-ordinates of every dead point scales as $\mathcal{O}\mathopen{}\mathclose{{\left(D^{3}}}\right)$ (with potentially nasty scaling factors for clustering in particular exploration strategies). This means that NS implementations that store every dead point become memory bound at $\mathcal{O}\mathopen{}\mathclose{{\left(500}}\right)$ dimensions.

We stress, furthermore, that users should be aware of ways in which NS may fail:

1. 1.

   NS may fail to successfully draw independent samples from the constrained prior. Consequently NS results including error estimates may be faulty and anticipated properties of NS, such as convergence, won’t hold. In 3 Results, we discussed cross-checks on NS results, such as ref. [89], that may identify this issue.

2. 2.

   Owing to item 1, NS may miss modes. multi-modal problems pose challenges in Bayesian computation, particularly in MCMC, where the chain must make a sequence of unlikely steps between modes. Unfortunately, it is in general impossible to know whether all modes have been found. Broadly speaking, as NS does not depend on slow transitions between modes, it is well-suited to multi-modal problems. Once a mode is established, furthermore, it won’t be abandoned until the likelihood threshold exceeds that of the points that were in the mode.

3. 3.

   NS may sample inefficiently from the constrained prior. As discussed in 2.3 Exploration strategies, this may occur in high-dimensional problems with rejection sampling strategies.

Our check-list in includes checks of NS failures.

There are several ways NS runs may be optimized to make best use of computing resources (see also 2.4 Parallelization). First, one may utilize fast and slow parameters by breaking the likelihood function into fast and slow operations that involve subsets of the $D$ parameters. The parameters associated with fast and slow operations are referred to as fast and slow parameters, respectively; for example, if the likelihood function may be written as,

then $y,z\ldots$ are fast parameters and $x$ is a slow parameter. When selecting a new point, to minimize runtime we should where possible change the fast parameters and keep the slow parameters constant, allowing caching of the slow operation [218]. This is natural in the slice sampling exploration strategy discussed in 2.3 Exploration strategies, as we may pick slices along fast directions [65]. If $D_{s}$ parameters are slow, we require $D_{s}$ rather than $D$ slow operations per iteration. Similarly, it may be exploited in modified Metropolis algorithms [24, 219] in which the proposals change blocks of parameters at a time.

In addition, one may repartition the posterior. NS must compress exponentially through the entire prior volume, which could be slow for a diffuse prior with substantial compression to the posterior. This in practice limits NS to $\mathcal{O}\mathopen{}\mathclose{{\left(100s-1000s}}\right)$ of dimensions. This prohibits fitting complex hierarchical Bayesian models and deep neural networks, or reaching HMC like dimensionalities of $\mathcal{O}\mathopen{}\mathclose{{\left(10^{6}}}\right)$. In some cases, it may be possible to use a narrower prior and correct the evidence estimates from NS post-hoc.

This issue may occur in problems in which the prior is unrepresentative [220], i.e., where the observed data lies in the tail of the prior predictive distribution. This may be mitigated by Bayesian automatic prior repartitioning [221] in which one redefines the prior and likelihood while leaving their product unchanged. This allows one to keep the evidence integral in eq. 1 unchanged, but reduce the prior to posterior compression in eq. 16. For example, consider a Gaussian likelihood, $\mathcal{N}(x,\sigma^{2})$, for a parameter $x\sim\mathcal{U}(-L/2,L/2)$ for $\sigma\lll L$. The compression would be approximately $H\approx\log(\sigma/L)$. Repartitioning the prior and likelihood by swapping them, we obtain $H\approx 0$.

The possibility of improving the robustness and efficiency of NS by exploiting the intrinsic degeneracy between the ‘effective’ likelihood and prior in the formulation of Bayesian inference problems was also discussed in ref. [123], and posterior repartitioning can further be viewed as the vanilla case (when the importance weight function equals to $1$) of nested importance sampling proposed in ref. [18].

Lastly, rather than performing new NS runs, one may reuse runs for similar likelihoods and priors. If the prior or likelihood are modified, it may be possible to reweight the posterior weights and evidence integral [40],

where $P_{i}=w_{i}L_{i}$. This re-weighting may be interpreted as a pseudo-importance sampling in which the original estimated posterior plays the role of the kernel. This is particularly useful for investigating prior sensitivity in the context of Bayesian inference. The effective sample size of draws from the new posterior may be used to judged the reliability of this procedure; if the new and original posterior distributions differ substantially, a fresh NS run needs to be performed.

There may, furthermore, be cases in which we wish to investigate several similar likelihood functions at once; ref. [59] presents a collaborative version of NS that operates on more than one likelihood function at once and in which parts of the likelihood evaluation may be recycled.