We present DYNESTY a public opensource PYTHON package to estimate Bayesian posteriors and evidences marginal likelihoods using the dynamic nested sampling methods developed by Higson et al. By adaptively allocating samples based on posterior structure dynamic nested sampling has the benefits of Markov chain Monte Carlo MCMC algorithms that focus exclusively on posterior estimation while retaining nested samplings ability to estimate evidences and sample from complex multimodal distributions. We provide an overview of nested sampling its extension to dynamic nested sampling the algorithmic challenges involved and the various approaches taken to solve them in this and previous work. We then examine DYNESTYs performance on a variety of toy problems along with several astronomical applications. We find in particular problems DYNESTY can provide substantial improvements in sampling efficiency compared to popular MCMC approaches in the astronomical literature. More detailed statistical results related to nested sampling are also included in the appendix.
['methods: data analysis', 'methods: statistical', 'Astrophysics - Instrumentation and Methods for Astrophysics', 'Statistics - Computation']
DYNESTY a dynamic nested sampling package for estimating Bayesian posteriors and evidences
{'dynesty : A Dynamic Nested Sampling Package for Estimating Bayesian Posteriors and Evidences': 'Joshua S. Speagle 1 ; 2 ?  \nCenter for Astrophysics | Harvard & Smithsonian, 60 Garden St., Cambridge, MA, USA  \n1 2 NSF Graduate Research Fellow  \nAccepted XXX. Received YYY; in original form ZZZ', 'ABSTRACT': "We present dynesty , a public, open-source, Python package to estimate Bayesian posteriors and evidences (marginal likelihoods) using Dynamic Nested Sampling. By adaptively allocating samples based on posterior structure, Dynamic Nested Sampling has the benefits of Markov Chain Monte Carlo algorithms that focus exclusively on posterior estimation while retaining Nested Sampling's ability to estimate evidences and sample from complex, multi-modal distributions. We provide an overview of Nested Sampling, its extension to Dynamic Nested Sampling, the algorithmic challenges involved, and the various approaches taken to solve them. We then examine dynesty 's performance on a variety of toy problems along with several astronomical applications. We find in particular problems dynesty can provide substantial improvements in sampling efficiency compared to popular MCMC approaches in the astronomical literature. More detailed statistical results related to Nested Sampling are also included in the Appendix.  \nKey words: methods: statistical - methods: data analysis", '1 INTRODUCTION': 'to give  \nMuch of modern astronomy rests on making inferences about underlying physical models from observational data. Since the advent of large-scale, all-sky surveys such as SDSS (York et al. 2000), the quality and quantity of these data increased substantially (Borne et al. 2009). In parallel, the amount of computational power to process these data also increased enormously. These changes opened up an entire new avenue for astronomers to try and learn about the universe using more complex models to answer increasingly sophisticated questions over large datasets. As a result, the standard statistical inference frameworks used in astronomy have generally shifted away from Frequentist methods such as maximum-likelihood estimation (MLE; Fisher 1922) to Bayesian approaches to estimate the distribution of possible parameters for a given model that are consistent with the data and our current astrophysical knowledge (see, e.g., Trotta 2008; Planck Collaboration et al. 2016; Feigelson 2017).  \nIn the context of Bayesian inference, we are interested in estimating the posterior P \' GLYPH<2> j D ; M \' of a set of parameters GLYPH<2> for a given model M conditioned on some data D . This can be written into a form commonly known as Bayes Rule  \n- ? E-mail: jspeagle@cfa.harvard.edu  \nP \' GLYPH<2> j D ; M \' = P \' D j GLYPH<2> ; M \' P \' GLYPH<2> j M \' P \' D j M \' (1)  \nwhere P \' D j GLYPH<2> ; M \' is the likelihood of the data given the parameters of our model, P \' GLYPH<2> j M \' is the prior for the parameters of our model, and  \nP \' D j M \' = \' GLYPH<10> GLYPH<2> P \' D j GLYPH<2> ; M \' P \' GLYPH<2> j M \' d GLYPH<2> (2)  \nis the evidence (i.e. marginal likelihood) for the data given our model, where the integral is taken over the entire domain GLYPH<10> GLYPH<2> of GLYPH<2> (i.e. over all possible parameter combinations). Throughout the rest of the paper, we will refer to these using shorthand notation  \nP\' GLYPH<2> M \' = L\' GLYPH<2> M \' GLYPH<25> \' GLYPH<2> M \' Z M (3)  \nwhere P\' GLYPH<2> M \' GLYPH<17> P \' GLYPH<2> j D ; M \' is the posterior, L\' GLYPH<2> M \' GLYPH<17> P \' D j GLYPH<2> ; M \' is the likelihood, GLYPH<25> \' GLYPH<2> M \' GLYPH<17> P \' GLYPH<2> j M \' is the prior, Z M GLYPH<17> P \' D j M \' is the evidence, and the subscript M will subsequently be dropped if we are only considering a single model. Here, the posterior P\' GLYPH<2> M \' tells us about the parameter estimates from a given model M while Z M enables us to compare across models marginalized over any particular set of parameters using the Bayes factor:  \nR GLYPH<17> Z M 1 Z M 2 GLYPH<25> \' M 1 \' GLYPH<25> \' M 2 \' (4)  \nFigure 1. A schematic representation of the different approaches Markov Chain Monte Carlo (MCMC) methods and Nested Sampling methods take to sample from the posterior. While MCMC methods attempt to generate samples directly from the posterior, Nested Sampling instead breaks up the posterior into many nested \'slices\', generates samples from each of them, and then recombines the samples to reconstruct the original distribution using the appropriate weights.  \n<!-- image -->  \nwhere GLYPH<25> \' M i \' is the prior belief in model M i .  \nFor complicated data and models, the posterior P\' GLYPH<2> \' is often analytically intractable and must be estimated using numerical methods. These fall into two broad classes: \'approximate\' and \'exact\' approaches. Approximate approaches try to find an (analytic) distribution Q\' GLYPH<2> \' that is \'close\' to P\' GLYPH<2> \' using techniques such as Variational Inference (Blei et al. 2016). These techniques are not the focus of this work and will not be discussed further in this paper.  \nExact approaches try to estimate P\' GLYPH<2> \' directly, often by constructing an algorithm that allows us to generate a set of samples f GLYPH<2> 1 ; GLYPH<2> 2 ; : : : ; GLYPH<2> N g that we can use to approximate the posterior as a weighted collection of discrete points  \nP\' GLYPH<2> \' GLYPH<25> ˆ P\' GLYPH<2> \' = " N i = 1 p \' GLYPH<2> i \' GLYPH<14> \' GLYPH<2> i \' " N i = 1 p \' GLYPH<2> i \' (5)  \nwhere p \' GLYPH<2> i \' is the importance weight associated with each GLYPH<2> i and GLYPH<14> \' GLYPH<2> i \' is the Dirac delta function located at GLYPH<2> i .  \nThere is a rich literature (see, e.g., Chopin & Ridgway 2015) on the approaches used to generate these samples and their associated weights. The most popular method used in astronomy today is Markov Chain Monte Carlo (MCMC), which generates samples\'proportional to\'the posterior such that p i = 1 . While MCMC has had substantial success over the past few decades (Brooks et al. 2011; Sharma 2017), the most common implementations (e.g., Plummer 2003; Foreman-Mackey et al. 2013; Carpenter et al. 2017) tend to struggle when the posterior is comprised of widely-separated modes. In addition, because it only generates samples proportional to the posterior, it is difficult to use those samples to estimate the evidence Z M to compare various models.  \nNested Sampling (Skilling 2004; Skilling 2006) is an al-  \nternative approach to posterior and evidence estimation that tries to resolve some of these issues. 1 By generating samples in nested (possibly disjoint) \'shells\' of increasing likelihood, it is able to estimate the evidence Z M for distributions that are challenging for many MCMC methods to sample from. The final set of samples can also be combined with their associated importance weights p i to generate associated estimates of the posterior. 2  \nSince a large portion of modern astronomy relies on being able to perform Bayesian inference, implementing these methods often can serve as the primary bottleneck for testing hypotheses, estimating parameters, and performing model comparisons. As such, packages that implement these approaches serve an important role enabling science by bridging the gap between writing down a model and estimating its associated parameters. These allow users to perform sophisticated analyses without having to implement many of the aforementioned algorithms themselves. Several prominent examples include the MCMC package emcee (Foreman-Mackey et al. 2013) and the Nested Sampling packages MultiNest (Feroz et al. 2009, 2013) and PolyChord (Handley et al. 2015), which collectively have been used in thousands of papers.  \nWe present dynesty , a public, open-source, Python package that implements Dynamic Nested Sampling.  \ndynesty is designed to be easy to use and highly modular, with extensive documentation, a straightforward application programming interface (API), and a variety of sampling implementations. It also contains a number of \'quality of life\' features including well-motivated stopping criteria, plotting functions, and analysis utilities for post-processing results.  \nThe outline of the paper is as follows. In § 2 we give an overview of Nested Sampling and discuss the method\'s benefits and drawbacks. In § 3 we describe how Dynamic Nested Sampling is able to resolve some of these drawbacks by allocating samples more flexibly. In § 4 we discuss the specific approaches dynesty uses to track and sample from complex, multi-modal distributions. In § 5 we examine dynesty \'s performance on a variety of toy problems. In § 6 we examine dynesty \'s performance on several real-world astrophysical analyses. We conclude in § 7. For interested readers, more detailed results on many of the methods outlined in the main text are included in Appendix A.  \ndynesty is publicly available on GitHub as well as on PyPI. See https://dynesty.readthedocs.io for installation instructions and examples on getting started.', '2 NESTED SAMPLING': "The general motivation for Nested Sampling, first proposed by Skilling (2004) and later fleshed out in Skilling (2006), stems from the fact that sampling from the posterior P' GLYPH<2> ' directly is hard . Methods such as Markov Chain Monte Carlo (MCMC) attempt to tackle this single difficult problem directly . Nested Sampling, however, instead tries to break down this single hard problem into a larger number of simpler problems by:  \n- (i) 'slicing' the posterior into many simpler distributions,\n- (ii) sampling from each of those in turn, and\n- (iii) re-combining the results afterwards.  \nWe provide a schematic illustration of this procedure in Figure 1 and give a broad overview of this process below. For additional details, please see Appendix A.", '2.1 Overview': "Unlike MCMC methods, which attempt to estimate the posterior P' GLYPH<2> ' directly, Nested Sampling instead focuses on estimating the evidence  \nZ GLYPH<17> ' GLYPH<10> GLYPH<2> P' GLYPH<2> ' d GLYPH<2> = ' GLYPH<10> GLYPH<2> L' GLYPH<2> ' GLYPH<25> ' GLYPH<2> ' d GLYPH<2> (6)  \nAs this integral is over the entire multi-dimensional domain of GLYPH<2> , it is traditionally very challenging to estimate.  \nNested Sampling approaches this problem by refactoring this integral as one taken over prior volume X of the enclosed parameter space  \nZ = ' GLYPH<10> GLYPH<2> L' GLYPH<2> ' GLYPH<25> ' GLYPH<2> ' d GLYPH<2> = ' 1 0 L' X ' d X (7)  \nHere, L' X ' now defines an iso-likelihood contour (or multiple) defining the edge(s) of the volume X , while the prior volume  \nX ' GLYPH<21> ' GLYPH<17> ' GLYPH<10> GLYPH<2> : L' GLYPH<2> 'GLYPH<21> GLYPH<21> GLYPH<25> ' GLYPH<2> ' d GLYPH<2> (8)  \nMNRAS 000 , 1-28 (2019)  \nis the fraction of the prior where the likelihood L' GLYPH<2> ' GLYPH<21> GLYPH<21> is above some threshold GLYPH<21> . Since the prior is normalized, this gives X ' GLYPH<21> = 0 ' = 1 and X ' GLYPH<21> = 1' = 0 , which define the bounds of integration for equation (7).  \nAs a rough analogy, we can consider trying to integrate over a spherically-symmetric distribution in 3-D. While it is possible to integrate over d x d y d z directly, it often is significantly easier to instead integrate over differential volume elements d V = 4 GLYPH<25> r 2 as a function of radius r GLYPH<17> p x 2 + y 2 + z 2 :  \n' P' x ; y ; z ' d x d y d z = ' P' V ' r '' d V ' r ' = ' P' r ' 4 GLYPH<25> r 2 d r  \nParameterizing the evidence integral this way allows Nested Sampling (in theory) to convert from a complicated D -dimensional integral over GLYPH<2> to a simple 1-D integral over X .  \nWhile it is straightforward to evaluate the likelihood at a given position L' GLYPH<2> ' , estimating the associated prior volume X ' GLYPH<2> ' and its differential dX ' GLYPH<2> ' is substantially more challenging. We can, however, generate noisy estimates of these quantities by employing the procedure described in Algorithm 1. We elaborate further on this procedure and how it works below.", '2.2 Generating Samples': "A core element of Nested Sampling is the ability to generate samples from the prior GLYPH<25> ' GLYPH<2> ' subject to a hard likelihood constraint GLYPH<21> . The most naive algorithm that satisfies this constraint is simple rejection sampling: at a given iteration i , generate samples GLYPH<2> i + 1 from the prior GLYPH<25> ' GLYPH<2> ' until L' GLYPH<2> i + 1 ' GLYPH<21> L' GLYPH<2> i ' .  \nIn practice, however, this simple procedure becomes progressively less efficient as time goes on since the remaining prior volume X i + 1 at each iteration of Algorithm 1 keeps shrinking. We therefore need a way of directly generating samples from the constrained prior:  \nGLYPH<25> GLYPH<21> ' GLYPH<2> ' GLYPH<17> ( GLYPH<25> ' GLYPH<2> 'GLYPH<157> X ' GLYPH<21> ' L' GLYPH<2> ' GLYPH<21> GLYPH<21> 0 L' GLYPH<2> ' < GLYPH<21> (9)  \nSampling from this constrained distribution is difficult for an arbitrary prior GLYPH<25> ' GLYPH<2> ' since the density can vary drastically from place to place. It is simpler, however, if the prior is standard uniform (i.e. flat from 0 to 1 ) in all dimensions so that the density interior to GLYPH<21> is constant then X behaves more like a typical volume V . We can accomplish this through the use of the appropriate 'prior transform' function T which maps a set of parameters GLYPH<8> with a uniform prior over the D -dimensional unit cube to the parameters of interest GLYPH<2> . 3 Taken together, these transform our original hard problem of sampling from the posterior P' GLYPH<2> ' directly to instead the much simpler problem of repeatedly sampling uniformly 4 within the transformed constrained prior  \nGLYPH<25> 0 GLYPH<21> ' GLYPH<8> ' GLYPH<17> ( 1 GLYPH<157> X ' GLYPH<21> ' L' GLYPH<2> = T' GLYPH<8> '' GLYPH<21> GLYPH<21> 0 otherwise (10)  \n3 In general, there is a uniquely defined prior transform T for any given GLYPH<25> ' GLYPH<2> ' ; see the dynesty documentation for additional details. 4 Technically this requirement is overly strict, as Nested Sampling can still be valid even if the samples at each iteration are correlated. See Appendix A for additional discussion.", 'Algorithm 1: Static Nested Sampling': "// Initialize live points.  \nDraw K 'live' points f GLYPH<2> 1 ; : : : ; GLYPH<2> K g from the prior GLYPH<25> ' GLYPH<2> ' .  \n// Main sampling loop.  \nwhile stopping criterion not met do  \nCompute the minimum likelihood L min among the current set of live points. Add the k th live point GLYPH<2> k associated with L min to a list of 'dead' points. Sample a new point GLYPH<2> 0 from the prior subject to the constraint L' GLYPH<2> 0 ' GLYPH<21> L min . Replace GLYPH<2> k with GLYPH<2> 0 .  \n// Check whether to stop.  \nEvaluate stopping criterion.", 'end': "// Add final live points.  \nwhile K > 0 do  \nCompute the minimum likelihood L min among the current set of live points. Add the k th live point GLYPH<2> k associated with L min to a list of 'dead' points. Remove GLYPH<2> k from the set of live points. Set K = K GLYPH<0> 1 .  \nend  \nThroughout the rest of the text we will henceforth assume GLYPH<25> ' GLYPH<2> ' is a unit cube prior unless otherwise explicitly specified.  \nBecause there is no constraint that this distribution is uni-modal, the constrained prior may define several 'blobs' of prior volume that we are interested in sampling from. While sampling from the blob(s) might be hard to do from scratch, because Nested Sampling samples at many different likelihood 'levels', structure tends to emerge over time rather than all at once as we transition away from the prior GLYPH<25> ' GLYPH<2> ' .", '2.3 Estimating the Prior Volume': "As shown in Appendix A, generating samples following the strategy in § 2.2 based on Algorithm 1 allows us to estimate the (change in) prior volume at a given iteration using the set of 'dead' points (i.e. the live points we replaced at each iteration). In particular, it leads to exponential shrinkage such that the (log-)prior volume at each iteration changes by  \nE GLYPH<2> GLYPH<1> ln ˆ X i GLYPH<3> = E GLYPH<2> ln ˆ X i GLYPH<0> ln ˆ X i GLYPH<0> 1 GLYPH<3> = GLYPH<0> 1 K (11)  \nwhere E »GLYPH<1>… is the expectation value (i.e. mean) and we have adopted the ˆ x notation to emphasize that we have a noisy estimator of the prior volume X . Using more live points K thus increases our volume resolution by decreasing the rate of this exponential compression. By default, dynesty uses K = 500 live points, although this should be adjusted depending on the problem at hand.  \nOnce some stopping criterion is reached and sampling terminates after N iterations, the remaining set of K live points are then distributed uniformly within the final prior volume X N (see Appendix A). These can be 'recycled' into the final set of samples by sequentially adding the live points to the list of'dead'points collected at each iteration in order of increasing likelihood. This leads to uniform shrinkage of the prior volume such that the (fractional) change in prior volume for the k th live point added this way is  \nE GLYPH<20> GLYPH<1> ˆ X N + k ˆ X N GLYPH<21> = E GLYPH<20> ˆ X N + k GLYPH<0> ˆ X N + k GLYPH<0> 1 ˆ X N GLYPH<21> = 1 K + 1 (12)  \nwhere ˆ X N is the estimating remaining prior volume at the final N th iteration.", '2.4 Stopping Criterion': "Since Nested Sampling is designed to estimate the evidence, a natural stopping criterion (see, e.g., Skilling 2006; Keeton 2011; Higson et al. 2017a) is to terminate sampling when we believe our set of dead points (and optionally the remaining live points) give us an integral that encompasses the vast majority of the posterior. In other words, at a given iteration i , we want to terminate sampling if  \nGLYPH<1> ln ˆ Z i GLYPH<17> ln GLYPH<16> ˆ Z i + GLYPH<1> ˆ Z i GLYPH<17> GLYPH<0> ln GLYPH<16> ˆ Z i GLYPH<17> < GLYPH<15> (13)  \nwhere GLYPH<1> ˆ Z i is the estimated remaining evidence we have yet to integrate over and GLYPH<15> determines the tolerance. If the final set of live points are excluded from the set of dead points, dynesty assumes a default value of GLYPH<15> = 10 GLYPH<0> 2 (i.e. glyph[lessorsimilar] 1% of the evidence remaining). If the final set of live points are included, dynesty instead uses the slightly more permissive GLYPH<15> = 10 GLYPH<0> 3 ' K GLYPH<0> 1 ' + 10 GLYPH<0> 2 .  \nWhile the remaining evidence GLYPH<1> ˆ Z i is unknown, we can in theory construct a strict upper bound on it by assigning  \nGLYPH<1> ˆ Z i GLYPH<20> L max X i (14)  \nwhere L max is the maximum-likelihood value across the entire domain GLYPH<10> GLYPH<2> and X i is the prior volume at the current iteration. This is equivalent to treating the remaining likelihood interior to the current sample ( X < X i ) as a uniform slab with amplitude L max .  \nUnfortunately, neither L max or X i is known exactly. However, we can approximate this upper bound by replacing both quantities with associated estimators to get the rough upper bound  \nGLYPH<1> ˆ Z i glyph[lessorsimilar] L max i ˆ X i (15)  \nwhere L max i is the maximum value of the likelihood among the live points at iteration i and ˆ X i is the estimated (remaining) prior volume.  \nWhile this rough upper bound works well in most cases,  \nFigure 2. An example highlighting the behavior of a Static Nested Sampling run in dynesty . See § 2 for additional details. Top: The number of live points as a function of prior volume X . Snapshots of their their distribution (purple) with respect to the current bounds (gray; see § 4.1) are highlighted in several insets. The number of live points remains constant until sampling terminates, at which point we add the final live points one-by-one to the samples. Top-middle: The (normalized) likelihood limit LGLYPH<157> L max associated with a the prior volume X 'L' in the top panel. This increases monotonically as we sample increasingly smaller regions of the prior. Bottom-middle: The importance weight PDF p ' X ' , roughly divided into regions dominated by the prior volume ( d X is large, L' X ' is small; yellow), posterior mass ( d X and L' X ' are comparable; orange), and likelihood density ( d X is small, L' X ' is large; red). The posterior mass is the most important for posterior estimation, while evidence estimation also depends on the prior volume. Bottom: The estimated evidence ˆ Z' X ' (blue line) and its 1, 2, and 3-sigma errors (blue shaded). The true value is shown in red.  \n<!-- image -->  \nbecause we only have access to the best likelihood L max i sampled by the K live points at a particular iteration there is always a chance that L max i GLYPH<28> L max and that we will terminate early. This can happen if there is an extremely narrow likelihood peak within the remaining prior volume that has not yet been discovered by the K live points.", '2.5 Estimating the Evidence and Posterior': 'Once we have a final set of samples f GLYPH<2> 1 ; : : : ; GLYPH<2> N g , we can estimate the 1-D evidence integral using standard numerical techniques. To ensure approximation errors on the numerical integration estimate are sufficiently small, dynesty uses the 2nd-order trapezoid rule  \nˆ Z = N + K GLYPH<213> i = 1 1 2 »L\' GLYPH<2> i GLYPH<0> 1 \' + L\' GLYPH<2> i \'… GLYPH<2> GLYPH<2> ˆ X i GLYPH<0> 1 GLYPH<0> ˆ X i GLYPH<3> GLYPH<17> N + K GLYPH<213> i = 1 ˆ p i (16)  \nwhere X 0 = X \' GLYPH<21> = 0 \' = 1 and  \nˆ p i GLYPH<17> »L\' GLYPH<2> i GLYPH<0> 1 \' + L\' GLYPH<2> i \'… GLYPH<2> GLYPH<2> ˆ X i GLYPH<0> 1 GLYPH<0> ˆ X i GLYPH<3> (17)  \nis the estimated importance weight. By default, dynesty uses the mean values of ˆ X i to compute the mean and standard deviation of ln ˆ Z following Appendix A, although these values can also be simulated explicitly.  \nWe can also estimate the posterior P\' GLYPH<2> \' from the same set of N + K dead points by using the associated importance weights derived above:  \nˆ P\' GLYPH<2> \' = " N + K i = 1 ˆ p \' GLYPH<2> i \' GLYPH<14> \' GLYPH<2> i \' " N + K i = 1 ˆ p \' GLYPH<2> i \' = ˆ Z GLYPH<0> 1 N + K GLYPH<213> i = 1 ˆ p \' GLYPH<2> i \' GLYPH<14> \' GLYPH<2> i \' (18)  \nBy default, dynesty uses the mean values of ˆ X i to compute this posterior estimate, although as with the evidence these values can also be simulated explicitly (see Appendix A).  \nAn illustration of a typical Nested Sampling run is shown in Figure 2.', '2.6 Benefits of Nested Sampling': "Because of its alternative approach to sampling from the posterior, Nested Sampling has a number of benefits relative to traditional MCMC approaches:  \n- (i) Nested Sampling can estimate the evidence Z as well as the posterior P' GLYPH<2> ' . MCMC methods generally can only constrain the latter (although see Lartillot & Philippe 2006; Heavens et al. 2017).\n- (ii) Nested sampling can sample from multi-modal distributions that tend to challenge many MCMC methods.\n- (iii) While most MCMC stopping criteria based on effective sample sizes can feel arbitrary, Nested Sampling possesses well-motivated stopping criteria focused on evidence estimation.  \n(iv) MCMC methods need to converge (i.e. 'burn in') to the posterior before any samples generated are valid. While optimization techniques can speed up this process, assessing this convergence can be challenging and time-consuming (Gelman & Rubin 1992; Vehtari et al. 2019). Nested Sampling doesn't suffer from similar issues because the method smoothly integrates over the posterior P' GLYPH<2> ' starting from the prior GLYPH<25> ' GLYPH<2> ' .", '2.7 Drawbacks': "While Nested Sampling has its fair share of benefits that have encouraged its rapid adoption in astronomical Bayesian analyses, it also suffers from a fair share of drawbacks. Most crucially, the standard Nested Sampling implementation outlined in Algorithm 1 focuses exclusively on estimating the evidence Z ; the posterior P' GLYPH<2> ' is entirely a by-product of the approach. This creates several immediate drawbacks relative to MCMC, which focuses exclusively on sampling the posterior P' GLYPH<2> ' .  \nFirst, because most Nested Sampling implementations rely on sampling from uniform distributions (see § 2.2), applying them to general distributions requires knowing the appropriate prior transform T . While these are straightforward to define when the prior can be decomposed into separable, independent components, they can be more difficult to derive when the prior involves conditional and/or jointly distributed parameters.  \nSecond, because the evidence depends on the amount of prior volume that needs to be integrated over, the overall expected runtime is sensitive to the relative size of the prior. In other words, while estimating the posterior mostly depends on generating samples close to where the majority of the distribution is located (i.e. the 'typical set'; Betancourt 2017), estimating the evidence requires generating samples in the extended tails of the distribution. Using less informative (broader) priors will increase the expected runtime even if the posterior is largely unchanged.  \nFinally, because the number of live points K is constant, the rate GLYPH<1> ln X at which we integrate over the posterior P' GLYPH<2> ' is the same regardless of where we are. This means that increasing the number of like points K , which increases the overall runtime, always improves the accuracy of both the posterior ˆ P' GLYPH<2> ' and evidence ˆ Z estimates. In other words, Nested Sampling does not allow users to prioritize between estimating the posterior or the evidence, which is not ideal for many analyses that are mostly interested in using Nested Sampling for either option. We focus on improving this behavior in § 3.  \nAs with any sampling method, we strongly advocate that Nested Sampling should not be viewed as being strictly 'better' or 'worse' than MCMC, but rather as a tool that can be more or less useful in certain problems. There is no 'One True Method to Rule Them All', even though it can be tempting to look for one.", '3 DYNAMIC NESTED SAMPLING': "In our overview of Nested Sampling in § 2, we highlighted three main drawbacks of basic implementations:  \n- (i) They generally require a prior transform.\n- (ii) Their runtime is sensitive to the size of the prior.\n- (iii) Their rate of posterior integration is always constant.  \nWhile the first two drawbacks are essentially inherent to Nested Sampling as sampling strategy, the last is not. Instead, the inability of Algorithm 1 to 'prioritize' estimating the evidence Z or posterior P' GLYPH<2> ' is a consequence of the fact that the number of live points K remains constant throughout an entire run, which sets the rate of integration GLYPH<1> ln X .", 'Algorithm 2: Dynamic Nested Sampling': "```\n// Initialize first set of live points. Draw K 'live' points f GLYPH<2> 1 ; : : : ; GLYPH<2> K g from the prior GLYPH<25> ' GLYPH<2> ' . // Main sampling loop. Set L min = 0 and K 0 = K . while stopping criterion not met do // Get current number of live points. Compute the previous number of live points K and the current number of live points K 0 . if K 0 GLYPH<21> K then // Add in new live points. while K 0 > K do Sample a new point GLYPH<2> 0 from the prior subject to the constraint L' GLYPH<2> 0 ' GLYPH<21> L min . Add GLYPH<2> 0 to the set of live points. Set K = K + 1 . end // Replace worst live point. Compute the minimum likelihood L min among the current set of K live points. Add the k th live point GLYPH<2> k associated with L min to a list of 'dead' points. Replace GLYPH<2> k with GLYPH<2> 0 . else // Iteratively remove live points. while K 0 < K do Compute the minimum likelihood L min among the current set of K = K 0 live points. Add the k th live point GLYPH<2> k associated with L min to a list of 'dead' points. Remove GLYPH<2> k from the set of live points. Set K = K GLYPH<0> 1 . end end // Check whether to stop. Evaluate stopping criterion. end // Add final live points. while there are live points remaining do Compute the minimum likelihood L min among the current set of live points. Add the k th live point GLYPH<2> k associated with L min to a list of 'dead' points. Remove GLYPH<2> k from the set of live points. end\n```  \nAs a result, we will henceforth call this procedure 'Static' Nested Sampling.  \nTo address this issue, Higson et al. (2017b) proposed a deceptively simple modification: let the number of live points vary during runtime. This gives a new 'Dynamic' Nested Sampling algorithm whose basic implementation is outlined in Algorithm 2. This simple change is transformative, allowing Dynamic Nested Sampling to focus on sampling the posterior P' GLYPH<2> ' , similar to MCMC approaches, while retaining all the benefits of (Static) Nested Sampling to estimate the evidence Z and sample from complex, multi-modal distributions. It also possesses well-motivated new stopping criteria for posterior and evidence estimation.  \nIt is important to note that we cannot take advantage of the flexibility offered by Dynamic Nested Sampling, however, without implementing appropriate schemes to specify exactly how live points should be allocated, when to terminate sampling, etc. While dynesty tries to implement a number of reasonable default choices, in practice this inevitably leads to many more tuning parameters that can affect the behavior of a given Dynamic Nested Sampling run.  \nWe provide an illustration of the overall approach in  \nFigure 3 and give a broad overview of the basic algorithm below. For additional details, please see Appendix A.", '3.1 Allocating Live Points': "The singular defining feature of the Dynamic Nested Sampling algorithm is the scheme we use for determining how the number of live points K i at a given iteration i should vary. Naively, we would like K i to be larger where we want our resolution to be higher (i.e. a slower rate of integration GLYPH<1> ln X i ) and smaller where we are interested in traversing the current region of prior volume more quickly. This allows us to prioritize adding samples in regions of interest.  \nIn general, we would like the number of live points K ' X ' as a function of prior volume X to follow a particular importance function I' X ' such that  \nK ' X ' / I' X ' (19)  \nWhile this function can be completely general, since most users are interested in estimating the posterior P' GLYPH<2> ' and/or evidence Z more generally, dynesty by default follows Hig-  \nFigure 3. An example highlighting different schemes for live point allocation between Static and Dynamic Nested Sampling run in dynesty with a fixed number of samples . See § 3 for additional details. Top panels: As Figure 2, but now highlighting the number of live points (upper) and evidence estimates (lower) for a Static Nested Sampling run (black) and Dynamic Nested Sampling runs focused entirely on estimating the posterior (blue), entirely on estimating the evidence (green), and with an 80%/20% posterior/evidence mixture (the default in dynesty ; red). Bottom panels: The distribution of samples from the targeted 3-D correlated Gaussian distribution in the Static (left), posterior-focused (middle), and evidence-focused (right) runs. Points are color-coded based on their important weight p i . The posterior-oriented run allocates points almost exclusively around the bulk of the posterior mass, while the evidence-oriented run preferentially allocates them in prior-dominated regions.  \n<!-- image -->  \nson et al. (2017b) and considers a function of the form:  \nI' X ' = f P I P ' X ' + ' 1 GLYPH<0> f P 'I Z ' X ' (20)  \nwhere f P is the relative amount of importance placed on estimating the posterior.  \nWe define the posterior importance function as  \nI P ' X ' GLYPH<17> p ' X ' (21)  \nwhere p ' X ' is the now the probability density function (PDF) of the importance weight defined in § 2.5. This choice just means that we want to allocate more live points in regions where the posterior mass / L' X ' d X is higher.  \nWe define the evidence importance function as  \nI Z ' X ' GLYPH<17> 1 GLYPH<0> Z' X 'GLYPH<157>Z fl 1 0 ' 1 GLYPH<0> Z' X 'GLYPH<157>Z' d X (22)  \nwhere Z' X ' is the evidence integrated up to X . This means that we want to allocate more live points when we believe we have not integrated over much of the posterior (i.e. in the prior volume-dominated regime at larger values of X ) and fewer as we integrate over larger portions of the posterior mass and become more confident in our estimated value of Z (see Figure 2).", 'Algorithm 3: Iterative Dynamic Nested Sampling': "- // Baseline Nested Sampling run.  \nRun Static Nested Sampling (Algorithm 1) with:  \n- (a) K live points\n- (b) sampled uniformly from the prior GLYPH<25> ' GLYPH<2> '\n- (c) until the default Static Nested Sampling stopping criterion is met.\n- // Main sampling loop.  \nwhile stopping criterion not met do  \n// Find region where new samples should be allocated.  \nCompute relative importance f ˆ I' ˆ X i 'g over all dead points f GLYPH<2> i g .  \nUse f ˆ I i g to assign lower L low = L' ˆ X high ' and upper L high = L' ˆ X low ' likelihood bounds.  \n// Batch Nested Sampling run.  \nRun Static Nested Sampling (Algorithm 1) with:  \n- (a) K 0 live points\n- (b) sampled uniformly from the constrained prior GLYPH<25> ' GLYPH<2> ' based on the lower likelihood bound GLYPH<21> = L low  \nGLYPH<21>  \n- (c) until the likelihood L' GLYPH<2> ' of the last dead point exceeds the upper likelihood bound L high .\n- // Merge samples from batch.  \nMerge new batch of dead points f GLYPH<2> 0 i g into the previous set of dead points f GLYPH<2> i g .  \n// Check whether to stop.  \nEvaluate stopping criterion.  \nend", '3.2 Iterative Dynamic Nested Sampling': 'As in § 2.4, we unfortunately do not have access to X or I\' X \' directly. We thus need to use noisy estimators to approximate them, which are only available after we have already generated samples from the posterior . In practice then, Dynamic Nested Sampling works as an iterative modification to Static Nested Sampling. We outline this \'Iterative\' Dynamic Nested Sampling approach, first proposed in Higson et al. (2017b) and implemented in dynesty , in Algorithm 3. It has five main steps:  \n- (i) Sample the distribution with Static Nested Sampling to lay down a \'baseline run\' to get a sense where the posterior mass P\' X \' d X is located.\n- (ii) Evaluate our importance function I\' X \' over the existing set of samples.\n- (iii) Use the computed importances I i to decide where to allocate additional live points/samples.\n- (iv) Add a new\'batch\'of samples in the region of interest using Static Nested Sampling.\n- (v) \'Merge\' the new batch of samples into the previous set of samples.  \nWe then repeat steps (ii) to (v) until some stopping criterion is met. By default, dynesty uses K base = K batch = 250 points for each run, although this should be adjusted depending on the problem at hand.  \nAllocating points using an existing set of samples is a two-step process. First, we evaluate a noisy estimate of our importance function over the samples:  \nˆ I i = f P ˆ p i " N i = 1 ˆ p i + \' 1 GLYPH<0> f P \' 1 GLYPH<0> ˆ Z i GLYPH<157>\' ˆ Z N + GLYPH<1> ˆ Z N \' " N i = 1 1 GLYPH<0> ˆ Z i GLYPH<157>\' ˆ Z N + GLYPH<1> ˆ Z N \' (23)  \nwhere we are now using the noisy importance weight ˆ p i to estimate the posterior and the rough upper limit GLYPH<1> ˆ Z N GLYPH<24> L max N ˆ X N to estimate the remaining evidence. Then, we use these values to define new regions of prior volume to sample. By default, dynesty only samples from a single contiguous  \nrange of prior volume \' X low ; X high … which define an associated (flipped) range in iteration » i low ; i high \' and likelihood »L low ; L high \' defined by the simple heuristic  \ni low = min GLYPH<2> min \'f i g\' GLYPH<0> n pad ; 0 GLYPH<3> i high = max GLYPH<2> max \'f i g\' + n pad ; N GLYPH<3> (24) 8 i 2 » 0 ; N … s : t : ˆ I i GLYPH<21> f max GLYPH<2> max \'f ˆ I i g\'  \nwhere f max serves as a threshold relative to the peak value and n pad pads the starting/ending iteration. In other words, we compute the importance values ˆ I i over the existing set of samples, compute the minimum i low and maximum i high iterations where the importance is above a threshold f max relative to the peak, and shift the final values by n pad . By default, dynesty assumes f P = 0 : 8 (80% posterior vs 20% evidence), f max = 0 : 8 (80% thresholding), and n pad = 1 .  \nOnce we have computed » i low ; i high … , we can then just start a new Static Nested Sampling run that samples from the constrained prior between »L low ; L high \' . In the case where L low = 0 , this is just the original prior GLYPH<25> \' GLYPH<2> \' and our Static Nested Sampling run is identical to Algorithm 1 except with stopping criteria L\' GLYPH<2> \' GLYPH<21> L high . If L low > 0 , however, then we are instead starting interior to the prior and thus not fully integrating over it. So while those new samples will improve the relative posterior resolution GLYPH<1> ln X i and thus the posterior estimate ˆ P\' GLYPH<2> \' , they will not actually improve the evidence estimate ˆ Z .  \nFinally, we need to \'merge\' our new set of N 0 samples f GLYPH<2> 0 1 ; : : : ; GLYPH<2> 0 N 0 g into our original set of samples f GLYPH<2> j g . This process is straightforward and can be accomplished following the procedure outlined in Appendix A. We are then left with a combined set of samples f GLYPH<2> 1 ; : : : ; GLYPH<2> N + N 0 g with new associated prior volumes f X 1 ; : : : ; X N + N 0 g and a variable number of live points f K 1 ; : : : ; K N + N 0 g at every iteration.', '3.3 Estimating the Prior Volume': 'As shown in Appendix A, we can reinterpret the results from § 2.3 as a consequence of the two different ways Nested Sampling traverses the prior volume. In the first case, where the number of live points K i GLYPH<21> K i GLYPH<0> 1 increases or stays the same, we know that we have (possibly) added live points and then replaced the one with the lowest likelihood L min . In this case, the prior volume experiences exponential shrinkage such that  \nE GLYPH<2> GLYPH<1> ln ˆ X i GLYPH<3> = GLYPH<0> 1 K i (25)  \nIn the second case, where the number of live points K j + 1 < K j strictly decreases, we know that we have removed the live point(s) with the lowest likelihood L min . For each of the k iterations where this continues to occur, the prior volume experiences uniform shrinkage such that  \nE " GLYPH<1> ˆ X j + k ˆ X j # = 1 K j + 1 (26)  \nIn Static Nested Sampling, these two regimes are cleanly divided, with the main set of dead points traversing the prior volume exponentially and the final set of \'recycled\' live points traversing it uniformly. In Dynamic Nested Sampling, however, we are constantly switching between exponential and uniform shrinkage as we increase or decrease the number of live points at a given iteration.', '3.4 Stopping Criterion': "The implementation of Static Nested Sampling outlined in Algorithm 1 generally exclusively targets evidence estimation. This gives a natural stopping criterion (see § 2.4) to terminate sampling once we believe that we have integrated over a majority of the posterior P' GLYPH<2> ' such that additional samples will no longer improve our evidence estimate ˆ Z .  \nIn the Dynamic Nested Sampling case, however, we are no longer just interested in computing the evidence. Because we now have the flexibility to vary the number of live points K i over time, we are also interested in the properties of our integral (and the samples that comprise the integrand) in addition to the question of whether our integral has converged.  \nThis flexibility necessitates the introduction of more complex stopping criteria to assess whether those alternative properties are behaving as expected. Similar to § 3.1, we consider a stopping criteria of the form:  \nS = s P S P + ' 1 GLYPH<0> s P 'S Z < GLYPH<15> (27)  \nwhere GLYPH<15> is our tolerance, S P is the posterior stopping criterion, S Z is the evidence stopping criterion, and s P is the relative amount of weight given to S P over S Z .  \nWe define our stopping criterion to be the amount of fractional uncertainty in the current posterior ˆ P' GLYPH<2> ' and evidence ˆ Z estimates. For the posterior P' GLYPH<2> ' , we start by defining 'posterior noise' to be the Kullback-Leibler (KL) divergence  \nH ' ˆ P 0 j j ˆ P' GLYPH<17> E ˆ P 0 GLYPH<2> ln ˆ P 0 GLYPH<0> ln ˆ P GLYPH<3> (28)  \n= ' GLYPH<10> GLYPH<2> ˆ P 0 ' GLYPH<2> ' ln ˆ P 0 ' GLYPH<2> ' d GLYPH<2> GLYPH<0> ' GLYPH<10> GLYPH<2> ˆ P 0 ' GLYPH<2> ' ln ˆ P' GLYPH<2> ' d GLYPH<2> (29)  \nbetween the posterior estimate ˆ P 0 ' GLYPH<2> ' from a random hypothetical Nested Sampling run with the same setup and our current estimate ˆ P' GLYPH<2> ' . This can be interpreted as the 'information loss' due to random noise in our posterior estimate ˆ P' GLYPH<2> ' . Our proposed posterior stopping criteria is then  \nS P GLYPH<17> 1 GLYPH<15> P GLYPH<27> GLYPH<2> H ' ˆ P 0 j j ˆ P' GLYPH<3> E GLYPH<2> H ' ˆ P 0 j j ˆ P' GLYPH<3> (30)  \nwhere GLYPH<15> P normalizes the posterior deviation to a desired scale. For the evidence Z , this is just the estimated fractional scatter between the evidence estimates ˆ Z 0 from random hypothetical Nested Sampling runs with the same setup. Following Higson et al. (2017b), we opt to compute this in logspace for convenience:  \nS Z GLYPH<17> 1 GLYPH<15> Z GLYPH<27> h ln ˆ Z 0 i (31)  \nwhere GLYPH<15> Z normalizes the evidence deviation to a desired scale.  \nUnsurprisingly, we do not have access to the distribution of all hypothetical Nested Sampling runs with the same setup to compute these exact estimates. However, as with § 2.4 and § 3.2, we do have access to noisy estimates of these quantities via procedures described in Higson et al. (2017a) and outlined in Appendix A for simulating Nested Sampling errors. dynesty uses M simulated values of these noisy estimates to estimate the stopping criteria as:  \nˆ S = s P GLYPH<15> P GLYPH<27> GLYPH<2> f ˆ H 1 ; : : : ; ˆ H M g GLYPH<3> E GLYPH<2> f ˆ H 1 ; : : : ; ˆ H M g GLYPH<3> + ' 1 GLYPH<0> s P ' GLYPH<15> Z GLYPH<27> h f ln ˆ ˆ Z 1 ; : : : ; ln ˆ ˆ Z M g i (32)  \nwhere the ˆ ˆ Z notation just emphasizes that we are constructing a noisy estimator of our already-noisy estimate ˆ Z . By default, dynesty assumes s P = 1 (100% focused on reducing posterior noise), GLYPH<15> = 1 , GLYPH<15> P = 0 : 02 , GLYPH<15> Z = 0 : 1 , and M = 128 .", '4 IMPLEMENTATION': 'Now that we have outlined the basic algorithm and approach behind Dynamic Nested Sampling, we now turn our attention to the problem of generating samples from the constrained prior. dynesty approaches this problem in two parts:  \n- (i) constructing appropriate bounding distributions that encompass the remaining prior volume over multiple possible modes and\n- (ii) proposing new live points by generating samples conditioned on these bounds.  \ndynesty contains several options for both constructing bounds and sampling conditioned on them. We provide an broad overview of each of these in turn.', '4.1 Bounding Distributions': 'In general, dynesty tries to use the distribution of the current set of live points to try and get a rough idea of the shape and size of the various regions of prior volume that  \nFigure 4. An example highlighting the various bounding distributions implemented in dynesty . These include the entire unit cube (left), a single ellipsoid (left-middle), multiple overlapping ellipsoids (middle), overlapping spheres (right-middle), and overlapping cubes (right). The current set of live points are shown in purple while draws from the bounding distribution are shown in grey. A schematic representation of each bounding distribution is shown in the bottom-right-hand corner of each panel. See § 4.1 for additional details.  \n<!-- image -->  \nUnit Cube (no bound)  \n<!-- image -->', 'Bounding Distributions': "<!-- image -->  \n<!-- image -->  \nSingle Ellipsoid  \n<!-- image -->  \nMultiple Ellipsoids  \nOverlapping  \nOverlapping Cubes  \nwe are currently sampling. These are then used to condition various sampling methods to try and improve the efficiency. There are five bounding methods currently implemented in dynesty :  \n- GLYPH<15> no bounds (i.e. the unit cube),\n- GLYPH<15> a single ellipsoid,\n- GLYPH<15> multiple ellipsoids,\n- GLYPH<15> many overlapping balls, and\n- GLYPH<15> many overlapping cubes.  \nIn general, single ellipsoids tend to perform reasonably well at estimating structure when the likelihood is roughly Gaussian and uni-modal. In more complex cases, however, decomposing the live points into separate clusters with their own bounding ellipsoids works reasonably well at locating and tracking structure. In low ( D glyph[lessorsimilar] 5 ) dimensions, allowing the live points themselves to define emergent structure through many overlapping balls or cubes can perform better provided the L' GLYPH<2> ' spans similar scales in each of the parameters. Finally, using no bounds at all is only recommended as an option of last resort and is mostly relevant when performing systematics checks or if the number of live points K GLYPH<28> D 2 GLYPH<157> 2 is small relative to the number of possible parameter covariances.  \nIn addition to these various options, dynesty also tries to increase the volume of all bounds by a factor GLYPH<11> to be conservative about the size of the constrained prior. While this is generally assumed to take a constant value of GLYPH<11> = 1 : 25 , it can also be derived 'on the fly' using bootstrapping methods following the approach outlined in Buchner (2016). Deriving accurate volume expansion factors are extremely important when sampling uniformly but are less relevant for other sampling schemes that are more robust to the exact sizes of the bounds (see § 4.2).  \nBy default, dynesty uses multiple ellipsoids to construct the bounding distribution. A summary of the various bounding methods can be found in Figure 4. We describe these each in turn below.", '4.1.1 Unit Cube': "The simplest case of using the entire unit cube (i.e. simple rejection sampling over the entire prior GLYPH<25> ' GLYPH<2> ' with no limits) can be useful in a few edge cases where the number of live points K is small compared to the number of dimensions D , or where users are interested in performing tests to verify sampling behavior.", '4.1.2 Single Ellipsoid': 'As shown in (Mukherjee et al. 2006), a single bounding ellipsoid can be effective if the posterior is unimodal and roughly Gaussian. dynesty uses a scaled version of the empirical covariance matrix C 0 = GLYPH<13> C centered on the empirical mean GLYPH<22> of the current set of live points to determine the size and shape of the ellipsoid, where GLYPH<13> is set so the ellipsoid encompasses all available live points.', '4.1.3 Multiple Ellipsoids': 'By default, dynesty does not assume the posterior is unimodal or Gaussian and instead tries to bound the live points using a set of (possibly overlapping) ellipsoids. These are constructed using an iterative clustering scheme following the algorithm outlined in Shaw et al. (2007) and Feroz & Hobson (2008) and implemented in the online package nestle . 5 In brief, we start by constructing a bounding ellipsoid over the entire collection of live points. We then initialize 2 k -means clusters at the endpoints of the major axes, optimize their positions, assign live points to each cluster, and construct a new pair of bounding ellipsoids for each new cluster of live points. The decomposition is accepted if the combined volume of the subsequent pair of ellipsoids is substantially smaller. This process is then performed recursively until no decomposition is accepted.', '12 J. S. Speagle': "By default, dynesty tries to be substantially more conservative when decomposing live points into separate clusters and bounding ellipsoids than alternative approaches used in MultiNest (Feroz & Hobson 2008; Feroz et al. 2013). This algorithmic choice, which can substantially reduce the overall sampling efficiency, is made in order to avoid 'shredding' the posterior into many tiny islands of isolated live point clusters. As shown in Buchner (2016), that behavior can lead to biases in the estimated evidence ˆ Z and posterior ˆ P' GLYPH<2> ' .", '4.1.4 Overlapping Balls': 'An alternate approach to using bounding ellipsoids is to allow the current set of live points themselves to define emergent structure. The simplest approach used in dynesty follows Buchner (2016, 2017) by assigning a D -dimensional ball (sphere) with radius r to each live point, where r is set using bootstrapping and/or leave-one-out techniques to encompass GLYPH<21> 1 other live points. One benefit to this approach over using multiple ellipsoids (which can depend sensitively on the clustering schemes) is that it is almost entirely free of tuning parameters, with the overall behavior only weakly dependent on the number of bootstrap realizations.', '4.1.5 Overlapping Cubes': "As with the set of overlapping balls, dynesty also implements a similar algorithm based on Buchner (2016, 2017) involving overlapping cubes with half-side-length ' . As § 4.1.4, ' is derived using either bootstrapping and/or leave-one-out techniques so that the cubes encompass GLYPH<21> 1 other live points.", '4.2 Sampling Methods': "Once a bounding distribution has been constructed, dynesty generates samples conditioned on those bounds. In general, this follows a strategy of  \nf ' s C b ; GLYPH<2> ' ! GLYPH<2> 0 (33)  \nwhere C b is the covariance associated with a particular bound b (e.g., an ellipsoid), GLYPH<2> is the starting position, GLYPH<2> 0 is the final proposed position, and s GLYPH<24> 1 is a scale-factor that is adaptively tuned over the course of a run to ensure optimal acceptance rates.  \ndynesty implements four main approaches to generating samples:  \n- GLYPH<15> uniform sampling,\n- GLYPH<15> random walks,\n- GLYPH<15> multivariate slice sampling, and\n- GLYPH<15> Hamiltonian slice sampling.  \nThese each are designed for different regimes. Uniform sampling can be relatively efficient in lower dimensions where the bounding distribution can approximate the prior volume better but struggles in higher dimensions since it is extremely sensitive to the size of the bounds. Random walks are less sensitive to the size of the bounding distribution and so tend to work better than uniform sampling in moderate dimensional spaces but still struggle in high-dimensional  \nspaces because of the exponentially increasing amount of volume it needs to explore. Multivariate and Hamiltonain slice sampling often performs better in these high-dimensional regimes by avoiding sampling directly from the volume and taking advantage of gradients, respectively.  \nIn addition to each method's performance in various regimes, there is also a fundamental qualitative difference between uniform sampling and the other sampling approaches outlined above. Uniform sampling, by construction, can only sample directly from the bounding distribution. This makes it uniquely sensitive to the assumption that the bounds entirely encompass the current prior volume at a given iteration, which is never fully guaranteed (Buchner 2016). By contrast, the other sampling methods are MCMCbased: they generate samples by 'evolving' a current live point to a new position. This allows them to generate samples outside the bounding distribution, making them less sensitive to this assumption.  \nBy default, dynesty resorts to uniform sampling when the number of dimensions D < 10 , random walks when 10 GLYPH<20> D GLYPH<20> 20 , and Hamiltonian/multivariate slice sampling when D > 20 if a gradient is/is not provided. A summary of the various sampling methods can be found in Figure 5. We describe these each in turn below.", '4.2.1 Uniform Sampling': 'If we assume that our bounding distribution B \' GLYPH<2> \' encloses the constrained prior GLYPH<25> GLYPH<21> \' GLYPH<2> \' , the most direct approach to generating samples from the bounds is to sample from them uniformly. This procedure by construction produces entirely independent samples between each iteration i , and tends to work best when the volume of the bounds X B \' GLYPH<21> \' is roughly the same order of magnitude as the current prior volume X \' GLYPH<21> \' (leading to glyph[greaterorsimilar] 10% acceptance rates).  \nIn general, the procedure for generating uniform samples from overlapping bounds is straightforward (see, e.g., Feroz & Hobson 2008; Buchner 2016):  \n- (i) Pick a bound b at random with probability p b / X b proportional to its volume X b .\n- (ii) Sample a point GLYPH<2> b uniformly from the bound.\n- (iii) Accept the point with probability 1 GLYPH<157> q , where q GLYPH<21> 1 is the number of bounds GLYPH<2> b lies within.  \nThis approach ensures that any proposed sample will be drawn from the bounding distributing B \' GLYPH<2> \' comprised of the union of all bounds, which has a volume X B GLYPH<20> " N b b = 1 X b that is strictly less than or equal to the sum of the volumes of each individual bound.  \nGenerating samples uniformly from the bounds in § 4.1 falls into two cases: cubes and ellipsoids. Generating points from an D -cube centered at GLYPH<2> b with half-side-length \' is trivial and can be accomplished via:  \n- (i) Generate D iid uniform random numbers U = f U 1 ; : : : ; U D g from »GLYPH<0> \'; \' … .\n- (ii) Set GLYPH<2> 0 = U + GLYPH<2> b .  \nGenerating points from an ellipsoid centered at GLYPH<2> b with covariance C b with matrix square-root C 1 GLYPH<157> 2 b is also straightforward but slightly more involved:', 'Random Walk': 'Figure 5. A schematic illustration of the different sampling methods implemented in dynesty . These include: uniform sampling from the bounding distribution (top-left), random walks proposals starting from a random live point based on the bounding distribution (top-right) with either fixed or variable scale-lengths for proposals, multivariate slice sampling proposals starting from a random live point (bottom-left) using either the principle axes or a random direction sampled from the bounding distribution, and Hamiltonian slice sampling away from a random live point forwards and backwards in time (bottom-right). See § 4.2 for additional details.  \n<!-- image -->  \n<!-- image -->  \nFixed Scale  \n<!-- image -->  \nVariable Scale', 'Hamiltonian Slice': "Principle Axes  \n<!-- image -->  \nRandom  \n<!-- image -->  \n- (i) Generate D iid standard normal random numbers Z = f Z 1 ; : : : Z D g .\n- (ii) Compute the normalized vector V GLYPH<17> Z GLYPH<157> j j Z j j .\n- (iii) Draw a standard uniform random number U and compute S GLYPH<17> U D V .\n- (iv) Set GLYPH<2> = C 1 GLYPH<157> 2 b S + GLYPH<2> b .  \nStep (ii) creates a random vector V that is uniformly distributed on the surface of the D -sphere. Step (iii) randomly moves V ! S to an interior radius r 2 ' 0 ; 1 ' based on the fact that the volume of a D -sphere scales as V ' r ' / r D . Finally, step (iv) adjusts the scale, shape, and center to match that of the bounding ellipsoid.", '4.2.2 Random Walks': "An alternative approach to sampling uniformly within the bounding distribution B ' GLYPH<2> ' is to instead to try and propose new positions by 'evolving' a given live point GLYPH<2> k ! GLYPH<2> 0 to a new position. Since L' GLYPH<2> k ' GLYPH<21> L min i at a given iteration by definition, this procedure also guarantees that we will be generating samples exclusively within the constrained prior GLYPH<25> ' GLYPH<2> ' .  \nGLYPH<21>  \nOne straightforward approach to 'evolving' a live point to a new position is to consider sampling from the con-  \nstrained prior using a simple Metroplis-Hastings (MH; Metropolis et al. 1953; Hastings 1970) MCMC algorithm:  \n- (i) Propose a new position GLYPH<2> 0 GLYPH<24> Q ' GLYPH<2> j GLYPH<2> k ' from the proposal distribution Q ' GLYPH<2> j GLYPH<2> k ' starting from GLYPH<2> k .\n- (ii) Move to GLYPH<2> 0 with probability A = GLYPH<25> GLYPH<21> ' GLYPH<2> 0 ' GLYPH<25> GLYPH<21> ' GLYPH<2> k ' Q ' GLYPH<2> k j GLYPH<2> 0 ' Q ' GLYPH<2> 0 j GLYPH<2> k ' . Otherwise, stay at GLYPH<2> k .\n- (iii) Repeat (i)-(ii) for N walks iterations.  \nSince the constrained prior is flat (see § 2.2), the ratio of the constrained prior values is by definition 1 . Likewise, if we choose a symmetric proposal distribution Q ' GLYPH<2> j GLYPH<2> k ' , then the ratio of the proposal distributions also evaluates to 1 . This procedure then reduces to simply accepting a new point if it is within the constrained prior with L' GLYPH<2> i ' GLYPH<21> GLYPH<21> and rejecting it otherwise. By default, dynesty takes N walks = 25 .  \ndynesty implements two forms of the proposal Q ' GLYPH<2> j GLYPH<2> k ' . The default option is to propose new positions uniformly from an associated ellipsoid centered on GLYPH<2> k with covariance C b , where C b is one of the bounding distributions that encompasses GLYPH<2> k (selected randomly). The second follows the same form as the first, except the covariance C b is re-scaled at each subsequent proposal t GLYPH<20> N walks by GLYPH<13> following the  \n<!-- image -->  \nprocedure outlined in Sivia & Skilling (2006):  \nGLYPH<11> ' t ' = 8 > > > < > > > e 1 GLYPH<157> N acc ' t ' GLYPH<2> GLYPH<13> ' t GLYPH<0> 1 ' N acc ' t ' t > f acc e GLYPH<0> 1 GLYPH<157> N rej ' t ' GLYPH<2> GLYPH<13> ' t GLYPH<0> 1 ' N acc ' t ' t < f acc GLYPH<13> ' t GLYPH<0> 1 ' N acc t = f acc (34)  \n:  \nwhere N acc ' t ' and N rej ' t ' is the total number of accepted and rejected proposals by iteration t , respectively, f acc is the desired acceptance fraction, and GLYPH<13> ' t = 0 ' = 1 . By default, dynesty targets f acc = 0 : 5 .", '4.2.3 Multivariate Slice Sampling': "In higher dimensions, rejection sampling-based methods such as the random walk proposals outlined in § 4.2.2 can become progressively more inefficient. To remedy this, dynesty includes slice sampling (Neal 2003) routines designed to sample from the constrained prior GLYPH<25> GLYPH<21> ' GLYPH<2> ' . These are based on the 'stepping out' method proposed in Neal (2003) and Jasa & Xiang (2012), which works as follows in the single-variable case starting from the position x k of the k th live point:  \n- (i) Draw a standard uniform random number U .\n- (ii) Set the left bound L = x k GLYPH<0> w U and the right as R = L + w where w is the starting 'window'.  \n(iii) While L' L ' GLYPH<21> GLYPH<21> , extend the position of the left bound L by w . Repeat this procedure for R .  \n- (iv) Sample a point x 0 GLYPH<24> Unif ' L ; R ' uniformly on the interval from L to R .\n- (v) If L' x 0 ' > GLYPH<21> , accept x 0 . Otherwise, reassign the corresponding bound to be x 0 ( L if x 0 < x and R otherwise) and repeat steps (iv)-(v).  \nWhen sampling in higher dimensions, the single-variable update outlined above can be interpreted as a Gibbs sampling update (Geman & Geman 1987) where instead of drawing GLYPH<2> directly we instead update each component in turn  \nGLYPH<2> 0 GLYPH<24> GLYPH<25> GLYPH<21> ' GLYPH<2> ' ) 8 > > > > < > > > > GLYPH<2> 1 GLYPH<24> GLYPH<25> GLYPH<21> ' GLYPH<2> 1 j GLYPH<2> n 1 ' : : : GLYPH<2> D GLYPH<24> GLYPH<25> GLYPH<21> ' GLYPH<2> D j GLYPH<2> n D ' (35)  \n:  \nwhere GLYPH<2> n i are the set of D GLYPH<0> 1 parameters excluding GLYPH<2> i . We then repeat this procedure for N slices iterations. By default dynesty takes N slices = 5 .  \nThis procedure is generally robust, although it can introduce longer correlation times if there are strong covariances between parameters. To mitigate this, dynesty by default executes single-variable slice sampling updates along the principle axes V b GLYPH<17> f v 1 ; b ; : : : ; v D ; b g associated with the covariance C b from a given bound b . This allows us to automatically set both the direction v i ; b and associated scale j j v i ; b j j of the window while trying to reduce the correlations among sets of parameters.  \nAlternately, instead of executing a full Gibbs update by rotating through the entire set of parameters in turn, we can sample along a random trajectory v 0 through the prior instead. This procedure is similar to that implemented in PolyChord (Handley et al. 2015), except that rather than 'whitening' the set of live points using the associated C b we instead draw v 0 from the surface of the corresponding bound with covariance C b . Provided a suitable number of N slices GLYPH<24> D , this procedure also can generate suitably independent new positions GLYPH<2> 0 .  \nOver the past two decades, sampling methods have increasingly attempted to incorporate gradients to improve their overall performance, especially in high-dimensional spaces. The most common class of methods are based on Hamiltonian Monte Carlo (HMC; Neal 2012; Betancourt 2017), whereby a particle at a given position x is assigned a mass matrix M and some momentum p and allowed to sample from the joint distribution  \nP ' x ; p j M ' / exp »GLYPH<0>H ' x ; p j M '… (36)  \nwhere  \nH' x ; p j M ' GLYPH<17> U ' x ' + K ' p j M ' GLYPH<17> GLYPH<0> ln » GLYPH<25> ' x 'L' x '… + 1 2 p T M GLYPH<0> 1 p (37)  \nis the Hamiltonian of the system with a 'potential energy' U ' x ' and 'kinetic energy' K ' p j M ' , and T is the transpose operator. Typically, proposals are generated by sampling the momentum from the corresponding multivariate Normal (Gaussian) distribution  \np GLYPH<24> N' 0 ; M ' ; (38)  \nwith mean 0 and covariance M , evolving the system via Hamilton's equations from H' x ; p ' ! H' x 0 ; p 0 ' , and then accepting the new position based on the MH acceptance criteria outlined in § 4.2.2. In other words, at each iteration we randomly assign a given particle some mass and velocity and then have it explore the potential defined by the (log-)posterior.  \nAs with the previous methods, this approach simplifies dramatically when sampling over the constrained prior GLYPH<25> GLYPH<21> ' GLYPH<2> ' . In that case, since the distribution is flat, the momentum remains unchanged until the particle hits the hard likelihood boundary, at which point it reflects so that  \np 0 = p GLYPH<0> 2 h p GLYPH<1> h j j h j j 2 (39)  \nwhere h is the gradient at the point of reflection. This version of the algorithm is referred to elsewhere as Galilean Monte Carlo (Skilling 2012; Feroz & Skilling 2013) or reflective slice sampling (Neal 2003).  \nIn practice, since we have to evolve the system discretely, there are a few additional caveats to consider. Most importantly, the use of discrete time-steps means reflection will not occur right at the boundary of the constrained prior but slightly beyond it, which does not guarantee reflections will end up back inside the constrained prior. This behavior, which arises from larger time-steps, 'terminates' the particle's trajectory in that particular direction and leads to inefficient sampling that isn't able to explore the full parameter space.  \nOn the other hand, using extremely small time-steps means spending the vast majority of time evaluating positions along a straight line, which is also non-optimal. dynesty by default attempts to compromise between these two behaviors by optimizing the time-step so that f move GLYPH<24> 0 : 9 of total steps are spent moving forward passively instead of reflecting or terminating. In addition, dynesty by default caps the total number of time-steps to N move = 100 to prevent trajectories from being evolved indefinitely.  \nSimilar to algorithms such as the No U-Turn Sampler (NUTS; Hoffman & Gelman 2011), dynesty also considers trajectories evolved forwards and backwards in time to  \n<!-- image -->  \n<!-- image -->  \nFigure 6. Illustration of dynesty 's performance using multiple bounding ellipsoids and uniform sampling over 2-D Gaussian shells (highlighted in Figure 4) meant to test the code's bounding distributions. Left : A smoothed corner plot showing the exact 1-D and 2-D marginalized posteriors of the target distribution. Middle : As before, but now showing the final distribution of weighted samples. Right: The volume of the bounding distribution when using a single ellipsoid (blue) versus multiple ellipsoids (orange) over the course of the run. Since a single ellipsoid is a poor model for this distribution, its volume quickly saturates as it becomes unable to accurately capture the distribution of live points. Allowing the bounding distribution to be modeled by multiple ellipsoids allows for dynesty to capture the more complex structure as the live points move increasingly into organized rings.  \n<!-- image -->  \nbroaden the range of possible positions explored in a given proposal. While these roughly double the number of overall time-steps, they substantially improve overall behavior by exploring larger regions of the constrained prior.  \ndynesty employs two additional schemes to try and further mitigate discretization effects on the sampling procedure described above. First, the time-step used at a given iteration is allowed to vary randomly by up to 30% following recommendations from Neal (2012). This helps to suppress resonant behavior that can arise from poor choices of timesteps without substantially impacting overall performance. Second, rather than merely accepting positions at the end of a trajectory, dynesty instead tries to sample uniformly from the entire trajectory by treating it as a set of slices defined by ' GLYPH<2> i L ; GLYPH<2> i ; GLYPH<2> i R ' left-inner-right position tuples. New samples are then proposed via the following scheme:  \n- (i) Compute the length ' i of each line segment ' GLYPH<2> i L ; GLYPH<2> i R ' .\n- (ii) Selecting a line segment i at random proportional to its length.  \n(iii) Sample a point GLYPH<2> 0 uniformly on the line segment defined by ' GLYPH<2> i L ; GLYPH<2> i R ' .  \n(iv) If L' GLYPH<2> 0 ' > GLYPH<21> , accept GLYPH<2> 0 . Otherwise, reassign the corresponding bound to be GLYPH<2> 0 ( GLYPH<2> i L if GLYPH<2> 0 is on the line segment » GLYPH<2> i L ; GLYPH<2> i ' and GLYPH<2> i R otherwise) and repeat steps (i)-(iv).  \nWhile there are a variety of possible approaches to applying HMC-like methods to Nested Sampling other than the basic procedure outlined above, we defer any detailed comparisons between them to possible future work.", '5 TESTS': "Here, we examine dynesty 's performance on a variety of toy problems designed to stress-test various aspects of the code. Additional tests can also be found online.", '5.1 Gaussian Shells': 'One standard problem that tests the efficiency of the ability of bounding distributions to transition between a flat surface to separated, elongated structures is the D -dimensional \'Gaussian shells\' from Feroz & Hobson (2008). The likelihood of the distribution is defined as  \nL\' GLYPH<2> \' = circ \' GLYPH<2> j c 1 ; r 1 ; w 1 \' + circ \' GLYPH<2> j c 2 ; r 2 ; w 2 \' (40)  \nwhere  \ncirc \' GLYPH<2> j c ; r ; w \' = 1 p 2 GLYPH<25> w 2 exp " GLYPH<0> 1 2 \'jj GLYPH<2> GLYPH<0> c j j GLYPH<0> r \' 2 w 2 # (41)  \nFollowing Feroz et al. (2013), we take the centers c 1 and c 2 of the two positions to be GLYPH<0> 3 : 5 and 3 : 5 in the first dimension and 0 in all others, respectively, the radius r = 2 , and the width w = 0 : 1 . Our prior is defined to be uniform from »GLYPH<0> 6 ; 6 … to encompass the majority of the likelihood and ensure a smooth transition between the uni-modal starting distribution and the multi-modal target distribution.  \nWe illustrate dynesty \'s performance in the 2-D case in Figure 6. The default configuration options in dynesty (multiple ellipsoid bounds with uniform sampling) lead to a roughly 10% sampling efficiency over the course of GLYPH<24> 20 k iterations when using Dynamic Nested Sampling and lead to excellent posterior estimates. We also see that the multiellipsoidal decomposition algorithm works as expected, with the total volume of the bounding distribution decreasing dramatically as the live points begin to organize themselves within the two shells.', '5.2 Eggbox': "Another distribution we consider to test the ability of dynesty to track and evolve multiple modes is the 2-D'Eggbox' likelihood from Feroz & Hobson (2008), which we defined as  \nL' x ; y ' = exp ( GLYPH<20> 2 + cos GLYPH<18> 5 GLYPH<25> ' x GLYPH<0> 1 ' 2 GLYPH<19> sin GLYPH<18> 5 GLYPH<25> ' y GLYPH<0> 1 ' 2 GLYPH<19> GLYPH<21> 5 ) (42)  \nFigure 7. Illustration of dynesty 's performance using multiple bounding ellipsoids and overlapping balls with uniform sampling over the 2-D 'Eggbox' distribution meant to test the code's bounding distributions. Top left : The true log-likelihood surface of the Eggbox distribution. Top right : A smoothed corner plot showing the 1-D and 2-D marginalized posteriors of the final distribution of weighted samples from a posterior-oriented Dynamic Nested Samplig run. Bottom : The importance weight PDF p ' X ' (top) and corresponding evidence estimate ˆ Z with 1, 2, and 3-sigma uncertainties (bottom) from two independent evidence-oriented Dynamic Nested Sampling runs using multiple ellipsoids (blue) and overlapping balls (red) as bounding distributions.  \n<!-- image -->  \nThis distribution is periodic over the 2-D unit cube, with 13 localized modes contained within a given period. We take our prior to be standard uniform in x and y to limit sampling to one period.  \nThe resulting posterior and evidence estimates from several posterior-oriented and evidence-oriented Dynamic Nested Sampling runs are shown in Figure 7. dynesty is able to sample from this distribution quite effectively, with average sampling efficiencies ranging from 20 GLYPH<0> 40% when sampling uniformly from the multiple ellipsoids or overlapping balls.", '5.3 Exponential Wave': "We next apply dynesty to a signal reconstruction problem with multiple modes and periodic boundary conditions. Our model is a transformed periodic single from 0 to 2 GLYPH<25> :  \ny ' x ' = exp » n a sin ' f a x + p a ' + n b sin ' f b x + p b '… (43)  \nwhere we observe noisy data points drawn from  \nˆ y ' x ' GLYPH<24> N GLYPH<16> y ' x ' ; GLYPH<27> 2 GLYPH<17> (44)  \nThe likelihood for this model is Gaussian over the corresponding observed datapoints such that  \nln L' GLYPH<2> ' = GLYPH<0> 1 2 N GLYPH<213> i = 1 ln GLYPH<16> 2 GLYPH<25> GLYPH<27> 2 GLYPH<17> + » ˆ y i GLYPH<0> y ' x i j GLYPH<2> '… 2 GLYPH<27> 2 (45)  \nFigure 8. Illustration of dynesty 's performance using multiple bounding ellipsoids and multivariate slice sampling over principle axes to model an 'Exponential Wave' signal meant to test the code's bounding distributions and incorporation of periodic boundary conditions. Left : Trace plots showing the 1-D positions of samples (dead points) over the course of the run, colored by their estimated importance weight PDF p ' X ' . The true model parameters are shown highlighted in red. We see that even though the underlying structure of the distribution spans many different scales and emerges in different stages, dynesty is able to confidently identify the final two modes and converge to the underlying model parameters. Middle : A corner plot showing the 1-D and 2-D marginalized posteriors from the distribution of the final weighted samples. The true model parameter values are shown in red. The 2.5%, 50%, and 97.5% percentiles (i.e. the 2-sigma credible region) are shown as vertical dashed lines. Top right : The noisy data (gray crosses) and underlying model (red points).  \n<!-- image -->  \nand has seven free parameters: two controlling the relevant amplitudes ( n a , n b ), two controlling the frequencies ( f a , f b ), two controlling the phases ( p a , p b ), and one controlling the scatter GLYPH<27> .  \nWe take our true model parameters to be f a = 1 : 05 , f b = 4 : 2 , n a = 0 : 8 , n b = 0 : 3 , p a = 0 : 1 , p b = 2 : 4 , and GLYPH<27> = 0 : 2 so that a solution is close to the boundary. We assign our prior to be uniform or log-uniform in all parameters with log n a 2 »GLYPH<0> 2 ; 2 ' , log f a 2 »GLYPH<0> 2 ; 2 ' , p a 2 » 0 ; 2 GLYPH<25> ' , log n b 2 »GLYPH<0> 2 ; 2 ' , log f b 2 »GLYPH<0> 2 ; 2 ' , p b 2 » 0 ; 2 GLYPH<25> ' , and log GLYPH<27> 2 »GLYPH<0> 2 ; 0 ' , where the priors in p a and p b are periodic.  \nWe illustrate dynesty 's performance on this problem in Figure 8. We find dynesty is able to robustly recover both modes in this problem, including the solution near the boundary.", '5.4 200-D Gaussian': "We next examine dynesty 's behavior in higher dimensions by testing its performance on a 200-D multivariate Gaussian likelihood with mean GLYPH<22> = 0 and covariance C = I where I is the identity matrix. We assign an identical prior (iid  \nGaussian with GLYPH<22> = 0 and C = I ), such that the posterior will also be iid Gaussian with mean GLYPH<22> = 0 but with covariance C = ' 1 GLYPH<157> 2 ' I .  \nWe sample from this distribution using Hamiltonian Slice Sampling with the analytic log-likelihood gradient. To further highlight the efficiency of these proposals to explore the posterior, we use a small ( K = 50 ) number of live points so that we are highly undersampled relative to the 200-D space. Since dynesty by default uses the empirical covariance (i.e. the MLE estimate) to construct any bounding ellipsoids, this process is dominated by shot noise that can substantially affect the covariance. We consequently impose no bounding distribution (which happens to also be optimal for this problem).  \nAs shown in Figure 9, we find dynesty is able to achieve unbiased recovery of the mean, covariance, and evidence under these conditions. The typical sampling efficiency we achieve for this problem is roughly 0 : 1% (i.e. 1000 likelihood calls per iteration), which translates to roughly 5 per dimension.", '18 J. S. Speagle': '<!-- image -->  \n<!-- image -->  \n<!-- image -->', 'Sampling terminates': "Figure 9. Illustration of dynesty 's performance sampling from a 200-D Gaussian using Hamiltonian Slice Sampling ( § 4.2.4) with gradients and no bounding distribution with only K = 50 live points. Top : Offsets in the recovered mean (left, black), variance (center, red), and covariance cross-terms (right, blue) relative to an expected mean of GLYPH<22> = 0 and covariance of C = ' 1 GLYPH<157> 2 ' I . Bottom : The estimated evidence ˆ Z (red line) along with the 1, 2, and 3-sigma errors (shaded). The true value is shown in black, along with the location where sampling terminates (dotted red vertical line).  \n<!-- image -->", '5.5 Comparison to MCMC': "Nested Sampling and MCMC sampling are different tools designed for different types of problems. Here we perform a limited comparison to highlight the advantages/disadvantages of each methodology.  \nWe consider a simple linear regression problem where our model is  \ny ' x ' = mx + b  \n(46)  \nand we observe noisy data from  \nˆ y i GLYPH<24> N GLYPH<16> y ' x i ' ; GLYPH<27> 2 i + » f y ' x i '… 2 GLYPH<17> (47)  \nwhere GLYPH<27> 2 i is the measured variance and f corresponds to an additional fractional systematic uncertainty that we would like to infer in addition to m and b . The likelihood is again Gaussian:  \nln L' m ; b ; f ' = GLYPH<0> 1 2 N GLYPH<213> i = 1 ln h 2 GLYPH<25> ' GLYPH<27> 2 i + f 2 ' mx i + b ' 2 ' i + » ˆ y i GLYPH<0> ' mx i + b '… 2 GLYPH<27> 2 + f 2 ' mx i + b ' 2 (48)  \nThis problem is unimodal and only has three parameters, making it very tractable to both Nested Sampling and MCMC methods.  \nWe choose our priors to be uniform so that m 2 »GLYPH<0> 5 ; 0 : 5 ' , b 2 » 0 ; 10 ' , and ln f 2 »GLYPH<0> 10 ; 1 … , which are substantially  \nbroader than the likelihood distribution but not so broad that the runtime of dynesty will be dominated merely integrating over the prior.  \nWe run dynesty in three configurations to sample from this posterior distribution, using the default settings whenever possible to highlight performance in a'typical'use case. First, we set the weight function to give the posterior 100% of the importance when allocating live points in order to imitate MCMC-like behavior. Then, we revert to the default 80%/20% posterior/evidence weighting scheme to see how much our posterior estimate degrades as we spend a larger fraction of runtime trying to improve our evidence estimates. Finally, we switch out the default sampling mode (uniform sampling) for random walks to forcibly decrease the overall sampling efficiency.  \nWe compare these results to two MCMC alternatives. The first is emcee (Foreman-Mackey et al. 2013), which is a common MCMC sampler used in astronomical analyses today. We opt to run it in its default configuration, which uses the 'stretch move' from (Goodman & Weare 2010) to make proposals, with K = 50 walkers. We initialize the walkers around the maximum-a-posteriori (MAP) solution based on the estimated covariance. We remove the first 300 samples from the chain to account for burn-in but do not count these 'wasted' samples when computing the overall sampling efficiency.  \nThe second alternative is a standard MH MCMC sam-  \n<!-- image -->  \nFigure 10. Comparison between dynesty and common MCMC alternatives inferring the slope m , intercept b , and (log-)fractional uncertainty ln f in a simple linear regression problem. See § 5.5 for additional details. Left : A corner plot showing the 1-D and 2-D marginalized posteriors for the slope m , intercept b , and (log-)fractional uncertainty ln f , with their true values in red. The 2.5%, 50%, and 97.5% percentiles (i.e. the 2-sigma credible region) are shown as vertical dashed lines. We see the posterior is well-constrained and roughly Gaussian. Right : The posterior sampling efficiency (i.e. the fraction of independent posterior samples generated per likelihood call) for dynesty , emcee , and simple MH MCMC plotted as a function of the total number of likelihood function calls. The predicted efficiency for a fixed effective sample size is shown in gray. We see that dynesty optimized for posterior estimation using uniform sampling (blue) can be up to 10x more efficient than emcee or MH MCMC at generating independent samples from the posterior. As expected, decreasing the emphasis on posterior vs evidence estimation to 80% (red) or using a less efficient but more flexible sampling method such as random walks (right) also reduces the overall efficiency.  \n<!-- image -->  \npler with a Gaussian proposal distribution. We take the covariance to be the same as that of the posterior distribution determined from the final set of weighted dynesty samples to create a relatively optimal proposal distribution. We then run with an identical setup to emcee (i.e. K = 50 chains initialized around MAP solution) to maintain consistency between approaches.  \nThe metric we use to compare between methods is the overall 'sampling efficiency', which we define to be the ratio of the estimated effective sample size (ESS) N ESS relative to the number of likelihood calls N call :  \nf samp GLYPH<17> N ESS N call (49)  \nFor dynesty , since the samples are all independent but assigned varying importance weights, we choose to estimate the ESS by counting the number of unique samples after using systematic resampling to redraw a set up equallyweighted samples. 6  \nFor the MCMC approaches, we use the standard definition of ESS as  \nN ESS = N (50)  \nGLYPH<28>  \nwhere GLYPH<28> is the auto-correlation averaged over all the chains. Since GLYPH<28> is computed for each parameter, to be conservative we set the value used to compute the ESS to be the maximum value. These choices tend to decrease the ESS by GLYPH<24> 25% relative to more optimistic ones but does not affect our overall conclusions.  \nWe compare the five different cases above and summarize the results from 25 independent trials in Figure 10. In all cases, we try to generate enough samples to give similar ESS between each approach based on dynesty 's default stopping criterion, which gives N ESS GLYPH<24> 17000 . We see that dynesty with uniform sampling within multiple bounding ellipsoids is roughly an order of magnitude more efficient at generating independent samples in this problem than MH MCMC and emcee . dynesty using random walks (i.e. running MH MCMC internally) gives efficiencies that are much more comparable to the two MCMC implementations.  \nAs discussed earlier, all methods experience some amount of overhead transitioning from the prior-dominated to posterior-dominated region. While this leads to glyph[lessorsimilar] 5% of samples being discarded for burn-in for the MCMC cases, it  \nFigure 11. Galaxy SED for object AEGIS 17 from the 3D-HST survey modeled with Prospector using dynesty . Left : A corner plot showing the 1-D and 2-D marginalized posteriors for the 14-parameter galaxy model. The 2.5%, 50%, and 97.5% percentiles (i.e. the 2-sigma credible region) are shown as vertical dashed lines. The posterior includes a bi-modal solution for the gas-phase metallicity. Top right : The modeled galaxy SED marginalized over the posterior. The 1-sigma (16-84% credible region) is also shown, along with the error-normalized residuals. The underlying model provides a reasonable fit to the observed data. Right middle : The median reconstructed star formation history as a function of look-back time along with the associated 16-84% credible region.  \n<!-- image -->  \nleads to a reduction in the ESS of GLYPH<24> 25% for dynesty . The fact that dynesty performs well even in this case illustrates how important Dynamic Nested Sampling is for ensuring samples are efficiently allocated during runtime.  \nThis result highlights the basic argument first outlined in § 2, illustrating that using Nested Sampling to sample from many simpler distributions in turn can sometimes be more effective than trying to sample from the posterior distribution directly with MCMC. In general, Nested Sampling performs well in cases like these where the likelihood varies smoothly in a given region and the prior has reasonable  \nbounds. In other cases where the prior is large or fewer samples from the posterior are needed, MCMC methods are more than sufficient.", '6 APPLICATIONS': "In addition to the toy problems in § 5, dynesty has also been applied in several packages and ongoing studies and shown to perform well when applied to real astronomical analyses. These include applications analyzing gravitational waves  \nFigure 12. Line-of-sight dust extinction (reddening) model for a sight-line in the Chameleon molecular cloud estimated with dynesty . Left : A corner plot showing the 1-D and 2-D marginalized posteriors for the 6-parameter line-of-sight model. The 16%, 50%, and 84% percentiles (i.e. the 1-sigma credible region) are shown as vertical dashed lines. The posterior includes a bi-modal solution for the cloud distance GLYPH<22> C as well as an extended tail for the foreground dust reddening f . Top right : The line-of-sight model from the estimated posterior. Individual distance-extinction posteriors for stars used in the fit as shown in grayscale, with most probable distance and extinction shown as a red cross. The blue line shows the typical extinction profile inferred for the sightline. The range of distance estimates is shown as the inverted blue histogram at the top of each panel, with the median cloud distance marked via the vertical blue line and yellow arrow and the 16-84% credible ranges marked via the vertical blue dashed lines. The horizontal blue lines show the estimated 1-sigma scatter in extinction behind the cloud.  \n<!-- image -->  \n(Ashton et al. 2018), exoplanets (Diamond-Lowe et al. 2018; Espinoza et al. 2018; Gunther et al. 2019), transients (Guillochon et al. 2018), galaxies (Leja et al. 2018a,b), and 3-D dust mapping (Zucker et al. 2018, 2019). We highlight two of these applications below that the author has been personally involved in.  \nIn Leja et al. (2018b), the authors modeled roughly 60k galaxy spectral energy distributions (SEDs) from the 3DHST survey (Brammer et al. 2012) over a redshift range of 0 : 5 < z < 2 : 5 . To conduct this analysis, they used the  \nBayesian SED fitting code Prospector (Johnson et al. in prep.), utilizing dynesty as their primary sampler, to sample from a 14-parameter model involving stellar mass, a non-parametric star formation history, stellar and gas metallicites, dust properties, and contributions from possible Active Galactic Nuclei. Compared to previous studies where emcee had been used to sample from the posterior (Leja et al. 2017, 2018c), the authors found that dynesty provided over an order of magnitude more efficient sampling and was  \nable to characterize a wide variety of posteriors. The results for a typical galaxy are shown in Figure 11.  \nIn Zucker et al. (2019), the authors used a combination of distance and reddening estimates to nearby stars from SED modeling (Speagle et al. in prep.) and Gaia parallax measurements (Gaia Collaboration et al. 2018) to derive distances to dozens of local molecular clouds. The distances to these clouds, however, are sensitive to the number and distribution of stars immediately in front of them as these stars help constrain the location of the 'jump' in dust extinction associated with the cloud. In cases where there are only a small number of foreground stars, this constraint can be quite weak, leading to extended posteriors with multi-modal solutions. This, along with the overall performance illustrated in Figure 10, motivated the use of dynesty to sample from the 6-parameter cloud distance model used in the analysis. We highlight one such multi-modal case in Chameleon in Figure 12.  \nThese examples, along with others listed earlier, are large-scale professional applications of dynesty that illustrate dynesty can work well in theory and in practice.", '7 CONCLUSION': "With Bayesian inference techniques now a large part of modern astronomical analyses, it has become increasingly important to develop and provide tools to the community that can help to 'bridge the gap' between writing the underlying model and estimating the corresponding posterior P' GLYPH<2> ' . Tools such as emcee , MultiNest , and PolyChord , which provide Markov Chain Monte Carlo and Nested Sampling implementations, have been heavily used and highly cited.  \nIn this paper we presented an overview of dynesty , a public, open-source, Python package that implements Dynamic Nested Sampling to enable flexible Bayesian inference over complex, multi-modal distributions. Building on previous work in the literature, we described the basics behind the Dyamic Nested Sampling approaches employed in the code, how we implement them, and how we use a variety of bounding and sampling methods to enable efficient inference. We then showcased dynesty 's performance on several toy problems as well as real astronomical application, highlighting its ability to estimate challenging posterior distributions both in theory and in practice.  \nWhile we have shown dynesty can perform similarly or better than existing MCMC approaches in one simple case, the real test for any package is based on users applying it to their analysis problems. We hope that dynesty will prove useful to the community and help facilitate exciting new science over the coming years.", 'ACKNOWLEDGEMENTS': "JSS is grateful to Rebecca Bleich for her support and patience.  \nThis project is the culmination of many individual efforts, not all of whom can be thanked here. First and foremost, JSS would like to thank Daniel Eisenstein, Charlie Conroy, and Doug Finkbeiner for their patience while he pursued this project, Catherine Zucker for her constant  \nstream of feedback during development, and Johannes Buchner for incredibly insightful and inspiring conversations. JSS would also like to thank Johannes Buchner, Hannah Diamond-Lowe, Daniel Eisenstein, Daniel Foreman-Mackey, Will Handley, Ben Johnson, Joel Leja, Locke Patton, and Catherine Zucker for feedback on early drafts that substantially improved the quality of this work. JSS would further like to thank Johannes Buchner, Phil Cargile, Ben Cook, James Guillochon, and Ben Johnson for their direct and indirect contributions to the dynesty codebase, as well as Kyle Barbary and collaborators for their contributions to nestle (upon which dynesty was initially based). JSS is also grateful to many beta-testers who provided invaluable feedback during dynesty 's development and suffered through many bugfixes, including (but not limited to) Gregory Ashton, Ana Bonaca, Phil Cargile, Tansu Daylan, Hannah DiamondLowe, Philipp Eller, Jonathan Fraine, Maximilian Gunther, Daniela Huppenkothen, Joel Leja, Sandro Tacchella, Ashley Villar, Catherine Zucker, and Joe Zuntz.  \nThis work has benefited from several software packages including numpy (van der Walt et al. 2011), scipy (Oliphant 2007), matplotlib (Hunter 2007), and corner (ForemanMackey 2016).", 'REFERENCES': "Ashton G., et al., 2018, arXiv e-prints, p. arXiv:1811.02042 Betancourt M., 2017, arXiv e-prints, p. arXiv:1701.02434  \nBlei D. M., Kucukelbir A., McAuliffe J. D., 2016, arXiv e-prints, Blitzstein J., Hwang J., 2014, Introduction to Probability. Chapman & Hall/CRC Texts in Statistical Science, CRC Press/Taylor & Francis Group, https://books.google.com/  \nbooks?id=ZwSlMAEACAAJ  \nBorne K., et al., 2009, in astro2010: The Astronomy and Astrophysics Decadal Survey. p. P6 ( arXiv:0909.3892 )  \n- Brammer G. B., et al., 2012, The Astrophysical Journal Supplement Series, 200, 13\n- Brewer B. J., P'artay L. B., Cs'anyi G., 2009, arXiv e-prints, p. arXiv:0912.2380\n- Brooks S., Gelman A., Jones G., Meng X.-L., 2011, Handbook of Markov Chain Monte Carlo. CRC press\n- Buchner J., 2016, Statistics and Computing, 26, 383\n- Buchner J., 2017, arXiv e-prints, p. arXiv:1707.04476\n- Carpenter B., et al., 2017, Journal of Statistical Software, Articles, 76, 1\n- Chopin N., Ridgway J., 2015, preprint, ( arXiv:1506.08640 )\n- Chopin N., Robert C. P., 2010, Biometrika, 97, 741\n- Diamond-Lowe H., Berta-Thompson Z., Charbonneau D., Kempton E. M. R., 2018, AJ, 156, 42\n- Douc R., Capp'e O., Moulines E., 2005, arXiv e-prints, p. cs/0507025\n- Efron B., 1979, Ann. Statist., 7, 1\n- Espinoza N., Kossakowski D., Brahm R., 2018, arXiv e-prints, p. arXiv:1812.08549  \nFeigelson E. D., 2017, in Brescia M., Djorgovski S. G., Feigelson E. D., Longo G., Cavuoti S., eds, IAU Symposium Vol. 325, Astroinformatics. pp 3-9 ( arXiv:1612.06238 ), doi:10.1017/S1743921317003453  \nFeroz F., Hobson M. P., 2008, MNRAS, 384, 449  \n- Feroz F., Skilling J., 2013, in von Toussaint U., ed., American Institute of Physics Conference Series Vol. 1553, American Institute of Physics Conference Series. pp 106-113 ( arXiv:1312.5638 ), doi:10.1063/1.4819989  \nFeroz F., Hobson M. P., Bridges M., 2009, MNRAS, 398, 1601  \nFeroz F., Hobson M. P., Cameron E., Pettitt A. N., 2013, preprint, ( arXiv:1306.2144 )  \nFisher R. A., 1922, Philosophical Transactions of the Royal Society of London Series A, 222, 309  \nForeman-Mackey D., 2016, The Journal of Open Source Software, 24  \nForeman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306  \nGaia Collaboration et al., 2018, A&A, 616, A1  \nGelman A., Rubin D. B., 1992, Statistical Science, 7, 457  \nGeman S., Geman D., 1987, in Fischler M. A., Firschein O., eds, , Readings in Computer Vision. Morgan Kaufmann, San Francisco (CA), pp 564 - 584, doi:https://doi.org/10.1016/B9780-08-051581-6.50057-X, http://www.sciencedirect.com/ science/article/pii/B978008051581650057X  \nGoodman J., Weare J., 2010, Communications in Applied Mathematics and Computer Science, 5, 65  \nGuillochon J., Nicholl M., Villar V. A., Mockler B., Narayan G., Mandel K. S., Berger E., Williams P. K. G., 2018, The Astrophysical Journal Supplement Series, 236, 6  \nGunther M. N., et al., 2019, arXiv e-prints, p. arXiv:1903.06107 Handley W. J., Hobson M. P., Lasenby A. N., 2015, MNRAS, 450, L61  \nHastings W., 1970, Biometrika, 57, 97  \nHeavens A., Fantaye Y., Mootoovaloo A., Eggers H., Hosenie Z., Kroon S., Sellentin E., 2017, arXiv e-prints, p. arXiv:1704.03472  \nHigson E., Handley W., Hobson M., Lasenby A., 2017a, arXiv e-prints, p. arXiv:1703.09701  \nHigson E., Handley W., Hobson M., Lasenby A., 2017b, arXiv e-prints, p. arXiv:1704.03459  \n- Higson E., Handley W., Hobson M., Lasenby A., 2019, MNRAS, 483, 2044\n- Hoffman M. D., Gelman A., 2011, arXiv e-prints, p. arXiv:1111.4246  \nHol J. D., Schon T. B., Gustafsson F., 2006, in 2006 IEEE Nonlinear Statistical Signal Processing Workshop. pp 79-82, doi:10.1109/NSSPW.2006.4378824  \nHunter J. D., 2007, Computing in Science Engineering, 9, 90 Jasa T., Xiang N., 2012, Acoustical Society of America Journal, 132, 3251  \nKeeton C. R., 2011, MNRAS, 414, 1418  \nLartillot N., Philippe H., 2006, Systematic Biology, 55, 195  \n- Leja J., Johnson B. D., Conroy C., van Dokkum P. G., Byler N., 2017, ApJ, 837, 170\n- Leja J., Carnall A. C., Johnson B. D., Conroy C., Speagle J. S., 2018a, arXiv e-prints, p. arXiv:1811.03637  \nLeja J., et al., 2018b, arXiv e-prints, p. arXiv:1812.05608  \nLeja J., Johnson B. D., Conroy C., van Dokkum P., 2018c, ApJ, 854, 62  \nMetropolis N., Rosenbluth A. W., Rosenbluth M. N., Teller A. H., Teller E., 1953, J. Chem. Phys., 21, 1087  \nMukherjee P., Parkinson D., Liddle A. R., 2006, ApJ, 638, L51  \nNagaraja H. N., 2006, Order Statistics from Independent Exponential Random Variables and the Sum of the Top Order Statistics. Birkhauser Boston, Boston, MA, pp 173-185, doi:10.1007/0-8176-4487-3 11, https://doi.org/ 10.1007/0-8176-4487-3\\_11  \nNeal R. M., 2003, Ann. Statist., 31, 705  \nNeal R. M., 2012, arXiv e-prints, p. arXiv:1206.1901  \nOliphant T. E., 2007, Computing in Science Engineering, 9, 10 Planck Collaboration et al., 2016, A&A, 594, A20  \nPlummer M., 2003, in Proceedings of the 3rd International Workshop on Distributed Statistical Computing.  \nSalomone R., South L. F., Drovandi C. C., Kroese D. P., 2018, arXiv e-prints, p. arXiv:1805.03924  \nSharma S., 2017, Annual Review of Astronomy and Astrophysics,  \n55, 213  \nShaw J. R., Bridges M., Hobson M. P., 2007, MNRAS, 378, 1365 Sivia D., Skilling J., 2006, Data analysis: a Bayesian tutorial. Oxford science publications, Oxford University Press, https: //books.google.com/books?id=6O8ZAQAAIAAJ  \nSkilling J., 2004, in Fischer R., Preuss R., Toussaint U. V., eds, American Institute of Physics Conference Series Vol. 735, American Institute of Physics Conference Series. pp 395-405, doi:10.1063/1.1835238  \nSkilling J., 2006, Bayesian Anal., 1, 833  \nSkilling J., 2012, in Goyal P., Giffin A., Knuth K. H., Vrscay E., eds, American Institute of Physics Conference Series Vol. 1443, American Institute of Physics Conference Series. pp 145-156, doi:10.1063/1.3703630  \nTrotta R., 2008, Contemporary Physics, 49, 71  \nVehtari A., Gelman A., Simpson D., Carpenter B., Burkner P.-C., 2019, arXiv e-prints, p. arXiv:1903.08008  \nYork D. G., et al., 2000, AJ, 120, 1579  \nZucker C., Schlafly E. F., Speagle J. S., Green G. M., Portillo S. K. N., Finkbeiner D. P., Goodman A. A., 2018, ApJ, 869, 83 Zucker C., Speagle J. S., Schlafly E. F., Green G. M., Finkbeiner D. P., Goodman A. A., Alves J., 2019, arXiv e-prints, p. arXiv:1902.01425  \nvan der Walt S., Colbert S. C., Varoquaux G., 2011, Computing in Science Engineering, 13, 22", 'APPENDIX A: DETAILED NESTED SAMPLING RESULTS': 'While we presented a broad overview of Nested Sampling in the main text, we glossed over much of the statistical background. We include more detailed results and discussion below.  \nThe outline of these results are as follows. In § A1 we outline the basic setup for Nested Sampling. In § A2 we derive statistical properties in the single live point case. In § A3 we discuss the process of utilizing multiple live points. In § A4 we derive properties in the many live point case. In § A5 we extend these results to encompass varying numbers of live points. Finally, in § A6 we discuss various error properties of Nested Sampling as well as schemes to estimate them.', 'A1 Setup': 'Following Skilling (2006), Feroz et al. (2013), and others, we start by (re-)defining Bayes Rule  \nP\' GLYPH<2> \' = L\' GLYPH<2> \' GLYPH<25> \' GLYPH<2> \' Z (A1)  \nwhere P\' GLYPH<2> \' is the posterior, L\' GLYPH<2> \' is the likelihood, GLYPH<25> \' GLYPH<2> \' is the prior, and  \nZ M = \' GLYPH<10> GLYPH<2> L\' GLYPH<2> \' GLYPH<25> \' GLYPH<2> \' d GLYPH<2> (A2)  \nis the evidence.  \nTo evaluate this integral, Nested Sampling seeks to transform it from one over position GLYPH<2> to one over prior volume X where  \nX \' GLYPH<21> \' GLYPH<17> \' GLYPH<10> GLYPH<2> s : t : L\' GLYPH<2> \' >GLYPH<21> GLYPH<25> \' GLYPH<2> \' d GLYPH<2> GLYPH<17> \' GLYPH<10> GLYPH<2> GLYPH<25> GLYPH<21> \' GLYPH<2> \' d GLYPH<2> (A3)  \ndefines the prior volume within a given iso-likelihood contour of level GLYPH<21> , assuming our priors are integrable, and  \nGLYPH<25> GLYPH<21> \' GLYPH<2> \' GLYPH<17> ( GLYPH<25> \' GLYPH<2> \'GLYPH<157> X \' GLYPH<21> \' L\' GLYPH<2> \' GLYPH<21> GLYPH<21> 0 L\' GLYPH<2> \' < GLYPH<21> (A4)  \nis the constrained prior. Note that X 2 \' 0 ; 1 … since the integral over the entire prior is x \' GLYPH<21> = 0 \' = 1 while the value as GLYPH<21> !1 should approach 0 if the maximum-likelihood value L max is a singular point.  \nSince GLYPH<21> 2 » 0 ; 1\' , this allows us to redefine the evidence integral as  \nZ = \' 1 0 X \' GLYPH<21> \' d GLYPH<21> (A5)  \nProvided the inverse L\' X \' of X \'L\' GLYPH<2> \' = GLYPH<21> \' exists (i.e. there are no flat \'slabs\' of likelihood anywhere, only contours), we can rewrite this integral in terms of the prior volume associated with a particular iso-likelihood contour:  \nZ = \' 1 0 L\' X \' d X (A6)  \nThis is now a 1-D integral over X that we can approximate using a discrete set of N points using, e.g., a Riemann sum  \nˆ Z = N GLYPH<213> i = 1 L\' GLYPH<2> i \' GLYPH<2> \' X i GLYPH<0> X i GLYPH<0> 1 \' GLYPH<17> N GLYPH<213> i = 1 p \' GLYPH<2> i \' (A7)  \nwhere X 0 = 1 and p \' GLYPH<2> i \' is the (un-normalized) importance weight. These values can also be used to approximate the posterior:  \nˆ P\' GLYPH<2> \' = " N i = 1 p \' GLYPH<2> i \' GLYPH<14> \' GLYPH<2> i \' " N i = 1 p \' GLYPH<2> i \' (A8)', 'A2 Using a Single Live Point': "Unfortunately, the exact value of X ' GLYPH<21> ' at a given likelihood level GLYPH<21> = L' GLYPH<2> ' is unknown. We can, however, construct an estimator ˆ X with a known statistical distribution. Looking back at the definition of the prior volume X 'L' , we see that it defines a cumulative distribution function (CDF) over L . We can then define the associated probability density function (PDF) for L as  \nP 'L' GLYPH<17> d X 'L' d L = d d L ' GLYPH<10> GLYPH<2> GLYPH<25> L ' GLYPH<2> ' d GLYPH<2> (A9)  \nAssuming we can sample L from its PDF P 'L' , we can use the Probability Integral Transform (PIT) to subsequently constrain the distribution of X 'L' . In other words:  \nL 0 GLYPH<24> P 'L' ) X 'L 0 ' GLYPH<24> Unif (A10)  \nwhere X GLYPH<24> f ' X ' notation implies the random variable X is drawn from f ' X ' and Unif is the standard Uniform distribution (i.e. flat from 0 to 1 ). This can be directly extended to cases where we are interested in sampling relative to a given threshold GLYPH<21> as  \nL 0 GLYPH<24> P 'LjL > GLYPH<21> ' ) X 'L 0 ' X ' GLYPH<21> ' GLYPH<24> Unif (A11)  \nWhile this does not appear to make things any easier, it actually helps us out enormously. That's because, at fixed GLYPH<21> , X ' GLYPH<21> ' is actually a CDF over the constrained prior GLYPH<25> GLYPH<21> ' GLYPH<2> ' . That means we can bypass GLYPH<21> and P 'L' altogether and just sample from GLYPH<25> GLYPH<21> ' GLYPH<2> ' directly to satisfy the PIT:  \nGLYPH<2> 0 GLYPH<24> GLYPH<25> GLYPH<21> ' GLYPH<2> ' ) X 'L' GLYPH<2> 0 '' X ' GLYPH<21> ' GLYPH<24> Unif (A12)  \nVarious methods for sampling from the constrained prior GLYPH<25> GLYPH<21> ' GLYPH<2> ' subject to a suitable prior transform T (see § 2.2) are outlined in § 4.  \nBefore moving on, we want to quickly note that while the above scheme is sufficient for generating values of L 0 GLYPH<24> P 'L' it is by no means necessary . As a counter-example, we can imagine a function f ' t ' ! GLYPH<2> t that traces out a singular path through the distribution with support over L' f ' 2 »L min ; L max … . Let us furthermore assume that we construct f ' t ' such that we spend more 'time' t where the likelihood PDF is higher so that the PDF P ' t ' / P 'L' GLYPH<2> t '' . Finally, let's define the constrained function f GLYPH<21> ' t ' to simply be the portion of the path with L' GLYPH<2> t ' > GLYPH<21> . While this path by no means encompasses the prior, it is clear that  \nt 0 GLYPH<24> f GLYPH<21> ' t ' ) X 'L' GLYPH<2> t 0 '' X ' GLYPH<21> ' GLYPH<24> Unif (A13)  \nThis result proves we can in theory satisfy the PIT for Nested Sampling using correlated samples provided they probe enough of the local portion of the prior to obtain sufficient coverage over the range of possible likelihoods (see also Salomone et al. 2018). It also provides support for why Nested Sampling works so well in practice even when samples are not fully independent.  \nFor a given prior volume X i GLYPH<0> 1 associated with a given likelihood level GLYPH<21> i GLYPH<0> 1 = L' GLYPH<2> i GLYPH<0> 1 ' after i GLYPH<0> 1 iterations of this procedure, this implies the current prior volume X i will be  \nˆ X i = U i ˆ X i GLYPH<0> 1 = i GLYPH<214> j = 1 U j (A14)  \nwhere  \nU 1 ; : : : ; U i iid GLYPH<24> Unif (A15)  \nare independent and identically distributed (iid) random variables drawn from the standard Uniform distribution and we have taken X 0 GLYPH<17> X ' GLYPH<21> = 0 ' = 1 . As we do not actually know the values of U 1 ; : : : ; U j , we consider ˆ X i to be a noisy estimator of X i .  \nWhile sampling, we obviously need to assign a value for ˆ X i to determine, e.g., whether to stop. While we can easily simulate random values of U 1 ; : : : ; U i , if we want these values to be consistent then a reasonable choice is the expectation value (arithmetic mean):  \nE GLYPH<2> ˆ X i GLYPH<3> = E 2 6 6 6 6 i GLYPH<214> j = 1 U j 3 7 7 7 7 = i GLYPH<214> j = 1 E GLYPH<2> U j GLYPH<3> = GLYPH<18> 1 2 GLYPH<19> i (A16)  \n4  \n5  \nAlternately, we might also be interested in the expectation value of ln ˆ X i (geometric mean):  \nE GLYPH<2> ln ˆ X i GLYPH<3> = i GLYPH<213> j = 1 E GLYPH<2> ln U j GLYPH<3> GLYPH<24> GLYPH<0> i GLYPH<213> j = 1 E GLYPH<2> E j GLYPH<3> = GLYPH<0> i (A17)  \nwhere we have used the fact that  \nU GLYPH<24> Unif )GLYPH<0> ln U GLYPH<24> Expo  \n(A18)  \nwhere Expo is the standard Exponential distribution and  \nE 1 ; : : : ; E i iid GLYPH<24> Expo (A19)  \nVarious stopping criteria are discussed in the main text ( § 2.4 and § 3.4) and so are not discussed further here.", 'A3 Combining Live Points': "Following Higson et al. (2017a), let's consider the case where we have two independent live points following the basic sampling approach described above. These each form a set of samples with increasing likelihood  \nL » 1 … N 1 > GLYPH<1> GLYPH<1> GLYPH<1> > L » 1 … 1 > 0 L » 2 … N 2 > GLYPH<1> GLYPH<1> GLYPH<1> > L » 2 … 1 > 0  \nwhere the »GLYPH<1>… superscript notation indicates the index of the associated live point. We now want to 'merge' these two sets of ordered samples together to get a single hypothetical ordered list:  \nL » 1 … N 1 > L » 2 … N 2 > GLYPH<1> GLYPH<1> GLYPH<1> > L » 1 … 2 > L » 2 … 2 > L » 2 … 1 > L » 1 … 1 > 0 !L N > GLYPH<1> GLYPH<1> GLYPH<1> > L 1 > 0  \nwhere N = N 1 + N 2 .  \nIndependently, we know that the prior volume at a given iteration for each live point is just  \nˆ X » j … i = i GLYPH<214> n = 1 U » j … n  \nWhat we want to know, however, is the distribution of ˆ X i of the merged list. Considering each sample independently implies X 2 = X » 2 … 1 and X 1 = X » 1 … 1 follow the same distribution (i.e. the first sampled prior volume for each run is similarly distributed). However, considering them together (based on the merged list) implies X 2 = X » 2 … 1 is strictly less than X 1 = X » 1 … 1 since L 2 > L 1 . This tells us that ˆ X i cannot follow the same distribution from the associated independent runs that comprise it.  \nWith this finding in hand, we now consider an approach for sampling from the prior volume using two live points. At each iteration i , we remove the one with the lowest likelihood GLYPH<21> = L min i and replace it with a new point sampled from the constrained prior GLYPH<25> GLYPH<21> ' GLYPH<2> ' . After N iterations, we will end up with a sorted list of likelihoods L N > GLYPH<1> GLYPH<1> GLYPH<1> > L 1 > 0 . If, however, we look at each live point individually (i.e. ignoring the other live point), we would find that each live point's evolution would comprise a list of independent samples with ordered likelihoods that would each be identical to L » 1 … N 1 > GLYPH<1> GLYPH<1> GLYPH<1> > L » 1 … 1 > 0 and L » 2 … N 2 > GLYPH<1> GLYPH<1> GLYPH<1> > L » 2 … 1 > 0 , respectively! Therefore, we see that this procedure for sampling with two live points is identical to combining two sets of independent samples derived using one live point each.  \nThe above procedure can be immediately generalized to K live points, producing the (Static) Nested Sampling procedure outlined in Algorithm 1. We will return to this duality between K independent Nested Sampling runs and a single Nested Sampling run with K live points in § A6.", 'A4 Using Many Live Points': 'Now that we have established a procedure for running Nested Sampling with K live points, we need to characterize how this affects our estimates ˆ X i of the prior volume. At any given iteration i , we know that the current set of prior volumes f X » 1 … i ; : : : ; X » K … i g associated with our K live points  \nare uniformly distributed within the prior volume from the previous iteration X i GLYPH<0> 1 so that  \nX » j … i = U » j … X i GLYPH<0> 1 (A20)  \nwhere  \nU » 1 … ; : : : ; U » K … iid GLYPH<24> Unif (A21)  \nWe are now want to replace the live point with the lowest likelihood L min i corresponding to the largest prior volume. This means we are now interested in the ordered list of prior volumes  \nX \' j \' i = U \' j \' X i GLYPH<0> 1 (A22)  \nwhere \' j \' now indicates the position in the ordered list (from smallest to largest) rather than the live point index » j … and  \nU \' j \' = min j GLYPH<16> n U » 1 … ; : : : ; U » K … oGLYPH<17> (A23)  \nis the j th standard uniform order statistic, where min j selects the j th smallest point (so j = 1 is the smallest and j = K is the largest).  \nUsing the Reny\'ı Representation, it can be shown that we can represent the joint distribution of our K standard uniform order statistics f U \' 1 \' ; : : : ; U \' K \' g such that (Nagaraja 2006):  \nU \' j \' = " j n = 1 E n " K + 1 n = 1 E n (A24)  \nwhere  \nE 1 ; : : : ; E K + 1 iid GLYPH<24> Expo (A25)  \nThe marginal distribution for U \' j \' is then (Blitzstein & Hwang 2014):  \nU \' j \' GLYPH<24> Beta \' j ; K + 1 GLYPH<0> j \' (A26)  \nwhere Beta \' GLYPH<11>; GLYPH<12> \' is the Beta distribution.  \nUsing these results, we see that the prior volume based on K live points at iteration i evolves as  \nˆ X i = i GLYPH<214> j = 1 U \' K \' j (A27)  \nwhere U \' K \' 1 ; : : : ; U \' K \' i are iid draws of the K th standard uniform order statistic with marginal distribution Beta \' K ; 1 \' . The arithmetic mean is  \nE GLYPH<2> ˆ X i GLYPH<3> = i GLYPH<214> j = 1 E h U \' K \' j i = GLYPH<18> K K + 1 GLYPH<19> i (A28)  \nThe geometric mean is  \nE GLYPH<2> ln ˆ X i GLYPH<3> = i GLYPH<213> j = 1 E h ln U \' K \' j i = GLYPH<0> i K (A29)  \nAs discussed in § 2.3, after we terminate sampling we can add the final set of K live points to our set of N samples. These will then just follow the final set of f U \' 1 \' ; : : : ; U \' K \' g standard uniform order statistics relative to ˆ X N with an arithmetic mean  \nE GLYPH<2> ˆ X N + k GLYPH<3> = GLYPH<18> K + 1 GLYPH<0> k K + 1 GLYPH<19> GLYPH<18> K K + 1 GLYPH<19> N (A30)  \nand geometric mean  \nE GLYPH<2> ln ˆ X N + k GLYPH<3> = GLYPH<0> N K GLYPH<0> » \' K + 1 \' GLYPH<0> \' K + 1 GLYPH<0> k \'… (A31)  \nwhere \'GLYPH<1>\' is the digamma function.', 'A5 Using a Varying Number of Live Points': 'As discussed in § 3, there\'s no inherent reason why the number of number of live points must remain constant from iteration to iteration. Indeed, we can interpret adding the final set of live points to the list of samples from § A4 as simply allowing the nested sampling run to continue while continually decreasing the number of live points. From this viewpoint, we have K 1 = GLYPH<1> GLYPH<1> GLYPH<1> = K N = K live points over iteration i = 1 to N , but only K N + k = K + 1 GLYPH<0> k live points at iteration i = N + k .  \nThe change in the number of live points also changes the overall behavior of the Nested Sampling run before and after adding the final set of live points. We can highlight these by rewriting the results from § A4 as:  \nln E GLYPH<2> ˆ X N + k GLYPH<3> = N GLYPH<213> i = 1 ln GLYPH<18> K K + 1 GLYPH<19> | {z } Exponential Shrinkage + k GLYPH<213> j = 1 ln GLYPH<18> K + 1 GLYPH<0> k K + 2 GLYPH<0> k GLYPH<19> | {z } Uniform Shrinkage (A32)  \nThis neatly decomposes the two \'modes\' in which Nested Sampling can traverse the prior. While \'replacing\' the worst live point (i.e. K i = K i GLYPH<0> 1 ), the prior volume shrinks exponentially by a constant factor at each iteration. However, when \'removing\' live points (i.e. K i < K i GLYPH<0> 1 ), we instead shrink uniformly by a variable factor at each iteration.  \nWe can now generalize this behavior to the case where K i is allowed to vary at each iteration (Higson et al. 2017b). This now generates two distinct classes of behavior. When K i GLYPH<21> K i GLYPH<0> 1 , we add K i GLYPH<0> K i GLYPH<0> 1 GLYPH<21> 0 live points to our existing set of live points, after which we replace the one with the worst likelihood L min i . This then gives a distribution for the prior volume shrinkage of Beta \' K i ; 1 \' .  \nIf K i < K i GLYPH<0> 1 , on the other hand, we instead have removed K i GLYPH<0> 1 GLYPH<0> K i live points from the previous set of live points. The expected shrinkage is then based on the associated K i standard uniform order statistic U \' K i \' from the initial set of K i GLYPH<0> 1 values. Although in theory we should consider cases where the number of live points can decrease by an arbitrary amount, in practice when following iterative schemes such as the one outlined in Algorithm 3 we only need to consider the case where K i GLYPH<0> 1 GLYPH<0> K i = 1 .  \nTaken together, these two types of behavior then give a mean estimate of:  \nln E GLYPH<2> ˆ X j GLYPH<3> = n 1 GLYPH<213> i = 1 ln GLYPH<18> K i K i + 1 GLYPH<19> + n 2 GLYPH<213> i = 1 ln GLYPH<18> K N 1 + 1 GLYPH<0> i K N 1 + 2 GLYPH<0> i GLYPH<19> + n 3 GLYPH<213> i = 1 ln GLYPH<18> K N 2 + i K N 2 + i + 1 GLYPH<19> + GLYPH<1> GLYPH<1> GLYPH<1> + n M GLYPH<0> 1 GLYPH<213> i = 1 ln GLYPH<18> K N M GLYPH<0> 2 + i + 1 GLYPH<0> i K N M GLYPH<0> 2 + i + 1 GLYPH<19> (A33) + n M GLYPH<213> i = 1 ln GLYPH<18> K N M GLYPH<0> 1 + 1 GLYPH<0> i K N M GLYPH<0> 1 + 2 GLYPH<0> i GLYPH<19>  \nwhere n m is the number of contiguous samples for which either exponential or uniform shrinkage dominates, N m = " i k = 1 n k is the total number of iterations that have occurred  \nup to that point, and M is the number of contiguous regions prior to iteration j = N M where one mode of shrinkage dominates. Note that for illustrative purposes here we have assumed the final samples are experiencing uniform shrinkage.  \nTo summarize, varying the number of live points at each iteration simply involves dynamically switching between exponential and uniform shrinkage over the course of a Nested Sampling run. While this adds additional bookkeeping, it remains straightforward to estimate the prior volume ˆ X i at any particular iteration.', 'A6 Nested Sampling Errors': 'We now turn our attention to characterizing various error properties of Nested Sampling, following the basic approach of Higson et al. (2017a, 2019). Similar to other sampling approaches, we expect some amount of \'sampling noise\' in our evidence ˆ Z and posterior estimates ˆ P\' GLYPH<2> \' arising from the fact that we are approximating a continuous distribution (and smooth integral) with a discrete set of N samples. We expect that as the number of live points at each iteration K i !1 such that change in prior volume X i GLYPH<0> X i GLYPH<0> 1 ! 0 and the total number of samples N !1 , these sampling errors will become negligible.  \nUnlike other sampling approaches such as Markov Chain Monte Carlo (MCMC), however, Nested Sampling, contains an additional source of noise arising from our use of noisy estimators X ! ˆ X i of the prior volume at a given iteration i (Skilling 2006). This \'statistical noise\' translates to a noisy estimator of the importance weight p \' GLYPH<2> i \' ! ˆ p \' GLYPH<2> i \' , which in turn gives noisy estimators for our previous evidence estimate  \nˆ Z = N GLYPH<213> i = 1 L\' GLYPH<2> i \' GLYPH<2> \' X i GLYPH<0> X i GLYPH<0> 1 \' GLYPH<25> N GLYPH<213> i = 1 L\' GLYPH<2> i \' GLYPH<2> \' ˆ X i GLYPH<0> ˆ X i GLYPH<0> 1 \' GLYPH<17> N GLYPH<213> i = 1 ˆ p \' GLYPH<2> i \' (A34)  \nand our previous posterior estimate  \nˆ P\' GLYPH<2> \' = " N i = 1 p \' GLYPH<2> i \' GLYPH<14> \' GLYPH<2> i \' " N i = 1 p \' GLYPH<2> i \' GLYPH<25> " N i = 1 ˆ p \' GLYPH<2> i \' GLYPH<14> \' GLYPH<2> i \' " N i = 1 ˆ p \' GLYPH<2> i \' (A35)  \nSimilar with the sampling noise, we also expect the statistical noise to become negligible as the number of live points at each iteration K i ! 1 such that our estimate ˆ X i ! X i and the total number of samples N !1 .  \nWe can highlight the decomposition of these two noise sources by considering trying to evaluate the expectation value of a target function f \' GLYPH<2> \' with respect to the posterior (Chopin & Robert 2010; Higson et al. 2017a):  \nE P » f … = \' GLYPH<10> GLYPH<2> f \' GLYPH<2> \'P\' GLYPH<2> \' d GLYPH<2> = 1 Z \' 1 0 ˜ f \' X \'L\' X \' d X (A36)  \nwhere  \n˜ f \' X \' = E GLYPH<25> » f \' GLYPH<2> \'jL\' GLYPH<2> \' = L\' X \'… (A37)  \nis the expectation value of f \' GLYPH<2> \' on the associated isolikelihood contour L\' GLYPH<2> \' = L\' X \' with respect to the prior  \nGLYPH<25> \' GLYPH<2> \' . Using the same Riemann sum approximation as § A1, Nested Sampling would approximate this integral as:  \nE P » f … GLYPH<25> N GLYPH<213> i = 1 ˜ f \' X i \' L\' X i \'\' X i GLYPH<0> X i GLYPH<0> 1 \' Z = N GLYPH<213> i = 1 ˜ f \' X i \' p \' X i \' (A38)  \nGLYPH<25> N GLYPH<213> i = 1 f \' GLYPH<2> i \' p \' X i \' (A39)  \nGLYPH<25> N GLYPH<213> i = 1 f \' GLYPH<2> i \' ˆ p \' GLYPH<2> i \' (A40)  \nWe can see the two error types enter in cleanly through the final two approximations. In equation (A39), we introduce sampling noise by replacing ˜ f \' X \' , which is averaged over the entire iso-likelihood contour, with the estimate f \' GLYPH<2> i \' evaluated at a single point. Then, in equation (A40), we replace the true importance weight p \' X i \' at a given prior volume with its noisy estimate ˆ p \' GLYPH<2> i \' based on our noisy estimators for the prior volume ˆ X i .', 'A6.1 Statistical Uncertainties': 'In § A2, § A4, and § A5, we derived the analytic distribution for our prior volume estimator ˆ X i at iteration i under a variety of assumptions. While the distribution for ˆ Z and ˆ P\' GLYPH<2> \' is not analytic, it is straightforward to draw from them. First, we simulate values of the prior volumes  \nˆ X 0 1 ; : : : ; ˆ X 0 N GLYPH<24> P \' ˆ X 1 ; : : : ; ˆ X N \' (A41)  \nby drawing a combination of Beta \' K i ; 1 \' -distributed random variables (when K i GLYPH<21> K i GLYPH<0> 1 ) and standard uniform order statistics (when K i < K i GLYPH<0> 1 ) and iteratively computing each ˆ X 0 i using the procedures outlined earlier. Then, we simply compute the corresponding evidence ˆ Z 0 and posterior ˆ P 0 \' GLYPH<2> \' estimates.  \nWhile we can simulate the prior volumes and trace their impact on ˆ Z and ˆ P\' GLYPH<2> \' explicitly, it is also helpful to derive a rough estimate of their impact. Since the posterior P\' GLYPH<2> \' can be arbitrarily complex, we will focus on the evidence Z for which this analysis is more tractable.  \nThere have previously been two main approaches for deriving the uncertainty, which focus either on trying to derive V h ˆ Z i (Keeton 2011) or V h ln ˆ Z i (Skilling 2006). Here we will focus on the latter, which gives a cleaner (if less precise) result.  \nWe first start with the Static Nested Sampling case using a constant number of live points K . To estimate the evidence Z , we must integrate over the unnormalized posterior P\' GLYPH<2> \' / GLYPH<25> \' GLYPH<2> \'L\' GLYPH<2> \' . This occurs after a certain number of iterations N have passed given a fixed stopping criterion.  \nThere are two factors that contribute to the overall N . The first is the rate of integration: at any given iteration i , the prior volume decreases by GLYPH<1> ln X GLYPH<25> 1 GLYPH<157> K . As a result, it must be the case that N / 1 GLYPH<157> K .  \nThe second is the total amount of prior volume that needs to be integrated over. This roughly scales as the Kullback-Leibler (KL) divergence (i.e. \'information gain\')  \nbetween the prior GLYPH<25> \' GLYPH<2> \' and posterior P  \nH \'Pjj GLYPH<25> \' GLYPH<17> H GLYPH<17> \' GLYPH<10> GLYPH<2> P\' GLYPH<2> \' ln P\' GLYPH<2> \' GLYPH<25> \' GLYPH<2> \' d GLYPH<2> (A42)  \n= 1 Z \' 1 0 L\' X \' ln L\' X \' d X GLYPH<0> ln Z (A43)  \nSince N is a discrete number that is typically large, it is reasonable to assume that it follows a Poisson distribution such that  \nE » N … = V » N … GLYPH<24> H GLYPH<1> ln X (A44)  \nThis leads to a rough uncertainty in ln ˆ Z of  \nGLYPH<27> h ln ˆ Z i GLYPH<24> GLYPH<27> GLYPH<2> ln ˆ X N GLYPH<3> GLYPH<24> GLYPH<27> » ln N … \' GLYPH<1> ln X \' GLYPH<24> p H \' GLYPH<1> ln X \' = r H K (A45)  \nWe now extend this result to encompass a variable number of live points K i at each iteration. We first rewrite our estimate of the variance as  \nV h ln ˆ Z i = V " N GLYPH<213> i = 1 GLYPH<16> ln ˆ Z i GLYPH<0> ln ˆ Z i GLYPH<0> 1 GLYPH<17> # GLYPH<17> V " N GLYPH<213> i = 1 GLYPH<1> ln ˆ Z i # GLYPH<25> N GLYPH<213> i = 1 V h GLYPH<1> ln ˆ Z i i (A46)  \nwhere the final approximation assumes the distribution of evidence updates is independent at each iteration i and ln ˆ Z 0 = 0 . If we further assume that the distribution of the actual evidence estimates ˆ Z i themselves are roughly independent at each iteration i and that the number of live points K i changes sufficiently slowly such that GLYPH<1> ln X i GLYPH<25> GLYPH<1> ln X i GLYPH<0> 1 , we find  \nV h GLYPH<1> ln ˆ Z i i GLYPH<25> V h ln ˆ Z i i GLYPH<0> V h GLYPH<1> ln ˆ Z i GLYPH<0> 1 i GLYPH<24> \' H i GLYPH<0> H i GLYPH<0> 1 \'\' GLYPH<1> ln X i \' GLYPH<17> \' GLYPH<1> H i \'\' GLYPH<1> ln X i \' (A47)  \nSubstituting this in to our original expression and taking GLYPH<1> ln X i GLYPH<25> 1 GLYPH<157> K i then gives a modified error estimate  \nGLYPH<27> h ln ˆ Z i GLYPH<24> v u t N GLYPH<213> i = 1 GLYPH<1> H i K i (A48)  \nWhile the modified estimator in equation (A48) is less reliable that our original estimate, it is somewhat reassuring that in the special case K 1 = GLYPH<1> GLYPH<1> GLYPH<1> = K N = K it reduces to the original estimator derived in equation (A45).', 'A6.2 Sampling Uncertainties': 'Unlike the statistical uncertainties on the prior volume estimator ˆ X i , we do not have analytic expression or ways to explicitly simulate from the distribution characterized by our sampling uncertainties. In fact, it is doubtful we will ever have access to these except in special cases, since they rely on having access to the distribution of all possible paths (or varying lengths) live points can take through the distribution over the course of a Nested Sampling run.  \nThis however, does not mean we cannot attempt to construct an estimate of this distribution. To do this, we follow Higson et al. (2017a) and turn to bootstrapping, which', '28 J. S. Speagle': "serves as a generic and robust tool for attempting to simulate the impact of sampling uncertainties with limited support (Efron 1979). Since in most cases we have many thousands of samples from our distribution and sample with K > 100 live points, Nested Sampling is almost always in a regime where bootstrapping should be viable.  \nNaively, we might expect to simply be able simulate values of, e.g., ˆ Z by just bootstrapping the underlying set of live points. However, this leads to three immediate complications:  \n- (i) This approach creates multiple samples at the same position. It is unclear how these points need to be ordered to assign them associated prior volumes.\n- (ii) This approach conserves the total number of samples N , which clearly must be allowed to change if we really want to simulate from all possible live point paths (with varying path-lengths).\n- (iii) This approach can leave out samples initially drawn from the prior. These points are crucial for establishing the normalization needed to estimate the evidence, and so removing them drastically distorts our evidence estimates.  \nWe address each of these in turn.  \nFirst, the ambiguous ordering, while at first glance a serious issue, is actually a non-concern since the impact on any derived quantity is actually completely insensitive to the ordering. For the evidence, since the likelihood L' GLYPH<2> ' is identical among the points, their contribution to ˆ Z will remain unchanged. Likewise, because they occupy the same position GLYPH<2> , their contribution to the posterior estimate ˆ P' GLYPH<2> ' is also unchanged. This implies that any ordering scheme (e.g., random) will suffice.  \nTo resolve the second issue, we now turn to the problem of simulating all possible live point paths along with their possibly varying path-lengths. Bootstrapping over all the samples by construction destroys this information by ignoring the paths of each individual live point. Analogous to the discussion in § A3, we can characterize these individual paths as being the collection of K lists of positions GLYPH<2> » j … 1 ! GLYPH<1> GLYPH<1> GLYPH<1> ! GLYPH<2> » j … N j traversed by each live point. Sampling from the space of all possible live point paths thus is equivalent to bootstrapping from these individual K 'strands' and then merging the K resampled strands f : : : ; f GLYPH<2> » j 0 … 1 ! GLYPH<1> GLYPH<1> GLYPH<1> !  \n»  \nj  \n0  \nN  \n…  \nj  \n0  \nUnfortunately, this procedure still can run afoul of the third issue when the number of live points K i is not constant. Going back to the discussion in § A3 and the Iterative Dynamic Nested Sampling scheme outlined in Algorithm 3, we see that increasing the number of live points at some iteration i > 1 means that those additional live points were sampled interior to the prior at some associated likelihood threshold L' GLYPH<2> i ' . Since these live points provide no information about the overall normalization (only the normalization relative to ˆ X i ), they are totally uninformative on their own when it comes to estimating the evidence ˆ Z .  \nTo account for this, we need to perform a stratified bootstrap over the set of K int 'interior' strands (i.e. strands with starting positions interior to the prior) and K anc 'anchor' strands (i.e. strands sampled directly from the prior that 'anchor' the interior strands). Once the set of K int interior strands and K anc strands have been resampled, we can then  \nGLYPH<2>  \ng  \ng  \ninto a new Nested Sampling run.  \n; : : :  \nmerge the new collection into a new Nested Sampling run. Following this scheme is then sufficient for simulating the evidence ˆ ˜ Z and posterior ˆ ˜ P' GLYPH<2> ' estimates, where we have used ˆ ˜ Z notation to indicate a we used bootstrapping rather than prior volume simulation.  \nNote that one interesting corollary of our bootstrap estimates is that we expect the total number of samples ˜ N to change. For a sufficient number of live points, this distribution is likely to be roughly Poisson. Assuming that the associated GLYPH<1> ˜ H i GLYPH<25> GLYPH<1> H i and ˜ K i GLYPH<25> K i from our bootstrapped Nested Sampling run are similar to the original, this immediately leads us to an estimate of GLYPH<27> h ln ˆ Z i identical to that in equation (A48). Although it has the exact same form, note that this error term is completely independent from the previous case.", 'A6.3 Combined Uncertainties': "The full uncertainties associated with a given Nested Sampling run involve both the statistical uncertainties described in § A6.1 and the sampling uncertainties described in § A6.2. Simulating from this combined error distribution is straightforward and can be done by the following procedure:  \n- (i) Resample the set of underlying K anc anchor and K int interior strands using stratified bootstrap resampling.\n- (ii) Merge the resampled strands into a single run.\n- (iii) Simulate the values of the prior volumes.  \nWe can then calculate the evidence ˆ ˜ Z 0 and posterior ˆ ˜ P 0 ' GLYPH<2> ' estimates accordingly. The combined uncertainty on the evidence that we estimate from both sources is then roughly  \nGLYPH<27> h ln ˆ Z i GLYPH<24> v u t 2 N GLYPH<213> i = 1 GLYPH<1> H i K i (A49)  \nbased on the identical error estimates derived in § A6.1 and § A6.2.  \nThis paper has been typeset from a T E X/L A T E X file prepared by the author."}