IMPORTANCE NESTED SAMPLING AND THE MULTINEST ALGORITHM

a (pseudo-)importance sample the totality of points collected by M ULTI N EST , including those previously discarded under the constrained likelihood sampling of the NS algorithm. We apply this technique to several challenging test problems and compare the accuracy of Bayesian evidences obtained with INS against those from vanilla NS.


INTRODUCTION
The last two decades in astrophysics and cosmology have seen the arrival of vast amounts of high quality data. To facilitate inference regarding the physical processes under investigation, Bayesian methods have become increasingly important and widely used (see e.g. Trotta 2008 for a review). In such applications, the process of Bayesian inference may be sensibly divided into two distinct categories: parameter estimation and model selection. Parameter estimation is typically achieved via MCMC sampling methods based on the Metropolis-Hastings algorithm and its variants, such as slice and Gibbs sampling (see e.g. Mackay 2003). Unfortunately, these methods can be highly inefficient in exploring multi-modal or degenerate distributions. Moreover, in order to perform Bayesian model selection (Clyde et al. 2007), estimation of the Bayesian 'evidence', or marginal likelihood, is needed, requiring a multi-dimensional integration over the prior density. Consequently, the computational expense involved in Bayesian model selection is typically an order of magnitude greater than that for parameter estimation, which has undoubtedly hindered its use in cosmology and astroparticle physics to-date. Nested sampling (NS; Skilling 2004Skilling , 2006Sivia and Skilling 2006) is a contemporary Monte Carlo (MC) method targeted at the efficient calculation of the evidence, yet which allows posterior inference as a by-product, providing a means to carry out simultaneous parameter estimation and model selection (and, where appropriate, model averaging).  and  have built on the NS framework by introducing the now-popular MULTINEST algorithm, which is especially efficient in sampling from posteriors that may contain multiple modes and/or degeneracies. This technique has already greatly reduced the computational cost of Bayesian parameter estimation and model selection and has successfully been applied to numerous inference problems in astrophysics, cosmology and astroparticle physics (see e.g. Feroz et al. , 2010Feroz et al. , 2011aBridges et al. 2009;Graff et al. 2012;White and Feroz 2010;Kipping et al. 2012;Karpenka et al. 2012Karpenka et al. , 2013Strege et al. 2013;Teachey and Kipping 2018).
In this paper, we discuss importance nested sampling (INS), an alternative summation of the draws from MULTINEST's exploration of the model parameter space with the potential to increase its efficiency in evidence computation by up to a order-ofmagnitude. Version (v3.0) of MULTINEST, which implements INS in addition to the vanilla NS scheme of previous versions, is available at https://github.com/farhanferoz/MultiNest.
The outline of this paper is as follows. We give a brief introduction to Bayesian inference in Sec. 2 and describe nested sampling along with the MULTINEST algorithm in Sec. 3. The INS technique is discussed in Sec. 4 and is applied to several test problems in Sec. 5. We summarize our findings in Sec. 6. Finally, in Appendix A we discuss the relationship between INS and other contemporary MC schemes, in Appendix B we give a detailed account of the convergence properties of INS within the MULTINEST algorithm, and in Appendix C we present a brief measure-theoretic commentary on vanilla NS. where Pr(Θ|D, H) ≡ P (Θ|D) is the posterior probability density of the model parameters, Pr(D|Θ, H) ≡ L(Θ) the likelihood of the data, and Pr(Θ|H) ≡ π(Θ) the parameter prior. The final term, Pr(D|H) ≡ Z (the Bayesian evidence), represents the factor required to normalize the posterior over the domain of Θ given by: Being independent of the parameters, however, this factor can be ignored in parameter inference problems which can be approximated by taking samples from the unnormalized posterior only, using standard MCMC methods (for instance).
Model selection between two competing models, H 0 and H 1 , can be achieved by comparing their respective posterior probabilities given the observed dataset as follows: Here Pr(H 1 )/ Pr(H 0 ) is the prior probability ratio for the two models, which can often be set to unity in situations where there is no strong a priori reason for preferring one model over the other, but occasionally requires further consideration (as in the Prosecutor's Fallacy; see also  for key astrophysical examples). It can be seen from Eq. (3) that the Bayesian evidence thus plays a central role in Bayesian model selection.
As the average of the likelihood over the prior, the evidence is generally larger for a model if more of its parameter space is likely and smaller for a model with large areas in its parameter space having low likelihood values, even if the likelihood function is sharply peaked. Thus, the evidence may be seen both as penalizing 'fine tuning' of a model against the observed data and as an automatic implementation of Occam's Razor.

NESTED SAMPLING AND THE MULTINEST ALGORITHM
Nested sampling estimates the Bayesian evidence by transforming the multi-dimensional evidence integral over the prior density into a one-dimensional integral over an inverse survival function (with respect to prior mass) for the likelihood itself. This is accomplished by considering the survival function, X(λ), for L(Θ), dubbed "the prior volume" here; namely, where the integral extends over the region(s) of parameter space contained within the iso-likelihood contour, L(Θ) = λ. Recalling that the expectation value of a non-negative random variable may be recovered by integration over its survival function (a result evident from integration by parts) we have (unconditionally): When L(X), the inverse of X(λ), exists (i.e., when L(Θ) is a continuous function with connected support; Chopin and Robert 2010) the evidence integral may thus be further rearranged as: Indeed, if L(X) were known exactly (and Riemann integrable 1 ), by evaluating the likelihoods, L i = L(X i ), for a deterministic sequence of X values, 0 < X N < · · · < X 2 < X 1 < X 0 = 1, as shown schematically in Fig. 1, the evidence could in principle be approximated numerically using only standard quadrature methods as follows: where the weights, w i , for the simple trapezium rule are given by w i = 1 2 (X i−1 −X i+1 ). With L(X) typically unknown, however, we must turn to MC methods for the probabilistic association of prior volumes, X i , with likelihood contours, L i = L(X i ), in our computational evidence estimation.

Evidence estimation
Under the default nested sampling algorithm the summation in Eq. (8) is performed as follows. First N live 'live' points are drawn from the prior, π(Θ), and the initial prior volume, X 0 , is set to unity. At each subsequent iteration, i, the point with lowest likelihood value, L i , is removed from the live point set and replaced by another point drawn from the prior under the constraint that its likelihood is higher than L i . The prior volume contained within this region at the i th iteration, is thus a random variable distributed as X i = t i X i−1 , where t i follows the distribution for the largest of N live samples drawn uniformly from the interval [0, 1] (i.e., Pr(t) = N live t N live −1 ). This sampling process is repeated until (effectively) the entire prior volume has been traversed; the live particles moving through nested shells of constrained likelihood as the prior volume is steadily reduced. The mean and standard deviation of log t, which governs the geometrical exploration of the prior volume, are: Since each draw of log t i is independent here, after i iterations the prior volume will shrink down as log X i ≈ −(i ± √ i)/N live . Thus, one may take X i ≈ exp(−i/N live ).

Stopping criterion
The NS algorithm should terminate when the expected evidence contribution from the current set of live points is less than a user-defined tolerance. This expected remaining contribution can be estimated (cautiously) as ∆Z i = L max X i , where L max is the maximum likelihood value amongst the current set of live points (with X i the expected value of remaining prior volume, as before).

Posterior inferences
Although the NS algorithm is designed specifically to estimate the Bayesian evidence, inferences from the posterior distribution can be easily obtained using the final live points and the full sequence of discarded points from the NS process, i.e., the points with the lowest likelihood value at each iteration of the algorithm. Each such point is simply assigned the importance weight, from which sample-based estimates for the key posterior parameter summaries (e.g. means, standard deviations, covariances and so on) may be computed 2 . (As a self-normalizing importance sampling estimator the asymptotic variance of these moments is of course dependent upon both the similarity between the NS path and the target and the accuracy ofẐ itself; cf. Hesterberg 1995.) Readers unfamiliar with importance sampling (IS) ideas may refer to Liu (2008) for an insightful overview of this topic and its application to diverse branches of modern science (including statistical physics, cell biology, target tracking, and genetic analysis).
3.4. Practical implementations of nested sampling The main challenge in implementing the computational NS algorithm is to draw unbiased samples efficiently from the likelihood-constrained prior. John Skilling, originally proposed to use the Markov Chain Monte Carlo (MCMC) method for this purpose (Skilling 2004(Skilling , 2006Sivia and Skilling 2006). One such implementation (Veitch and Vecchio 2010), with specific proposal distributions for the MCMC step, has been used successively in gravitational wave searches.
In astrophysics in particular, rejection sampling schemes have been successfully employed to draw samples from the likelihoodconstrained prior. It was first proposed in the COSMONEST package (Mukherjee et al. 2006) through the use of ellipsoidal rejection sampling scheme and was shown to work very well for uni-modal posterior distributions. This method was improved upon in COSMOCLUST package (Shaw et al. 2007) through the use of a clustering scheme to deal with multi-modal distributions. MULTINEST was then proposed with several innovations to make ellipsoidal rejection sampling more robust in dealing with multi-modal distributions. Other methods employing ellipoidal rejection sampling scheme within Nested Sampling framework include the DIAMONDS (Corsaro and De Ridder 2014) and DYNESTY (Speagle 2019) packages. One particular problem with rejection sampling schemes is the exponential reduction in sampling efficiency with increasing dimensionality of the problem. In order to address this issue, a slice sampling method has been employed to draw unbiased samples efficiently from the likelihood-constrained prior in the POLYCHORD (Handley et al. 2015a,b) package.
Another algorithm to increase the efficiency of Nested Sampling through the variable number of live points is the "Dynamic Nested Sampling" method (Higson et al. 2018) which has been used in the DYNESTY (Speagle 2019) package.

MULTINEST algorithm
The MULTINEST algorithm ) addresses this problem of drawing unbiased samples from the likelihood-constrained prior, through an ellipsoidal rejection sampling scheme. At each iteration, i, the full set of N live live points is enclosed within a set of (possibly overlapping) ellipsoids and the desired replacement point sought from within their union. The ellipsoidal decomposition of the live point set is performed through an expectation-minimisation algorithm such that the sum of volumes of the ellipsoids is minimised with the additional constraint that the total volume enclosed by the ellipsoids is at least X i /f . Again X i ≈ exp(−i/N live ) is the expected prior volume, while 0 < f ≤ 1 is a user defined value for the target efficiency (the ratio of points accepted to points sampled). Thus, f is analogous to the (inverse of the) "enlargement factor" introduced by Mukherjee et al. (2006) into their pioneering ellipsoid-based NS code; the larger the target f the faster the algorithm runs, but the greater the chance of some ellipsoids failing to cover the full L > L i volume (biasing the vanilla NS estimates, though not necessarily the INS estimates, as we discuss later).
The MULTINEST ellipsoidal decomposition algorithm thus allows substantial flexibility in the geometry of its posterior exploration; with bent and/or irregularly-shaped posterior modes typically broken into a relatively large number of small 'overlapping' ellipsoids and smooth, near-Gaussian posterior modes kept whole (or broken into relatively few ellipsoids), as shown in Fig. 2. It thereby automatically accommodates elongated, curving degeneracies while maintaining high efficiency for simpler problems. MULTINEST also specifically enables the identification of distinct modes by isolating non-overlapping subsets of the ellipsoidal decomposition; so identified, these distinct modes can then be evolved independently.
Once the ellipsoidal bounds have been created at a given iteration of the MULTINEST algorithm a new point is drawn uniformly from the union of these ellipsoids as follows. If there are L ellipsoids at iteration i, a particular ellipsoid is chosen with probability p l given as: where V tot = L l=1 V l , from which a single point is then drawn uniformly and checked against the constraint L > L i . If satisfied the point is accepted with probability 1/q, where q is the number of ellipsoids the new point lies in (in order to take into account the possibility of non-empty intersections), otherwise it is rejected (but saved for INS summation) and the process is repeated with a new random choice of ellipsoid.
In higher dimensions, most of the volume of an ellipsoid lies in its outer shells and therefore any overshoot of the ellipsoidal decomposition relative to the true iso-likelihood surface can result in a marked drop in sampling efficiency. In order to maintain the sampling efficiency for such high dimensional problems, MULTINEST can also operate in a 'constant efficiency mode'. In this mode, the total volume enclosed by the ellipsoids is no longer linked with the expected prior volume X i by requiring the total ellipsoidal volume to be at least X i /f , instead the total volume enclosed by the union of ellipsoids is adjusted such that the sampling efficiency is as close to the user defined target efficiency f as possible while keeping every live point enclosed in at least one ellipsoid. Despite the increased chance of the fitted ellipsoids encroaching within the constrained-likelihood volume (i.e., missing regions of parameter space for which L > L i ), past experience has shown (e.g. ) this constant efficiency mode may nevertheless produce reasonably accurate posterior distributions for parameter estimation purposes. The vanilla NS evidence values, however, cannot be relied upon in this mode, with the expectation being a systematic over-estimation of the model evidence. Interestingly, the same is not strictly true of the INS evidence estimates, which use the NS technique only for posterior exploration (not evidence summation); though we will note later some important caveats for its error estimation.
In the rest of this paper, we refer to the mode in which MULTINEST links the volume of the ellipsoidal decomposition with the expected prior volume as its 'default' mode, and we specifically highlight instances where 'constant efficiency mode' has been trialled.

IMPORTANCE NESTED SAMPLING
Though highly efficient in its approximation to the iso-likelihood contour bounding the live particle set at each iteration, the ellipsoidal rejection sampling scheme used by MULTINEST ultimately discards a significant pool of sampled points failing to satisfy the NS constraint, L > L i , for which the likelihood has nevertheless been evaluated at some computational cost. In order to redress this final inefficiency the technique of importance nested sampling (INS) has recently been proposed by Cameron and Pettitt (2013) as an alternative summation for the Bayesian evidence in this context. In particular, INS uses all the points drawn by MULTINEST, or any other ellipsoidal rejection sampling algorithm, at each iteration regardless of whether they satisfy the constraint L > L i or not. The relationship of INS to existing MC schemes is summarised in Appendix A.
4.1. Pseudo-importance sampling density One begins by defining the following pseudo-importance sampling density: where N iter is the total number of iterations (ellipsoidal decompositions) performed by MULTINEST, n i the number of points collected at the i th iteration (with total, N tot = Niter i=1 n i ), V tot,i the total volume enclosed in the union of ellipsoids at the i th iteration, and E i (Θ) an indicator function returning 1 when Θ lies in the i th ellipsoidal decomposition and 0 otherwise. We call g(Θ) here a pseudo-importance sampling density since it is of course defined only a posteriori to our sampling from it, with the consequence that all Θ ∼ E j>i (Θ) are to some (ideally negligible) extent dependent on all previous Θ ∼ E j≤i (Θ) (some important implications of which we discuss in Sec. 4.3 below). The heritage of this technique lies with the reverse logistic regression strategy of Geyer (1994) and the "biased sampling" framework of Vardi (1985). Another term that has been used in place of pseudo-importance sampling is "recycling of past draws" (e.g. Cornuet et al. 2012).
If at each iteration, the ellipsoidal decomposition would consist of only one ellipsoid then V tot,i is simply the geometric volume of the ellipsoid at iteration i. MULTINEST, however, may enclose its live points in a set of possibly overlapping ellipsoids. An analytical expression for calculating the volume in the overlapped region of ellipsoids is not available and therefore we estimate the volume occupied by the union of ellipsoids through the following MC method. Whenever an ellipsoidal decomposition is constructed, we draw M points (Θ m , m = 1, 2, . . . , M ) from it as follows: for each draw we first pick an ellipsoid with probability V l / L l=1 V l , where V l are the volumes of the L ellipsoids in the decomposition; a point Θ m is then drawn uniformly from the chosen ellipsoid and we calculate, q m , the number of ellipsoids it lies in. The volume in the union of ellipsoids is then: We note that this Monte Carlo procedure does not require any evaluations of the likelihood function, and thus is not computationally demanding.
4.2. Evidence estimation and posterior samples As an alternative to the canonical NS summation given by Eq. (8) the Bayesian evidence can instead be estimated with reference to the above pseudo-importance sampling density as: Moreover, each one of the N tot points collected by MULTINEST can be assigned the following estimator of its posterior probability density: Since the importance nested sampling scheme does not rely on the ellipsoidal decomposition fully enclosing the region(s) satisfying the constraint L > L i , it can also achieve accurate evidence estimates and posterior summaries from sampling done in the constant efficiency mode of MULTINEST. 3 However, as we discuss shortly, the utility of this feature is often limited by ensuing difficulties in the estimation of uncertainty for such constant efficiency mode evidence estimates. From a computational perspective we note that in a naïve application of this scheme it will be necessary to store N tot points, Θ k , along with the likelihood, L(Θ k ), and prior probability, π(Θ k ), for each, as well as all relevant information describing the ellipsoidal decompositions (centroids, eigen-values and eigen-vectors) at each iteration. Even with a Cholesky factorization of the eigen-vectors, storing the latter may easily result in excessive memory requirements. However, since in the MULTINEST algorithm the prior volume, and consequently the volume occupied by the bounding ellipsoids, shrinks at each subsequent iteration one can confidently assume E i (Θ) = 1 for all points drawn at iterations j > i. At a given iteration then, one needs only to check if points collected from previous iterations lie in the current ellipsoidal decomposition and add the contribution to g(Θ) coming from the current iteration as given in Eq. (12). This results in an enormous reduction in memory requirements as information about the ellipsoidal decomposition from previous iterations no longer needs to be stored.
At each iteration, MULTINEST draws points from the ellipsoidal decomposition, but in order to take account of the volume in the overlaps between ellipsoids, each point is accepted only with probability 1/q where q is the number of ellipsoids in which the given point lies. Rather than discarding all these rejected points, which would be wasteful, we include them by dividing the importance sampling weights as given in Eq. (12), in three components: Assuming that the point Θ was drawn at iteration i, g 1 , g 2 and g 3 are the contributions to importance weight for Θ coming from iteration i, iterations before i and iterations after i respectively. Thus, g 1 is calculated as follows: where q is the number of ellipsoids at iteration i in which point Θ lies, while g 2 is calculated as follows: where V tot,j is volume occupied by the union of ellipsoids at iteration j as given in Eq. (13). Here we have assumed that ellipsoids shrink at subsequent iterations and therefore points drawn at iteration i lie inside the ellipsoidal decompositions of previous iterations as discussed earlier. Finally, g 3 is calculated as follows: 4.3. Evidence error estimation As discussed by Skilling (2004) (and by  for the specific case of MULTINEST) repeated summation of the NS draws under random sampling of the associated X i (governed by t i ; cf. Sec. 3) allows one to estimate the error on the NS evidence approximation from just a single run (whereas many other MC integration techniques, such as thermodynamic integration, require repeat runs to achieve this). Provided that the parameter space has been explored with sufficient thoroughness (i.e., the N live point set has evolved through all the significant posterior modes), the reliability of this evidence estimate was demonstrated in . Importantly, such a single run error estimate can also be calculated for the INS scheme as described below.
Under ordinary (as opposed to pseudo-) importance sampling the unbiased estimator for the asymptotic variance of the evidence estimate here, Var[Ẑ], would be given as follows: withẐ given by Eq. (14). With the draws from MULTINEST representing our a posteriori constructed g(Θ) not in fact an independent, identically distributed sequence from this pseudo-importance sampling function, the above uncertainty estimate is, unfortunately, not strictly applicable here. In particular, with the placement of subsequent ellipses, E j>i , dependent on the position of the live particles drawn up to the present step, i, so too are the subsequently drawn Θ j>i . However, when MULTINEST is run in its default mode, such that we strongly govern the maximum rate at which the volume of the successive E i can shrink we can be confident that our sampling becomes ever more nearly independent and that the dominant variance component is indeed given in Eq. (20). Our reasoning behind this is explained in detail in Appendix B. On the other hand, when MULTINEST is being run in 'constant efficiency mode' we recommended for the user to check (via repeat simulation) that the INS evidence is stable (with respect to its error estimate) for reasonable variation in N live and/or f .

APPLICATIONS
In this section we apply the MULTINEST algorithm with INS described above to three test problems to demonstrate that it indeed calculates the Bayesian evidence much more accurately than vanilla NS. These test examples are chosen to have features that resemble those that can occur in real inference problems in astro-and particle physics.    In this section, we apply MULTINEST with and without INS to sample from a posterior containing multiple modes with pronounced (curving) degeneracies in relatively high dimensions. Our test problem here is the same one used in Allanach and Lester (2008); ; . The likelihood function of this problem is defined as, where In two dimensions, this distribution represents two well separated rings, centred on the points c 1 and c 2 respectively, each of radius r and with a Gaussian radial profile of width w (see Fig. 3).
We investigate the above distribution up to a 50-dimensional parameter space θ. In all cases, the centres of the two rings are separated by 7 units in the parameter space, and we take w 1 = w 2 = 0.1 and r 1 = r 2 = 2. We make r 1 and r 2 equal, since in higher dimensions any slight difference between these two values would result in a vast difference between the volumes occupied by the rings and consequently the ring with the smaller r value would occupy a vanishingly small fraction of the total probability volume, making its detection almost impossible. It should also be noted that setting w = 0.1 means the rings have an extremely narrow Gaussian profile. We impose uniform priors U(−6, 6) on all the parameters. For the two-dimensional case, with the parameters described above, the likelihood is shown in Fig. 3. Table 1 lists the total number of live points (N live ) and target efficiency (f ) used and the total number of likelihood evaluations (N like ) performed by MULTINEST in default and constant efficiency (ceff) modes. The volume of the parameter space increases exponentially with the dimensionality D, therefore we need to increase N live and/or decrease f with D, in order to get accurate estimates of log(Z). The true and estimated values of log(Z) are listed in Table 2.
It can be seen from Table 2 that log(Ẑ) values obtained by MULTINEST with and without INS and in both default and constant efficiency modes are consistent with the true log(Z) for D ≤ 20, the only exception being the log(Ẑ) from constant efficiency mode with INS which is ∼ 6σ away from the analytical log(Z). We attribute this to the heightened potential for underestimation of the INS uncertainties in constant efficiency mode discussed in Sec. 4 and Appendix B. For D ≥ 30 however, the log(Ẑ) values obtained by MULTINEST without INS start to become inaccurate, with constant efficiency mode again giving more inaccurate results as expected. These inaccuracies are caused by inadequate numbers of live points used to cover the region satisfying the constraint L > L i at each iteration i. However, with the same values for N live and f , and indeed with the same set of points,  has a large bias but as discussed in Sec. 4.3, these error estimates are reliable only when the importance sampling distribution is guaranteed to give non-vanishing probabilities for all regions of parameter space where posterior distribution has a non-vanishing probability as well. This is much more difficult to accomplish in a 50D parameter space. In addition to this, the approximations we have made to calculate the volume in the overlapped region of ellipsoids are expected to be less accurate in higher dimensions. Therefore, it is very encouraging that INS can obtain log(Z) to within 0.2 units for a very challenging 50D problem with just 500 live points. We should also notice the number of likelihood evaluations in constant efficiency mode starts to become significantly smaller than in the default mode for D ≥ 20.

Test problem 2: egg-box likelihood
We now demonstrate the application of MULTINEST to a highly multimodal two-dimensional problem, for which the likelihood resembles an egg-box. The un-normalized likelihood is defined as: and we assume a uniform prior U(0, 10π) for both x and y.
A plot of the log-likelihood is shown in Fig. 4 and the prior ranges are chosen such that some of the modes are truncated, making it a challenging problem for identifying all the modes as well as to calculate the evidence accurately. The true value of the log-evidence is log Z = 235.856, obtained by numerical integration on a very fine grid, which is feasible for this simple two-dimensional example.
It was shown in  that MULTINEST can explore the parameter space of this problem efficiently, and also calculate the evidence accurately. Here we demonstrate the accuracy of the evidence obtained with MULTINEST using the INS summation. For low-dimensional problems, results obtained with the constant efficiency mode of MULTINEST agree very well with the ones obtained with the default mode, we therefore only discuss the default mode results in this section.
We use 1000 live points with target efficiency f = 0.5. The results obtained with MULTINEST are illustrated in Fig. 4, in which the dots show the points with the lowest likelihood at successive iterations of the nested sampling process. MULTINEST required ∼ 20, 000 likelihood evaluations and obtained log(Ẑ) = 235.837 ± 0.008 (235.848 ± 0.078) with (without) INS, which compares favourably with the true value given above. In each case, the random number seed was the same, so the points sampled by MULTINEST were identical with and without INS. In order to check if the error estimates on log(Ẑ) are accurate, we ran 10 instances of MULTINEST in both cases, each with a different seed and found the mean and standard deviation of log(Ẑ) to be 235.835 ± 0.009 (235.839 ± 0.063) with (without) INS. In both cases, the standard error agrees with the error estimate from just a single run. There is, however, some indication of bias in the log(Ẑ) value evaluated with INS, which lies ∼ 2σ away from the true value. This is most likely due to the approximations used in calculating the volume in the overlapped region of ellipsoids, as discussed in Sec. 4. Nonetheless, the absolute value of the bias is very low (∼ 0.02), particularly compared with the accuracy (∼ 0.5) to which log-evidence values are usually required in practical applications. 5.3. Test problem 3: 16D Gaussian mixture model Our next test problem is the same as test problem 4 in Weinberg et al. (2013) which is a mixture model of four randomlyoriented Gaussian distributions with their centers uniformly selected from the hypercube [0.5 − 2σ, 0.5 + 2σ] D with D being the dimensionality of the problem and the variance σ 2 of all four Gaussians is set to 0.003. Weights of the Gaussians are distributed according to a Dirichlet distribution with shape parameter α = 1. We impose uniform priors U(0, 1) on all the parameters. The analytical posterior distribution for this problem, marginalized in the first two dimensions is shown in Fig. 5(a).
The analytical value of log(Z) for this problem is 0, regardless of D. We set D = 16 and used 300 live points with target efficiency f = 0.05. The marginalized posterior distribution in the first two dimensions, obtained with the default mode of MULTINEST with INS is shown in Fig. 5(b). The posterior distribution obtained from the constant efficiency mode is identical to the one obtained from the default and therefore we do not show it. In the default mode MULTINEST performed 208, 978 likelihood evaluations and returned log(Ẑ) = −0.03 ± 0.01 (0.39 ± 0.27) with (without) INS. In the constant efficiency mode, 158, 219 likelihood evaluations were performed and log(Ẑ) = 0.21 ± 0.01 (0.25 ± 0.27) with (without) INS.

Test problem 4: 20D Gaussian-LogGamma mixture model
Our final test problem is the same as test problem 2 in Beaujean and Caldwell (2013), in which the likelihood is a mixture model consisting of four identical modes, each of which is a product of an equal number of Gaussian and LogGamma 1D distributions, centred at θ 1 = ±10, θ 2 = ±10, θ 3 = θ 4 = · · · = θ D = 0 in the hypercube θ ∈ [−30, 30] D , where D is the (even) dimensionality of the parameter space. Each Gaussian distribution has unit variance. The LogGamma distribution is asymmetric and heavy-tailed; its scale and shape parameters are both set to unity. We impose uniform priors U(−30, 30) on all the parameters. The analytical marginalised posterior distribution in the subspace (θ 1 , θ 2 ) is shown in Fig. 6(a).
We set D = 20, for which the analytical value of the log-evidence is log(Z) = log(60 −20 ) = −81.887. To be consistent with test problem 3, which is of similar dimensionality, we again used 300 live points with a target efficiency f = 0.05 (note these values differ from those used in Beaujean and Caldwell (2013), who set N live = 1000 and f = 0.3 in the standard vanilla NS version of MULTINEST). The marginalized posterior in the first two dimensions, obtained in the default mode of MULTINEST with INS is shown in Fig. 6(b), and is identical to the corresponding analytical distribution, recovering all four modes with very close to equal weights. The posterior distribution obtained from the constant efficiency mode is identical to the one obtained from the default and therefore we do not show it. In the default mode MULTINEST performed 2,786,538 likelihood evaluations and returned log(Ẑ) = −81.958 ± 0.008 (−78.836 ± 0.398) with (without) INS. In both cases, we see that, for this more challenging problem containing multi-dimensional heavy-tailed distributions, the log-evidence estimates are substantially biased, with each being ∼ 8σ from the true value. Nonetheless, we note that the estimate using INS is much more accurate than that obtained with vanilla NS, and differs from the true value by only ∼ 0.1 units, which is much smaller than the accuracy required in most practical applications. As one might expect, however, the log-evidence estimates obtained in constant

SUMMARY AND DISCUSSION
With the availability of vast amounts of high quality data, statistical inference is increasingly playing an important role in cosmology and astroparticle physics. MCMC techniques and more recently algorithms based on nested sampling have been employed successfully in a variety of different areas. The MULTINEST algorithm in particular has received much attention in astrophysics, cosmology and particle physics owing to its ability to efficiently explore challenging multi-modal distributions as well as to calculate the Bayesian evidence.
In this paper we have discussed further development of the MULTINEST algorithm, based on the implementation of the INS scheme recently proposed by Cameron and Pettitt (2013). INS requires no change in the way MULTINEST explores the parameter space, but can calculate the Bayesian evidence at up to an order-of-magnitude higher accuracy than vanilla nested sampling. Moreover, INS also provides a means to obtain reasonably accurate evidence estimates from the constant efficiency mode of MULTINEST. This is particularly important, as the constant efficiency mode enables MULTINEST to explore higher-dimensional spaces (up to ∼ 50D) much more efficiently than the default mode. Higher evidence accuracy from INS could potentially allow users to use fewer live points N live or higher target efficiency f to achieve the same level of accuracy as vanilla nested sampling, and therefore to speed-up the analysis by several factors. We recommend that users should always check that their posterior distributions are stable with reasonable variation of N live and f . A slight drawback of INS is increased memory requirements. As the importance sampling distributions given in Eqs. 18 and 19 change for every point at each iteration, all the points need to be saved in memory. However, with N live ≤ 1000 the increased memory requirements should be manageable on most modern computers.
Finally, we give some recommendations for setting the number of live points N live and target efficiency f , which determine the accuracy and computational cost of running the MULTINEST algorithm, with or without INS. Generally, the larger the N live and lower the f , the more accurate are the posteriors and evidence values but the higher the computational cost. For multi-modal problems, N live is particularly important as it determines the effective sampling resolution. If it is too small, certain modes, in particular the ones occupying a very small prior mass, can be missed. Experience has shown that the accuracy of evidence is more sensitive to f than N live . In general, for problems where accuracy of evidence is paramount, we suggest f to be no larger than 0.3 in the 'default' mode. In 'constant efficiency mode', we suggest f to be no larger than 0.1 in all cases. Generally, a value of N live in lower hundreds is sufficient. For very low dimensional problems N live can even be in tens. However, for highly multi-modal problems, one may need to set N live to be in a few thousands. It is always advisable to increase N live and reduce f to check if the posteriors and evidence values are stable as function of N live and f .
We review here the heritage of INS amongst the wider family of pseudo-importance sampling, NS, and adaptive Monte Carlo algorithms, for which limited convergence proofs have yet been achieved.
As described in Cameron and Pettitt (2013) the initial idea for INS arose from the study of recursive marginal likelihood estimators, as characterised by Reverse Logistic Regression (RLR ;Geyer 1994;Chen and Shao 1997;Kong et al. 2003) and the Density of States (DoS; Habeck 2012;Tan et al. 2012). In these (equivalent; cf. Cameron and Pettitt 2013) schemes, the marginal likelihood is sought by pooling (or 'losing the labels' on) a series of draws from a pre-specified set of largely unnormalised importance sampling densities, bridging (at least crudely) the prior and posterior; after this a maximum-likelihood-based estimator is used to infer the relative normalisation of each bridging density in light of the 'missing' label information. As emphasised by Kong et al. (2003), these recursive algorithms may, in turn, be seen as deriving from the 'biased sampling' results of Vardi (1985) and collaborators (e.g. Gill et al. 1988), who give consistency and Central Limit Theorem proofs for this deterministic (i.e., non-adaptive), pseudo-importance sampling procedure under simple connectedness/non-separability conditions for the supports of the bridging sequence.
Developed in parallel to the recursive estimators described above, the Deterministic Multiple Mixture Sampling scheme (DMMS; Veach and Guibas 1995;Owen and Zhou 2000) applies much the same strategy, but for a given sequence of strictly normalised importance sampling proposal densities; hence, the motivation for 'losing the labels' here becomes simply the reduction of variance in the resulting estimator with respect to that achievable by allocating the same number of draws to ordinary importance sampling from each proposal separately. [We will discuss further the limiting variance of a simple 'losing the labels' estimator, as relevant to INS, in Appendix B.] At the expense of introducing an intractable (but asymptotically-diminishing) bias, Cornuet et al. (2012) have recently constructed a yet more efficient extension called Adaptive Multiple Importance Sampling (AMIS) in which the sequence of importance sampling proposal densities is refined adaptively towards the target at runtime. In particular, each proposal density after the first is chosen (from a given parametric family) in a manner dependent on the weighted empirical distribution of draws from all previous densities. As suggested by their given numerical examples this approach appears superior to other adaptive importance sampling schemes (e.g. the cross-entropy method; cf. Rubinstein and Kroese 2004) in which the past draws from sub-optimal proposals are ultimately discarded.
In our opinion, despite its genesis in the study of RLR/DoS, INS may perhaps most accurately be viewed as a 'descendent' of this AMIS methodology; the key difference being that INS builds an efficient mixture proposal density for marginal likelihood estimation via the NS pathway (Sec. 3), whereas AMIS aims to iterate towards such a form chosen from within a pre-specified family of proposal densities. In other words, in INS the proposal densities represented by our sequence of ellipsoidal decompositions should share (by design) a near-equal 'importance' in the final mixture, while in AMIS those draws from the earlier proposals are expected to become increasingly insignificant as the later proposals achieve refinement towards their target.
As acknowledged by Cornuet et al. (2012), the inherent dependence structure of the pseudo-importance weighted terms entering the final AMIS summation-owing entirely in this approach to the dependence structure of the corresponding sequence of proposal densities-renders intractable the demonstration of a general consistency for the algorithm. Indeed even the elegant and detailed solution presented by Marin et al. (2012) is reliant on key modifications to the basic procedure, including that the number of points drawn from each successive density grows significantly at each iteration (incompatible with INS; and seemingly at odds too with Cornuet et al.'s original recommendation of heavy sampling from the first proposal), as well as numerous assumptions regarding the nature of the target and proposal families. In light of this historical background we will therefore give particular attention in the following Appendix to the variance reduction benefits of the theoretically-problematic 'losing the labels' strategy as employed in INS, before sketching a rough proof of consistency thereafter (albeit under some strong assumptions on the asymptotic behaviour of the EM plus k-means algorithm employed for ellipsoidal decomposition with respect to the target density; which may be difficult, if not impossible, to establish in practice).
Finally, before proceeding it is worth mentioning briefly the heritage of INS with respect to vanilla NS (Skilling 2004(Skilling , 2006. As described in Sections 3 and 4, the original MULTINEST code ) was designed for estimation of Z via the NS pathway (Skilling 2006) with the challenge of constrained-likelihood sampling tackled via rejection sampling within a series of ellipsoidal decompositions bounding the evolving live point set. In contrast to AMIS (and INS) the convergence properties of the simple NS algorithm are well understood; in particular, Chopin and Robert (2010) have derived a robust CLT for nested sampling, and both Skilling (2006) and Keeton (2011) give insightful discussions. Despite the value of this ready availability of a CLT for vanilla NS, our experience (cf. Sec. 5 of the main text) is that by harnessing the information content of the otherwise-discarded draws from MULTINEST's ellipsoidal rejection sampling the INS summation does ultimately yield in practice a substantially more efficient approximation to the desired marginal likelihood, the reliability of which is moreover adequately estimable via the prescription given subsequently.

B. CONVERGENCE BEHAVIOUR OF INS
In this Appendix we discuss various issues relating to the asymptotic convergence of the INS marginal likelihood estimator (14), which we denote here byẐ INS , towards the true marginal likelihood, Z, as the sample size (controlled by the size of the live point set, N live ) approaches infinity. We begin by considering the intriguing role of pseudo-importance sampling for variance reduction within certain such schemes; this step, ironically, is itself primarily responsible for the intractable bias of the complete algorithm. With this background in mind we can at last outline a heuristic argument for the consistency of INS and consider a break down of its variance into distinct terms of transparent origin.
To be precise, we will investigate here the asymptotic convergence behaviour of the INS estimator with ellipsoidal decompositions almost exactly as implemented in MULTINEST, a detailed description of which is given in the main text (Sections 3 & 4).
For reference, we takeẐ Vtot,i , and c×N live a fixed stopping proxy for the total number of ellipsoidal decompositions required to meet our actual flexible stopping criterion (cf. Sec. 3.2). Each collection of Θ (i) k=1,...,ni here is assumed drawn uniformly from the corresponding ellipsoidal decomposition (of the live particle set), E i (·), with volume, V tot,i , until the discovery of a single point, say Θ This new constrained-likelihood point serves, of course, as the replacement to the L i−1 member of the NS live point set against which the next, E i+1 (·), decomposition is then defined.
The equality in (24) highlights the fact that having pooled our draws from each E i (·) into the pseudo-importance sampling function, g(·), we may proceed to 'lose the labels', (i) , on these as in, e.g., Reverse Logistic Regression or "biased sampling". Note also that we suppose E 1 (·) is fixed to the support of the prior itself (to ensure that the support of L(Θ)π(Θ) is contained within that of g(Θ)), and that we must sample an initial collection of N live live points from the prior as well to populate the original live particle set in advance of our first constrained-likelihood exploration.
Finally, we neglect in the ensuing analysis any uncertainty in our V tot,i since, although these are in fact estimated also via (simple MC) simulation, without the need for likelihood function calls in this endeavour we consider the cost of arbitrarily improving their precision effectively negligible. B.1. Motivation for 'losing the labels' on a normalised pseudo-importance sampling mixture The effectiveness of the so-called 'losing the labels' strategy for marginal likelihood estimation via the recursive pathway can be easily appreciated for the typical RLR/DoS case of multiple unnormalised bridging densities, since by allowing for, e.g., the use of tempered Monte Carlo sampling we immediately alleviate to a large extent the burdens of importance sampling proposal design (cf. Hesterberg 1995). However, its utility in cases of strictly normalised mixtures of proposal densities as encountered in DMMS and INS is perhaps surprising. Owen and Zhou (2000) give a proof that, under the DMMS scheme, the asymptotic variance of the final estimator will not be very much worse than that achievable under ordinary importance sampling from the optimal distribution alone. However, as the INS sequence of ellipsoids is not designed to contain a single optimal proposal, but rather to function 'optimally' as an ensemble we focus here on demonstrating the strict ordering (from largest to smallest) of the asymptotic variance for (I) ordinary importance sampling under each mixture component separately, (II) ordinary importance sampling under the true mixture density itself, and (III) pseudo-importance sampling from the mixture density (i.e., 'losing the labels').
Consider a grossly simplified version of INS in which, at the nth iteration (it is more convenient here to use n rather than i as the iteration counter), a single random point is drawn independently from each of n labelled densities, h k,n (·) (k = 1, 2, . . . , n), with identical supports matching those of the target, f (·). We denote the resulting set of n samples by Θ (n) k=1,2,...,n . The three key simplifications here are: (I) that the draws are independent, when in MULTINEST they are inherently dependent; (II) that the supports of the h k,n (·) match, when in fact the ellipsoidal decompositions, E n (·), of MULTINEST have generally nested supports (though one could modify them appropriately in the manner of defensive sampling; Hesterberg 1995); and (III) that a single point is drawn from each labelled density, when in fact the sampling from each E n (Θ) under MULTINEST follows a negative binomial distribution for one E n (Θ) ∩ {L(Θ) > L n−1 } 'success'. Suppose also now that the unbiased estimator, Z , for the normalizing constant belonging to f (·), namely Z = f (Θ)dΘ, in such single draw importance sampling from each of the specified h k,n (·) has finite (but non-zero) variance (cf. Hesterberg 1995), i.e., and that together ourẐ (n) k satisfy Lindeberg's condition such that the CLT holds for this triangular array of random variables (cf. Billingsley 1995). Now, if we would decide to keep the labels, k, on our independent draws from the sequence of h k,n (·) then supposing no prior knowledge of any σ 2 k,n (i.e., no prior knowledge of how close each proposal density might be to our target) the most sensible option might be to take as a 'best guess' for Z the (unweighted) sample mean of our individualẐ With a common mean and finite variances for each, this sum over a triangular array converges (in distribution) to a univariate Normal with mean, Z, and variance, here we use the abbreviation, s 2 n = n k=1 σ 2 k,n . On the other hand, if we would instead decide to lose the labels on our independent draws we might then follow Vardi's method and imagine each Θ (n) k to have come from the mixture distribution, g(Θ) = 1 n n j=1 h j,n (Θ), for which the alternative estimator,Ẑ may be derived. To see that theẐ unlabelled estimator so defined is in fact unbiased we letẐ For n iid samples drawn faithfully from the mixture density, g(·), we would expect (via the CLT) that the estimatorẐ unlabelled will converge (in distribution) once again to a univariate Normal with mean, Z, but with alternative variance, σ 1 n f (Θ) 2 /g(Θ)dΘ − Z 2 . However, for the pseudo-importance sampling from g(·) described above, in which we instead pool an explicit sample from each of its separate mixture components, the asymptotic variance ofẐ unlabelled is significantly smaller again. In particular, with equality achieved only in the trivial case that all mixture components are identical. The variance reduction here in the pseudoimportance sampling framework relative to the true importance sampling case derives of course from the effective replacement of multinomial sampling of the mixture components by fixed sampling from their expected proportions. Comparing now the asymptotic variances of our labelled and unlabelled estimators we can see that the latter is (perhaps surprisingly) always smaller than the former, i.e., (recalling that all densities here are strictly positive, of course); thus, we observe the ordering, This is, as has been remarked in the past, the paradox of the 'losing the labels' idea; that by throwing away information about our sampling process we appear to gain information about our target! In fact, however, all we are really doing by choosing to estimate Z withẐ unlabelled rather thanẐ labelled is to use the information we have extracted from f (·) in a more efficient manner, as understood (from the above error analysis) a priori to our actual importance sampling. The strict ordering shown here explains why we have selected a pseudo-importance sampling strategy for combining the ellipsoidal draws in MULTINEST as opposed to, e.g., modifying our sampling from the E n (·) to match (defensively) the support of π(Θ) and compiling estimatorŝ Z k separately-though the latter would simplify our convergence analysis the former should be (asymptotically) much more efficient.
B.2. Consistency of INS To establish a heuristic argument for consistency of the INS scheme we must first consider the nature of the limiting distribution for the sequence of ellipsoidal decompositions, {E i (·)}, as N live → ∞. To do so we introduce the following strong assumption: that for the constrained-likelihood contour corresponding to each X i ∈ [0, 1] there exists a unique, limiting ellipsoidal decomposition, E * Xi (·), to which the MULTINEST algorithm's E i (·) will converge (if not stopped early) almost surely for all N live for which X i = exp(−i/N live ) for some i ∈ {1, 2, . . . , c × N live }. In particular, we suppose that both the design of the EM plus k-means code for constructing our E i (·) and the nature of the likelihood function, L(Θ), are such for any given > 0 there is an N live large enough that thereafter Another supposition we make is that the limiting family of ellipsoidal decompositions, {E * Xi (·)}, is at least left or right 'continuous' in the same sense at every point of its rational baseline; i.e., for each X i and any > 0 there exists a δ > 0 such that X i − X j < δ and/or X j − X i < δ implies Various conditions for almost sure convergence of EM (Wu 1983) and k-means (Pollard 1981) algorithms have been demonstrated in the past, but we have an intractable dependence structure operating on the E k (·) for INS and it is not at all obvious how to clearly formulate such conditions here. The complexity of this task can perhaps most easily be appreciated by considering the limited availability of results for the convergence in volume of random convex hulls from uniform sampling within regular polygons in high-dimensions (e.g. Schneider and Wieacker 1980;Schneider 2008). On the other hand, we may suspect that the necessary conditions for the above are similar to those required in any case for almost sure 'coverage' of each constrainedlikelihood surface by its corresponding ellipsoidal decomposition; the latter being an often ignored assumption of rejection NS. That is, even for a generous dilation of the simple proposal ellipsoids, as suggested by Mukherjee et al. (2006), one can easily identify some family of (typically non-convex) L(Θ) for which the given dilation factor will be demonstrably insufficient; though whether such a 'pathological' L(Θ) would be likely to arise in standard statistical settings is perhaps another matter entirely! The necessity of these assumptions, and in particular our second regarding the 'continuity' of the limiting {E * Xi (·)}, is two-fold: to ensure that a limiting distribution exists (this echoes the requirement for there to exist an optimal proposal in the equivalent AMIS analysis of Marin et al. 2012), and to ensure that its form is such as to render irrelevant the inevitable stochastic variation and bias in our negative binomial sampling of n i points from E i (·). Important to acknowledge is that not only is the number of points drawn from each ellipsoidal decomposition, n i , a random variable, but the collection of n i − 1 draws in E i (·) ∩ {L(Θ) < L i−1 } plus one in E i (·)∩L(Θ) > L i−1 from a single realization cannot strictly be considered a uniform draw from E i (·), though we treat it as such in our summation for g(Θ). Indeed the expected proportion of these draws represented by the single desired L(Θ) > L i−1 point, namely E[ 1 ni ] = −pi log pi 1−pi , does not even match its fraction of π(Θ) by 'volume', here p i . Our argument must therefore be that (asymptotically) with more and more near-identical ellipsoids converging in the vicinity of each E * i (·) as N live → ∞ the pool of all such biased draws from our constrained-likelihood sampling within each of these nearby, near-identical ellipsoids ultimately approximates an unbiased sample, and that the mean number of draws from these will approach its long-run average, say n * i . With such convergence towards a limiting distribution, F * (Θ), defined by the set of pairings, {E * i (·), n * i }, supposed it is then straightforward to confirm that via the Strong Law of Large Numbers the empirical distribution function of the samples Θ k drawn under INS, F INS (Θ), will converge in distribution to this asymptotic form; the convergence-determining classes here being simply the (lower open, upper closed) hyper-cubes in the compact subset, [0,1] and thus and hence (with the corresponding Var[Ẑ INS ] → 0) the consistency ofẐ INS .
B.3. Variance breakdown of INS Given the evident dependence of the INS variance on three distinct sources-(I) the stochasticity of the live point set, and its decompositions, {E i (·)}; (II) the negative binomial sampling of the n i ; and (III) the importance sampling variance of the drawn f (Θ k )/g(Θ k )-it makes good sense to break these components down into their contributory terms via the Law of Total Variance as follows: Now the first term represents explicitly the variance contribution from the inherent randomness of the ellipsoidal decomposition sequence, {E i (·)}, which we might suppose negligble provided the geometric exploration of the posterior has been 'sufficiently thorough', meaning that the N live point set has evolved through all the significant posterior modes. The second and third terms represent the negative binomial sampling and 'ordinary' importance sampling variance contributions expected under the distribution of {E i (·)}. With the realised {E i (·)} being, of course, an unbiased draw from its parent distribution any unbiased estimator of these two additional variance components applied to our realised {n i } and {Θ k } could be considered likewise unbiased. However, no such estimators are readily available, hence we pragmatically suppose the second term also negligble and make do with the 'ordinary' importance sampling estimator, given by Eq. (20), for the third term.
Acknowledging the possibility for under-estimation of the INS variance in this way it becomes prudent to consider strategies for minimising our unaccounted variance contributions. The first, suggested by our preceeding discussion of asymptotic consistency for the INS, is to maximise the size of the live point set used. Of course, whether for INS or vanilla NS with MULTINEST we have no prescription for the requisite N live , and the range of feasible N live will often be strongly limited by the available computational resources. Hence we can give here only the general advice of cautious application; in particular it may be best to confirm a reasonable similarly between the estimated variance from Eq. (20) above and that observed from repeat simulation at an N live of manageable runtime prior to launching MULTINEST at a more expensive N live . The second means of reducing the variance in our two unaccounted terms is to stick with the default mode of MULTINEST, rather than opt for 'constant efficiency mode', since by bounding the maximum rate at which the ellipsoidal decompositions may shrink towards the posterior mode we automatically reduce the variance in the random variable, {E i (·)}, and that of {n i } andẐ conditional upon it! C. SOME MEASURE-THEORETIC CONSIDERATIONS When outlining in Sec. 3 the transformation from integration of L(Θ) over π(Θ)dΘ to integration of L(X) over dX (the prior mass cumulant) at the heart of the NS algorithm, we elected, in the interest of simplicity, to omit a number of underlying measure-theoretic details. The significance of these are perhaps only of particular relevance to understanding the use of the NS posterior weights, L i w i /Ẑ from Eq. (10), for inference regarding functions of Θ (e.g. its first, second, and higher-order moments) with respect to the posterior density. However, as this issue has been raised by Chopin and Robert (2010) and we feel that their Lemma 1 deserves clarification we give here a brief measure-theoretic formulation of NS with this in mind.
As with many Bayesian inference problems we begin by supposing the existence of two well-defined probability spaces: (I) that of the prior with (normalised) probability measure, P π , defined for events in the σ-algebra, Σ Θ , of its sample space, Ω Θ , i.e., (P π , Ω Θ , Σ Θ ), and (II) that of the posterior with measure P π defined on the same space, i.e., (P π , Ω Θ , Σ Θ ). Moreover, we suppose that each may be characterised by its Radon-Nikodym derivative with respect to a common σ-finite baseline measure on a complete, separable metric space; that is, P π (A ∈ Σ Θ ) = A π(Θ){dΘ} and P π (A ∈ Σ Θ ) = A π(Θ)L(Θ)/Z{dΘ}. NS then proposes to construct the induced measure, P X , on the σ-algebra, Σ X , generated by the Borel sets of the sample space, Ω X = [0, 1] ∈ R, and defined by the transformation, X : Ω Θ → Ω X with X(Θ ) = {Θ:L(Θ)>L(Θ )} π(Θ){dΘ}. The validity of which requires only the measurability of this transformation (i.e., X −1 B ∈ Σ Θ for all B ∈ Σ X ); e.g. in the metric space R k with reference Lebesgue measure, the almost everywhere continuity of L(Θ). However, for the proposed Riemann integration of vanilla NS to be valid we will also need the induced P X to be absolutely continuous with respect to the Lebesgue reference measure on [0, 1], such that we can write P X (B ∈ Σ X ) = B L(X)/Z{dX}. The additional condition for this given by Chopin and Robert (2010) is that L(Θ) has connected support. To state the objective of vanilla NS in a single line: if we can compute L(X) we can find Z simply by solving for holds whenf (X) is absolutely continuous. We agree that this is true and Chopin and Robert (2010) give a valid proof in their Appendix based on the Monotone Convergence Theorem. However, we suggest that given the already supposed validity of the measure P X , and its Radon-Nikodym derivative with respect to the reference Lebesgue measure, {dX}, upon which NS is based, the equality of the above relation is already true without absolute continuity via the change of variables theorem (Halmos 1950), in the sense that wherever one side exists the other exists and will be equal to it. One trivial example of a discontinuous f (X) for which both sides of 36 exist and are equal is that induced by the indicator function for L(Θ) > X * . To see that thef (X) corresponding to a given, measurable f (Θ) has the stated interpretation as a conditional expectation we observe that E π [f (Θ)|L(Θ) = L(X)] may be written as a function of X, using the interpretation of conditional distributions as derivatives (Pfanzagl 1979). Thus,