TL;DR
This paper investigates the sources of sampling errors in nested sampling parameter estimation, introduces a new diagrammatic representation, and proposes a novel algorithm to accurately measure these errors, validated through empirical tests.
Contribution
It identifies the main sources of sampling errors in nested sampling parameter estimation and introduces a new algorithm for accurately measuring these errors in single runs.
Findings
Identified two main sources of sampling errors in nested sampling.
Current methods cannot accurately measure errors in a single run.
Proposed a new algorithm validated through empirical verification.
Abstract
Sampling errors in nested sampling parameter estimation differ from those in Bayesian evidence calculation, but have been little studied in the literature. This paper provides the first explanation of the two main sources of sampling errors in nested sampling parameter estimation, and presents a new diagrammatic representation for the process. We find no current method can accurately measure the parameter estimation errors of a single nested sampling run, and propose a method for doing so using a new algorithm for dividing nested sampling runs. We empirically verify our conclusions and the accuracy of our new method.
| Repeated runs | |||
|---|---|---|---|
| Bootstrap / Repeats | |||
| Simulated / Repeats | |||
| Bootstrap estimate variation | |||
| Simulated estimate variation | |||
| Bootstrap | |||
| Bootstrap Mean coverage | |||
| Bootstrap coverage |
| Repeated runs | |||
|---|---|---|---|
| Bootstrap / Repeats | |||
| Simulated / Repeats | |||
| Bootstrap estimate variation | |||
| Simulated estimate variation | |||
| Bootstrap | |||
| Bootstrap Mean coverage | |||
| Bootstrap coverage |
| Repeats | |||
|---|---|---|---|
| Split into 20 / Repeats | |||
| Bootstrap / Repeats | |||
| Split into 20 estimate variation | |||
| Bootstrap estimate variation |
| Repeated runs | |||
|---|---|---|---|
| Bootstrap / Repeats | |||
| Simulated / Repeats | |||
| Bootstrap estimate variation | |||
| Simulated estimate variation | |||
| Bootstrap | |||
| Bootstrap Mean coverage | |||
| Bootstrap coverage |
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Code & Models
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Sampling errors in nested sampling parameter estimation
Edward Higsonlabel=e1][email protected] [
Will Handleylabel=e2][email protected] [
Mike Hobsonlabel=e3][email protected] [
Anthony Lasenbylabel=e4][email protected] [ Astrophysics Group, Battcock Centre, Cavendish Laboratory, JJ Thomson Avenue, Cambridge CB3 0HE, UK
Kavli Institute for Cosmology, Madingley Road, Cambridge, CB3 0HA, UK
(0000)
Abstract
Sampling errors in nested sampling parameter estimation differ from those in Bayesian evidence calculation, but have been little studied in the literature. This paper provides the first explanation of the two main sources of sampling errors in nested sampling parameter estimation, and presents a new diagrammatic representation for the process. We find no current method can accurately measure the parameter estimation errors of a single nested sampling run, and propose a method for doing so using a new algorithm for dividing nested sampling runs. We empirically verify our conclusions and the accuracy of our new method.
nested sampling,
parameter estimation,
doi:
0000
keywords:
††volume: 00††issue: 0
\sidecaptionvpos
figurec
, ,
and
t1 t2 t3 t4
1 Introduction
Nested sampling (Skilling, 2006) is a Monte Carlo method for Bayesian analysis which simultaneously calculates both Bayesian evidences and posterior samples. The early development of the algorithm was focused on evidence calculation, which is computationally expensive using variants of standard Markov chain Monte Carlo (MCMC) sampling based on the Metropolis-Hastings algorithm (MacKay, 2003).
Contemporary implementations such as MultiNest (Feroz and Hobson, 2008; Feroz et al., 2009, 2013) and PolyChord (Handley et al., 2015a, b) are now also extensively used for parameter estimation from posterior samples (see for example Planck Collaboration, 2016). Nested sampling compares favourably to MCMC-based parameter estimation for degenerate, multi-modal likelihoods as it has no “thermal” transition probability and exponentially compresses the prior distribution to the posterior. However, despite its increasing popularity, the sampling errors in nested sampling parameter estimation are poorly understood.
Correctly quantifying uncertainty is vital for identifying spurious results — in particular we find sampling errors often significantly affect estimates of credible intervals on parameters. Conversely, finding such errors are very small may imply an unnecessarily large amount of computational resource is being used for the calculation. This paper has two goals: to provide an explanation of the sources of these errors and an empirical technique for estimating them. One obvious method is to repeat the analysis some times, although this increases the computational cost by a corresponding factor. Interestingly, we find no current method can accurately estimate these errors on parameter estimates from a single analysis, and so we present a new method for doing this. Our approach uses a new algorithm for dividing a single nested sampling run into multiple valid nested sampling runs; these can then be recombined in different combinations using resampling techniques such as the bootstrap. We test our results and new method empirically.
The paper begins with background on sampling errors in parameter estimation from posterior samples, then describes the nested sampling algorithm and how it is currently used for parameter estimation in Section 2. We explain the two main sources of sampling errors in nested sampling parameter estimation in Section 3, and present a new diagrammatic representation of the process (illustrated in Figures 3(a), 3(c), 3(b), 3(d) and 3(e)). Section 4 describes our new method for measuring sampling errors from a single nested sampling run, using our new algorithm for division of such runs.
We empirically test our method’s accuracy in Section 5 with the help of analytical cases in the manner described by Keeton (2011). Here one can obtain uncorrelated samples from the prior space within some likelihood contour using standard techniques, and we term the resulting procedure perfect nested sampling. Our approach accurately quantifies uncertainties on parameter estimates from the stochasticity of the nested sampling algorithm, but software used for practical problems may produce additional errors from correlated samples within likelihood contours that are specific to a given implementation. We discuss implementation-specific errors in Section 6, including testing sampling error estimates from our method for PolyChord calculations. Our method gives superior performance to the current approach and can be easily included in nested sampling software; we are currently working on incorporating it into future versions of PolyChord.
Background: sampling errors in parameter estimation
Sampling can be used to represent a posterior distribution via a set of weighted samples
[TABLE]
where each is drawn from the posterior distribution with probability proportional to , and . Likelihoods are often computationally expensive functions, so the goal of parameter estimation is to sample the posterior distribution numerically with a limited number of likelihood calls.
Samples may be used to compute numerical results. For example, the posterior expectation of a function of the parameters can be estimated as
[TABLE]
In this case the sampling error is the difference between and the exact value of . Often the posterior distributions of parameters are of interest, and are estimated numerically from the samples by dividing the parameter space into cells or via kernel density estimation.
There have been many works on approximating MCMC sampling errors, including investigation of quantiles and the amount of computation required to reach some level of accuracy — see for example Doss et al. (2015); Flegal et al. (2008); Liu et al. (2016). In particular Sequential Monte Carlo samplers (Del Moral et al., 2006) have similarities with nested sampling, and their sampling errors are better understood. For some related methods such as the Tootsie Pop algorithm (Huber and Schott, 2014) and accelerated simulated annealing (Bezáková et al., 2008) the error distribution is known exactly, although these techniques are less widely used. This paper introduces empirically tested techniques for quantifying sampling errors from the nested sampling algorithm.
2 The nested sampling algorithm
Nested sampling (Skilling, 2006) is a numerical method computing Bayesian evidences
[TABLE]
and samples from the posterior distribution
[TABLE]
given some likelihood and prior .
Initially points, termed live points, are sampled randomly from the prior. At each iteration , the live point with the lowest likelihood is removed and replaced by a new live point sampled from the prior subject to the constraint that it has a likelihood higher than . Iterating until some termination condition is met generates a list of discarded samples known as dead points, which are used to estimate the evidence and make posterior inferences111The remaining live points at termination can also be used if required, but termination conditions can be chosen such that this makes a negligible difference to calculation results.. We refer to the completed nested sampling process as a run.
To compute the evidence, the many-dimensional integral (3) is reduced to a one-dimensional integral in terms of the fractional prior volume within an iso-likelihood contour. We define the fraction of the prior with likelihood greater than some value as , where
[TABLE]
and . Provided the inverse exists222A sufficient condition for to exist is for to be continuous and to have a connected support. See Chopin and Robert (2010) and Feroz et al. (2013, Appendix C) for a more detailed measure-theoretic discussion., the evidence (3) can be expressed as
[TABLE]
Given a set of dead points with likelihoods , the corresponding prior volumes are unknown but are modelled statistically as , where and each shrinkage ratio is independently distributed as the largest of random variables from the interval (Skilling, 2006). Hence:
[TABLE]
and the algorithm samples within an exponentially shrinking part of the prior. This exponential shrinkage is shown schematically in Figure 1.
2.1 Evidence estimation
Nested sampling therefore allows one to approximate the evidence (6) via a quadrature sum over the dead points
[TABLE]
where are the unknown set of shrinkage ratios for the iterations of the nested sampling process, and each is an independent random variable drawn from distribution (7). The shrinkage ratios define the prior volumes via , and the are appropriately chosen quadrature weights roughly corresponding to the volume of the “prior shell” to which a given dead point belongs. For example, using the trapezium rule: .
Given that the shrinkage ratios are a priori unknown, we may quantify our knowledge of by simulating sets of according to (7), and working with the distribution of the resulting set of evidences from (8) (Skilling, 2006). Typically one then computes and reports a mean value and error for from this distribution.
Several alternative methods for calculating evidence inferences are reported in the literature. Skilling (2006) also proposes an error calculation based on relative entropy, which demonstrates that the uncertainty of is dominated by the Poisson variability in the number of steps required to reach the bulk of the posterior mass. Keeton (2011) uses distribution moments and running totals which are updated with each nested sampling step. This method has been extended by Handley et al. (2015b) to allow the splitting of multi-modal likelihoods into different clusters and the treatment of variable numbers of live points. For a more detailed discussion of the convergence properties of nested sampling evidences, see Chopin and Robert (2010).
Thus, the dominant sampling error in the evidence estimate (8) from perfect nested sampling is from statistical variation in the unknown volumes of the prior “shells” that each point represents. The error from approximating the integral for with a sum can be safely neglected unless is very small333The trapezium rule error is , and if required other methods such as Simpson integration could be used. (Skilling, 2006). There is also some error from terminating the algorithm and truncating the sum, but this is can be made negligible with appropriate termination conditions.
2.2 Parameter estimation
One may also perform posterior inference from nested sampling by using the dead points to construct a set of posterior samples with weights proportional to their share of the posterior mass (Skilling, 2006):
[TABLE]
As before, is the set of prior shrinkage ratios and in the trapezium rule case .
The weights defined by (9) present a departure from traditional sampling approaches in that the are random variables, with their stochasticity determined by (7). When computing expectations (2) there is now an additional error associated with our lack of knowledge of the precise values . Nested sampling software packages such as MultiNest and PolyChord produce posterior files containing only the expected values
[TABLE]
To account for the stochasticity in the weights , Skilling (2006) suggests simulating the prior volume shrinkage ratios in the same manner as for evidence estimation, and using these simulations to calculate a set of values for estimators such as (2). The sampling error should then be estimated from the variation within this sample; we term this the simulated weights method. We believe this procedure is the only estimate of sampling errors in parameter estimation from a single nested sampling run proposed in the literature. However it is in general an underestimate, as can be seen in the numerical tests in Section 5. Section B of the supplementary material discusses this underestimation of errors in detail.
We now describe why the simulated weights method does not capture all sources of sampling errors, and in Section 4 we propose a new method for correctly computing these errors.
3 Sources of sampling errors in nested sampling parameter estimation
In order to understand why the simulated weights method underestimates sampling errors, we require a result from Chopin and Robert (2010). They show that the expectation integral (2) may be re-phrased in terms of the prior volume via:
[TABLE]
where is the prior expectation of on some iso-likelihood contour ,
[TABLE]
The simulated weights approach amounts to discretising the integral (11) as
[TABLE]
and, most importantly, further requiring that we may use as a proxy for at each point , In some special cases for all and this approach is valid, for example when (entropy computation), but in general it is not. This can cause significant inaccuracies as iso-likelihood contours often span wide ranges of different parameter values, as illustrated in Figure 2.
To summarise, the dominant sampling errors in estimating some parameter or function of parameters from perfect nested sampling typically come from two sources:
- (i)
approximating the unknown prior volumes with their expectation using (7); 2. (ii)
approximating the mean value of a function of parameters over an entire iso-likelihood contour with its value at a single point .
Errors from (i) are also present in evidence calculation; in the parameter estimation case they are typically smaller as results depend only on the relative weights of the samples. In contrast (ii) is only present in parameter estimation, where it is typically a significant or dominant source of sampling errors. The relative contributions of (i) and (ii) are empirically tested in Section A of the supplementary material, where they are calculated for analytical cases by using exact values for weights and by replacing with . The simulated weights method underestimates sampling errors in parameter estimation as it ignores errors from (ii).
We now introduce a new diagrammatic representation of nested sampling parameter estimation to illustrate the two different sources of sampling errors.
3.1 Diagrammatic representation
Nested sampling transforms evidence calculations of any dimension into a 1-dimensional problem444For practical nested sampling problems implementation-specific errors can differ for two likelihoods with the same . For example if one likelihood has a much higher dimension and a much larger number of modes than the other it may have larger errors from the implementation software. in which can be entirely represented on a diagram like Figure 1. An analogous diagram for parameter estimation must also illustrate sampling a single point on each iso-likelihood contour from the distribution .
We propose a generalisation of Figure 1 for visualising parameter estimation problems, and present it in Figures 3(a), 3(c), 3(b), 3(d) and 3(e). The top panel in each figure is similar to Figure 1 and shows the relative posterior mass at each value of . The lower central panel shows the probability distribution and its mean . The posterior distribution is shown on the left — this is equal to the distributions (the lower central panel) marginalised over in proportion to the posterior weight at each (the top panel).
For these example plots we use -dimensional spherical unit Gaussian likelihoods
[TABLE]
and -dimensional spherical unit Cauchy likelihoods
[TABLE]
with -dimensional co-centred spherical Gaussian priors
[TABLE]
We denote the first component of the vector as , although by symmetry the results will be the same for any component. and are the first and second moments of the posterior distribution of .
The form of the distribution as varies depends on the likelihood only through the shape of the iso-likelihood contours . Therefore the lower central panel of the diagrams for some is the same for any likelihoods with the same contours — this can be seen in Figures 3(a) and 3(b), where the differences in the posterior (left panel) are due only to the different posterior weights in (top panel).
These diagrams can be constructed for any nested sampling calculation using the posterior samples and kernel density estimation; this could provide insight into the nature of the calculation and the relative contributions from the different sources of sampling errors.
3.2 Transforming a parameter estimation problem into 2 dimensions
As illustrated by our diagrams, nested sampling parameter estimation is fundamentally a 2-dimensional problem in and . In fact a parameter estimation calculation for some given is equivalent to a 2-dimensional problem for given when
[TABLE]
for all . Any transformation satisfying (17) and (18) will leave our proposed diagram for the calculation unchanged. Parameter estimation can also be represented as a 1-dimensional problem in combined with a univariate stochastic process for each dead point with the distribution .
One way to express a general nested sampling calculation in 2 dimensions is to map it onto the unit square with uniform priors and a likelihood which is independent of and satisfies (17). In this case is as before the remaining fractional prior volume and parameterises each iso-likelihood contour. Using inverse transform sampling, for a general a corresponding satisfying (18) is
[TABLE]
where is the inverse of the cumulative distribution
[TABLE]
As an example let us consider -dimensional spherically symmetric likelihoods such as (14) or (15) with co-centred spherically symmetric priors such as (16). Then is a function only of the radial distance from the centre , and the iso-likelihood contours are hyperspherical shells of some radius . The probability distribution of a single parameter (a single component of ) on such an iso-likelihood contour is then
[TABLE]
can be sampled directly or used to calculate the inverse cumulative distribution which together with knowledge of the function allows the parameter estimation of a -dimensional Gaussian to be transformed into a 2-dimensional problem on the unit square.
Samples from (21) can be generated efficiently using the symmetry around and the change of variables to give a Beta distribution
[TABLE]
This technique is used for the numerical tests in Section 5, and allows the efficient sampling of high dimensional spherically symmetric distributions where only a few parameters are of interest without generating all the remaining uninteresting parameters.
4 Estimating sampling errors in nested sampling parameter estimation
Following the discussion of sources of sampling errors in Section 3, we seek a method for correctly calculating parameter estimation sampling errors from a single nested sampling run. As no additional samples are available, a natural starting point is to utilise resampling techniques such as the jackknife (Tukey, 1958), bootstrap (Efron, 1979) and Bayesian bootstrap (Rubin, 1981), which estimate the uncertainty on inferences from a set of samples by calculating the variation when samples are re-weighted.
However, as described in Section 2.2, the uncertainty in nested sampling weights produces additional sampling errors which are unique to the nested sampling process. These are not accounted for by naïvely applying jackknives and bootstraps to posterior samples produced by nested sampling, and these approaches fail when tested numerically. We instead require a method for dividing runs in a manner that preserves the statistical properties of nested sampling. No such method exists in the literature, so we present one in the remainder of this section.
4.1 Dividing runs into threads
Skilling (2006) describes how several nested sampling runs with live points may be combined simply by merging the dead points and sorting by likelihood value. The combined sequence of dead points is equivalent to a single nested sampling run with live points.
In fact, as we show now, the reverse procedure is also possible. A nested sampling run with points can be unwoven into a set of valid nested sampling runs, each with . We term these single live point runs threads. During nested sampling, each dead point is replaced by a new point sampled uniformly within its iso-likelihood contour . Starting from each initial live point that is generated, one may follow this sequence of replacements down the set of dead points. This sub-sequence of dead points is in fact a nested sampling run with . More formally:
A few points are worthy of note:
splitting a run by randomly selecting some fraction of the dead points will not produce threads (i.e. single point nested sampling runs); 2. 2.
one may split a given nested sampling run into separate runs with by first separating into threads, and then recombining threads as desired; 3. 3.
the algorithm can be easily adapted for varying numbers of live points by permitting it to select multiple points on contours where increases. This can result in constituent threads stopping or dividing into multiple threads part of the way through the run; 4. 4.
typically there is only one point which was sampled uniformly from the prior volume within each dead point ’s iso-likelihood contour — the point which replaced . A sufficient condition for a nested sampling run to only have one unique division into threads is that is an injective function; 5. 5.
in order for the threads to be true nested sampling runs, care must be taken with the termination conditions conditions used. See Section D of the supplementary material for a full discussion.
Given that threads represent independent nested sampling runs, one may apply standard resampling techniques to the set of threads and approximate the entire sampling error distribution without making assumptions about its form. This works as the values of the dead points from some run with live points form a Poisson process with rate , meaning the values of the dead points of a single thread are a Poisson process of rate 1. For typical problems with computationally expensive likelihoods the computational cost of even a large number of resampling replications is negligible.
Having introduced a framework for applying resampling to nested sampling parameter estimation we now present an example method using bootstrap resampling.
4.2 Bootstrap estimate of sampling errors
Given observations , the bootstrap (Efron, 1979) creates new data sets by drawing samples from with replacement. This corresponds to approximating the probability distribution of a single data point as
[TABLE]
where is the Dirac delta function (Ivezić et al., 2014).
As the form of the distribution of sampling errors for a general nested sampling parameter estimation problem is not known, we use the non-parametric bootstrap. In this case the uncertainty on a quantity calculated from the data can be estimated by calculating for a number of resampled data sets . For example the bootstrap estimate of the standard error on is
[TABLE]
There are many methods for calculating approximate credible intervals on from bootstrap replications — see Efron and Tibshirani (1986) for a detailed discussion. A simple approach from Johnson (2001) is to estimate the boundaries of the and credible regions555If the distribution of bootstrap replications is skewed then the implied probability distribution of is skewed in the opposite direction, as can be seen from (26) and (27). See Loredo (2012, Section 2) for a discussion. as
[TABLE]
where is the inverse cumulative distribution of the bootstrap samples . is typically sufficient for an estimate of the standard deviation of a parameter estimate due to sampling errors, but depending on the method used credible intervals on parameter estimates may require 1,000 bootstrap replications or more (Efron and Tibshirani, 1986).
When the bootstrap is applied to nested sampling each observation is a thread, and the number of observations is . Calculating the quantity involves first combining the set of threads into a single run using Skilling (2006)’s method (described in Section 4.1), then performing a standard nested sampling calculation including estimating the weight of each point statistically. Including the same thread multiple times does not cause problems — repeated dead points are simply assigned the weights and respectively.
The following algorithm provides a set of bootstrap replications and an estimate of the standard deviation of sampling errors.
We find that bootstrap resampling gives better results than jackknife resampling, which fails to calculate sampling errors on credible intervals of posterior distributions of parameters such as . The Bayesian bootstrap was not used as it gives each observation a non-integer weight, which requires modifying nested sampling’s use of dead points to statistically estimate prior volume shrinkages.
Resampling techniques such as the bootstrap can generate many simulated runs with the same number of live points as the original run. In comparison sampling error estimates from simply splitting a run into many smaller runs and assessing their variation perform poorly, as shown in Section C of the supplementary material.
5 Numerical tests
Following Keeton (2011) we first test our new method using analytic cases where uncorrelated samples can be easily obtained from the prior within an iso-likelihood contour, allowing us to perform perfect nested sampling. This ensures our results are not affected by imperfect implementation of the nested sampling algorithm by a specific software.
As discussed in Section 3 perfect nested sampling parameter estimation problems depend on the likelihood and prior only through the distribution of posterior mass and the distribution of parameters on iso-likelihood contours , both of which are functions of both and . We therefore empirically test our method using a wide range of distributions of posterior mass, and examine several functions of parameters in each case. We construct such tests using Gaussian likelihoods (14) and Cauchy likelihoods (15) of a variety of dimensions , each with a Gaussian prior (16). The different distributions of posterior mass for different are illustrated in Figure 4; the Cauchy distributions have extremely fat tails with significant sample weights throughout the range of values explored.
We use the termination conditions described by Handley et al. (2015b, Section 3.4), stopping when the estimated evidence contained in the live points is less than times the evidence contained in dead points (see Section D of the supplementary material for a discussion of termination conditions for nested sampling parameter estimation). Numerical calculations for high-dimensional cases are performed in two dimensions using the technique described in Section 3.2.
As in Section 3 we denote the first component of the vector as , although by symmetry the results will be the same for any component. and are the first and second moments of the posterior distribution of , and the one-tailed upper credible interval is the value for which .
5.1 3-dimensional Gaussian example
We first test our bootstrap approach to estimating sampling errors on a 3-dimensional Gaussian likelihood (14) — Figure 5 illustrates sampling errors on posterior distributions of parameters in this case. Unlike the simulated weights method, the mean estimates of sampling errors from our method are very close to measurements of sampling errors from repeated calculations — this shown in the second row of Table 1. Furthermore the fractional variation of estimates from single runs around the mean estimate is similar to that from the simulated weights method, as shown in the fourth and fifth rows, indicating our method will give a reasonable estimate of sampling errors when only a single nested sampling run is available.
The final two rows of Table 1 show the empirical coverage rates for bootstrap credible intervals are very close to their nominal values. Figure 6 shows estimates of the full sampling error distribution from a single run nested sampling run using the bootstrap and simulated weights methods; the bootstrap results are much closer to the sampling errors observed in repeated calculations, and give accurate estimates of the and credible intervals.
Section E of the supplementary material shows similar numerical tests for a 3-dimensional Cauchy likelihood (15). Even for this challenging, fat-tailed distribution our method performs similarly to the Gaussian case, giving accurate mean error estimates and estimates of credible intervals with measured coverage similar to their nominal coverage.
5.2 Sampling errors in different dimensions
We now verify the bootstrap method’s accuracy for Gaussian (14) and Cauchy (15) likelihoods of between 2 and 50 dimensions. Figure 7 shows bootstrap sampling error estimates accurately match the errors measured from repeated calculations, even for the challenging fat-tailed Cauchy distribution. In contrast the simulated weights method consistently underestimates the sampling errors in parameter estimation, although as expected it is accurate for errors on the evidence . See Section B of the supplementary material for a detailed discussion of the simulated weights method.
As the dimension increases, Figure 7 shows parameter estimation errors decreasing and the evidence errors increasing (with a constant number of live points ). This effect is due to the posterior being contained in a smaller fraction of the prior volume in higher dimensions. In the spherically symmetric cases considered, the range of to be explored increases approximately linearly with the dimension , as can be seen in Figure 4. With a constant number of live points, the number of samples is therefore also approximately proportional to .
In parameter estimation from posterior samples only points’ relative weights matter, so the increased number of samples in higher dimension problems typically increases accuracy as can be seen in Figures 7(a) and 7(b). However for high dimensional Cauchy likelihoods (15) the posterior mass is spread over a wide range of values, so errors in the relative weights of points become large in high dimensions666When the errors in points’ relative weights become dominant the simulated weights method captures the majority of the sampling error, as can be seen for high dimensional Cauchy distributions in Figure 7(b).
For the dominant error is in the absolute value of point weights, which is approximately proportional to the square root of the number of steps required to reach the posterior (Skilling, 2006). errors are therefore approximately proportional to when is constant, as can be seen in Figure 7(c).
6 Application to existing nested sampling software
Nested sampling software such as MultiNest and PolyChord can be easily modified to output information about the step at which dead points were sampled and give sampling error estimates using bootstrap resampling of threads. We are currently working on incorporating this into future releases of PolyChord.
Sampling error estimates from our approach will be accurate provided the software is performing nested sampling approximately correctly. However such software can only approximately sample randomly from the prior within iso-likelihood contours — this may result in additional errors which are specific to a given implementation and which may not be captured by general methods such as resampling threads. Tuning parameters such as “num_repeats” (the number of slice samples taken between dead points) in PolyChord allow reduced correlation between samples at a higher computational cost. Testing sampling error estimates from our method against those from repeated calculations can be used to detect implementation-specific errors and to select appropriate values for tuning parameters.
In addition our algorithm for dividing nested sampling runs could be used to detect implementation-specific errors by testing the difference in correlation between threads from the same run and from different runs. In principle this could allow estimates of sampling errors to be corrected for the effects of correlations between threads.
We now demonstrate our method’s application to nested sampling results produced with PolyChord.
6.1 Sampling errors on data fitting with PolyChord
We fit a set of points with normally distributed errors on the values using a sinusoid
[TABLE]
The likelihood is then
[TABLE]
where and we use a uniform prior for , and . Numerical tests use 40 data points sampled from with Gaussian noise of size added to the values; and the data points are shown in Figure 9 and the posterior distribution of given the data is shown in Figure 9. Posterior distributions on , and can be calculated with nested sampling — these are illustrated in Figure 10 along with example sampling errors.
Table 2 shows sampling errors from PolyChord with — the default value for a 3-dimensional problem. As in perfect nested sampling, our bootstrap estimates of the standard error agree with the variation in results observed, and the observed coverage of credible intervals is close to their nominal coverage. This implies that is sufficient for PolyChord to perform parameter estimation accurately in this case.
7 Conclusion
Sampling errors in nested sampling parameter estimation arise principally from two sources: uncertain sample weights , and approximating the average of a function of parameters on each iso-likelihood contour with a single sample . The latter error is not present in evidence calculation and has been previously ignored. The added stochasticity from sampling each iso-likelihood contour makes nested sampling parameter estimation a 2-dimensional problem, with a dependence on both the distribution of posterior mass and the distribution of parameter values on each iso-likelihood contour. We proposed a new diagram for representing both aspects of the calculation, and presented it in Figures 3(a), 3(c), 3(b), 3(d) and 3(e).
Estimating sampling errors is vital for interpreting the results of a nested sampling calculation, as well as for allocating computational resources — for example by choosing an appropriate number of live points. However the current approach (the simulated weights method) underestimates sampling errors as it does not account for approximating with a single sample . We proposed a new method for estimating sampling errors using our new algorithm (Algorithm 1) for dividing a nested sampling run into single live point runs (“threads”), which can then be resampled with techniques such as the bootstrap. This works as the values of the dead points from some nested sampling run with live points form a Poisson process with rate , meaning the values of the dead points of a single thread are a Poisson process of rate 1.
Our method shows accurate and robust estimation of sampling errors in parameter estimation in empirical tests, and compares favourably to the other methods discussed. The new method can be easily incorporated into existing nested sampling software, and will be reliable provided the implementation is performing the nested sampling algorithm accurately. We are currently working on including nested sampling run division and sampling error estimates from our method in future versions of PolyChord.
Supplementary material
Appendix A Relative contributions of different sources of parameter estimation sampling errors
The relative contributions of sampling errors from unknown prior weights of points and from taking a single sample on each iso-likelihood contour (discussed in Section 3) can be calculated by using exact values for weights and replacing with . For a Gaussian likelihood (14) and Gaussian prior (16) both and can be calculated analytically for each — sampling errors from calculations using this additional information are shown in Figure 11 and Table 3.
When calculating , as expected, using exact weights reduces uncertainty to the small trapezium rule error and using has no effect. However for parameter estimation significant error remains when using exact values777For the error increases when exact values are used. This is because the true weights are more variable than the expected ones and this reduces the information content (entropy) of the set of samples.. The relative contribution to sampling errors from estimating weights statistically is greatest when has a strong dependence on over the interval in containing the bulk of the posterior mass. In contrast when then for all , and the analysis using always gives the analytically correct answer of zero. In all cases, when both exact and samples from are used the sampling error is reduced to close to zero.
Appendix B Analysis of the simulated weights method
The simulation method underestimates sampling errors in nested sampling parameter estimation, as shown by the numerical tests in Tables 1 and 2 and Figure 7. This is because it assumes that for each dead point , neglecting the sampling errors from taking a single sample on each iso-likelihood contour which are described in Section 3. However some of this error is captured because repeatedly simulating points’ weights behaves like a resampling scheme, with similarities to the Bayesian bootstrap (Rubin, 1981). Resampling estimates the uncertainty on inferences from a set of samples by calculating its variation when data points are re-weighted, but the simulated weights method does so in a way that systematically underestimates sampling errors. This behavior has not been documented in the literature.
For example, consider the case with a Gaussian likelihood (14) and Gaussian prior (16) — here for all . If is used instead of there is no sampling error on estimates of regardless of any uncertainty in the weights of each point , as can be seen in Figure 11 and Table 3. However the simulated weights method gives a non-zero estimate which on average differs from the sampling errors measured by repeated calculations by a factor of very close to , as shown in the third row of Table 1.
Further numerical tests show that in special cases when is the same at all the ratio of sampling errors from the simulated weights method to the error observed in repeated calculations has a value close to . We give an analytical explanation for this result below. However we note that for practical problems is a priori unknown and likely varies in , meaning the true sampling error cannot be predicted by adjusting estimates from the simulated weights method.
B.1 Sampling error estimates for special cases when is constant for all
B.1.1 Variance of sampling error distribution
Nested sampling calculates the expected value of a function of parameters as . Here the sampling error is the difference between the exact value of from the posterior, and is distributed as
[TABLE]
The variance of this distribution provides a measure of the size of the sampling error. As the nested sampling estimator is unbiased, the variance of the sampling error distribution is equal to the variance of the results of repeated calculations:
[TABLE]
Expanding and dropping the explicit dependence of and on and for brevity gives
[TABLE]
where .
Each is an independent random variable from the distribution , so the expectation of products of are zero for all . Furthermore expectation of products and the covariance are zero for .
The weights have a dependence on , but in the case is the same for all the covariance terms are also zero for all . (32) therefore simplifies to
[TABLE]
Expanding the variance term on the right hand side when is constant and and are therefore independent gives
[TABLE]
B.1.2 Simulated weights method variance estimate
The simulated weights method corresponds to fixing the values while retaining the stochastic dependence of on . This means taking , , which combined with (34) gives
[TABLE]
Taking the expected values for and this becomes
[TABLE]
Using that by definition so ,
[TABLE]
In contrast the repeated runs method retains the sampling error of on and uses the expected values of the weight . Hence for a large number of trials , for all . Subbing into (34) gives
[TABLE]
Using that ,
[TABLE]
B.1.3 Ratio of simulated weights and repeated runs variance estimates
Combining equations (37) and (39) gives the ratio of the simulated weights method and repeated runs variances as
[TABLE]
If is the same for all this simplifies to
[TABLE]
By definition the normalised weights so
[TABLE]
Numerical results suggest that for a range of problems and are approximately independent, in which case
[TABLE]
Typical problems with a large often also have , in which case
[TABLE]
Keeton (2011) gives expressions for the weights as888These formulae omit the trapezium rule and for brevity take — this approximation has little effect on the results.
[TABLE]
For a general likelihood the summation in (45) cannot be found exactly. However one can estimate the ratio for each live point
[TABLE]
This supports the observation that the ratio of simulated weights method estimates of the standard deviation of stochastic errors to measurements from repeated runs is close to for special cases such as calculating the mean of a parameter for spherically symmetric likelihoods with spherically symmetric co-centred priors.
B.1.4 Ratio in the special case where and are constant for all
If the likelihood is constant999We assume has an infinitesimal slope to give direction to nested sampling’s inward iteration. throughout the parameter space and is the same for all then the likelihood terms in the numerator and denominator of (45) cancel and the summation can be found exactly. Furthermore estimates of are very precise in this case as there is no stochastic variation in , justifying the approximation (45). In this case the ratio is
[TABLE]
where the final step sums the geometric series and neglects terms from the truncation of the sum due to termination of the nested sampling run.
Appendix C Split runs method
Instead of spending all available computational resources on a single nested sampling run with live points, one might consider performing smaller runs with live points and estimating the sampling error from the variation of the smaller runs — for example as times their sample standard deviations. However this provides a limited number of sub-runs, and does not give accurate credible interval estimates. Furthermore while sampling errors in nested sampling are typically proportional to , this breaks down when the number of samples is small due to trapezium rule errors in sample weights which are . As a result multiple runs are best analysed by combining them into a single run (Skilling, 2006).
Sampling error estimates from taking the standard deviation of the results of sub-runs and multiplying by are shown in Table 4. The split runs method is inaccurate for the approximately log-normally distributed sampling errors in as well as for credible intervals on distribution tails such as , as can be seen in the third row of Table 4.
Appendix D Termination conditions
The sensitivity to termination conditions can be far higher for parameter estimation than for evidence calculation. This is both because parameter estimation can have much smaller sampling errors, and because the region close to the likelihood peak can have very high weight for some . For example for the Gaussian likelihood (14) an estimator such as may show significant errors due to termination conditions which were perfectly adequate for calculating . Numerical tests in this paper use the termination conditions described by Handley et al. (2015b, Section 3.4), stopping when the estimated evidence contained in the live points is less than times the evidence contained in dead points.
When splitting runs into their constituent threads (Section 4.1) then even in perfect nested sampling termination conditions must be chosen carefully to avoid causing differences between threads from different runs which terminate at different likelihoods. This typically happens when
termination conditions are worked out from the current set of dead points — e.g. estimating the evidence remaining as in Handley et al. (2015b, Section 3.4). This means some runs continue for longer than others; 2. 2.
the final point which violates the condition is kept. This means threads from small runs are much more likely to have final points far exceeding the termination condition than threads from large runs.
When comparing threads from different nested sampling runs, their equivalence can be maintained by using a termination condition which does not infer anything from the previous points, such as setting a fixed likelihood value for termination and discarding any point that exceeds it. As we do not mix threads from different runs in our numerical tests we do not need this approach.
Appendix E Additional numerical tests: 3-dimensional Cauchy likelihood
Table 5 shows numerical tests of sampling error estimates with a 3-dimensional Cauchy likelihood (15) with a Gaussian prior (16). As in the 3-dimensional Gaussian case shown in Table 1, the mean estimates of sampling errors from our bootstrap method are very close to measurements of sampling errors from repeated calculations — this can be seen in the second row of Table 5. Again the empirical coverage rates for bootstrap credible intervals are close to their nominal values, as shown in the final two rows.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Bezáková et al. (2008) Bezáková, I., Stefankovic, D., Vazirani, V. V., and Vigoda, E. (2008). “Accelerating simulated annealing for the permanent and combinatorial counting problems.” SIAM , 37(5): 1429–1454.
- 2Chopin and Robert (2010) Chopin, N. and Robert, C. P. (2010). “Properties of nested sampling.” Biometrika , 97(3): 741–755.
- 3Del Moral et al. (2006) Del Moral, P., Doucet, A., and Jasra, A. (2006). “Sequential Monte Carlo samplers.” Journal of the Royal Statistical Society. Series B: Statistical Methodology , 68(3): 411–436.
- 4Doss et al. (2015) Doss, C., Flegal, J., Jones, G., and Neath, R. (2015). “Markov chain monte carlo estimation of quantiles.” Electronic Journal of Statistics , 8: 2448–2478.
- 5Efron (1979) Efron, B. (1979). “Bootstrap Methods: Another Look at the Jackknife.” The Annals of Statistics , 7(1): 1–26.
- 6Efron and Tibshirani (1986) Efron, B. and Tibshirani, R. (1986). “Bootstrap Methods for Standard Errors, Confidence Intervals, and Other Measures of Statistical Accuracy.” Statistical Science , 1(1): 54–77.
- 7Feroz and Hobson (2008) Feroz, F. and Hobson, M. P. (2008). “Multimodal nested sampling: An efficient and robust alternative to Markov Chain Monte Carlo methods for astronomical data analyses.” Monthly Notices of the Royal Astronomical Society , 384(2): 449–463.
- 8Feroz et al. (2009) Feroz, F., Hobson, M. P., and Bridges, M. (2009). “Multi Nest: An efficient and robust Bayesian inference tool for cosmology and particle physics.” Monthly Notices of the Royal Astronomical Society , 398(4): 1601–1614.
