remarkRemark hypothesisHypothesis claimClaim Nested Sampling and Rare Event
EstimationJ. Latz, D. Schneider, and P. Wacker

# Nested Sampling for Uncertainty Quantification and Rare Event
Estimation††thanks: Submitted to the editors .

Jonas Latz Department of Mathematics, University of Manchester, UK ().
<EMAIL_ADDRESS>Doris Schneider Ostbayerische Technische
Hochschule Amberg-Weiden, Germany<EMAIL_ADDRESS>Philipp Wacker
School of Mathematics and Statistics, University of Canterbury, Christchurch,
New Zealand<EMAIL_ADDRESS>
###### Abstract

Nested Sampling is a method for computing the Bayesian evidence, also called
the marginal likelihood, which is the integral of the likelihood with respect
to the prior. More generally, it is a numerical probabilistic quadrature rule.
The main idea of Nested Sampling is to replace a high-dimensional likelihood
integral over parameter space with an integral over the unit line by employing
a push-forward with respect to a suitable transformation. Practically, a set
of active samples ascends the level sets of the integrand function, with the
measure contraction of the super-level sets being statistically estimated. We
justify the validity of this approach for integrands with non-negligible
plateaus, and demonstrate Nested Sampling’s practical effectiveness in
estimating the (log-)probability of rare events.

###### keywords:

Nested Sampling, Marginal Likelihood, Quadrature, Rare Event Estimation

65C05, 28A25, 62-08

##  1 Introduction

Throughout this work, we study the numerical approximation of expected values

(1)
$\mu(\mathcal{L}):=\int_{\Omega}\mathcal{L}(x)\mathop{}\\!\mathrm{d}\mu(x)$

where $(\Omega,\mathcal{A},\mu)$ with $\Omega\subset\mathbb{R}^{d}$ is a
probability space and $\mathcal{L}:\Omega\to\mathbb{R}$ is a measurable
function. We are thinking especially of cases where the evaluation of
$\mathcal{L}$ is computationally costly, e.g., requiring numerical simulations
of partial differential equations or stochastic processes. Here, our methods
are particularly useful, if

  * (i)

we are able to sample independently from $\mu$, but a usual Monte Carlo
approximation with samples from $\mu$ would be inefficient, e.g., where
$\mathcal{L}$ and $\mu$ have most of their masses concentrated in different
regions of $\Omega$ or

  * (ii)

we need to compute $\mu(\mathcal{L}_{\theta})$ for a large number of similar
$\mathcal{L}_{\theta}$, where $\Theta\ni\theta$ is an appropriate parameter
space. Similarity is here expressed in terms of the level sets of the
functions, assuming that
$\\{\\{\mathcal{L}_{\theta}=\alpha\\}:\alpha\in\mathbb{R}\\}\approx\\{\\{\mathcal{L}_{\zeta}=\alpha\\}:\alpha\in\mathbb{R}\\}$,
for $\theta,\zeta\in\Theta$, where we use the usual way of denoting preimages
$\\{\mathcal{L}\in A\\}:=\\{\omega\in\Omega:\mathcal{L}(\omega)\in A\\}$ for
any $A\subseteq\mathbb{R}$.

We give several examples for quadrature problems of the forms (i) and (ii) in
the next section.

###  1.1 Overview of this manuscript

To illustrate the importance of computationally efficient methods for
integrals of type (1), we now give examples of their application in
uncertainty quantification, statistics, physics, and engineering.

#### Model evidence

In Bayesian statistics, we encounter integrals of the described form
$\mu(\mathcal{L})$ when computing the _(model) evidence_ or _marginal
likelihood_. Here, $\mathcal{L}(x)=L(y|x)$, where $L$ describes the likelihood
of the data set $y$ given the parameter $x\in\Omega$; $\mu$ is the prior
distribution. The model evidence is a vital tool in Bayesian model selection
[23], but hard to compute, if the likelihood is highly concentrated. A closely
related problem in statistical physics is that of estimating the partition
function of a system in the thermodynamic equilibrium.

#### Rare event estimation

We consider an event $G\in\mathcal{A}$ with $\mu(G)$ small, say, $\approx
10^{-9}$. Of course, $\mu(G)=\mu(\mathcal{L}),$ in the setting above, when
choosing $\mathcal{L}$ to be the indicator or characteristic function
$\chi_{G}$; defined by $\chi_{G}(\omega)=1,$ if $\omega\in G$ and
$\chi_{G}(\omega)=0,$ otherwise. Rare event estimation is a central tool in
reliability analysis, where failures of physical system need to be estimated,
but also in finance. A straightforward Monte Carlo approximation of a rare
event requires a number of samples of order $p^{-1}$ – unachievable in most
applications; especially when physical simulations are necessary.

#### Moments and the moment-generating function

The moment generating function of $\mathcal{L}$ is defined as
$M_{H}(\theta)=\int_{\Omega}\exp(\theta H(x))\mathrm{d}\mu(x)$, where
$H:\Omega\rightarrow\mathbb{R}$ is a random variable and
$\theta\in\Theta=\mathbb{R}$. The function $M_{H}$ can be used to accurately
represent the probability distribution of the random variable $H(X)$, where
$X\sim\mu$. Here, $\mathcal{L}_{\theta}:=\exp(\theta H)$ satisfies the
condition mentioned above under (ii). Either through the moment generating
function or using the immediate condition, we are also able to compute all
moments of $H(X)$. Here, we can choose $\mathcal{L}_{\theta}=H^{\theta}$, for
$\theta\in\mathbb{N}:=\\{1,2,\ldots\\}$.

#### Cumulative distribution functions and quantiles

The estimation of the cumulative distribution function (cdf)
$F:=\mu(Y\leq\cdot)$ of a random variable $Y:\Omega\rightarrow\mathbb{R}$ is
closely related to the rare event estimation problem, with the difference that
we are now interested in the probability of the event $\\{Y\leq y\\}$ for a
range of values $y\in\mathbb{R}$ (with $y\to-\infty$ corresponding to
increasingly rare events). A particularly interesting problem is the
estimation of quantiles. Here, we need to find $q\in\mathbb{R}$ where
$F(q)=\alpha$ for a pre-determined $\alpha>0$. To find $q$, a bisection method
may require several evaluations of $\mu(Y\leq y)$ with similar sets of the
form $\\{Y\leq y\\}$.

For all these integration tasks, we propose application of the celebrated
_Nested Sampling_ methodology. Nested Sampling was introduced by [36] as a
specialized algorithm for the computation of the Bayesian model evidence and
has been applied successfully in astronomy and computational physics ([3, 13,
14, 25, 26, 29, 30, 39, 40]), biomathematics ([1, 9, 31]), and other fields.
From a bird’s eye view, Nested Sampling is a method of computing the integral
$\int_{\Omega}\mathcal{L}\mathop{}\\!\mathrm{d}\mu$ by converting this high-
dimensional integral into an integral on the unit line $[0,1]$:111The point of
view of [36] is slightly different, interpreting
$\int_{\Omega}\mathcal{L}(x)\mathop{}\\!\mathrm{d}\mu(x)=\int_{0}^{1}\tilde{\mathcal{L}}(X)\mathop{}\\!\mathrm{d}X$,
where $\tilde{\mathcal{L}}$ is an “overloaded” form of the likelihood
$\mathcal{L}$ but we will work with (2), since it allows us to generalise the
methodology considerably.

(2)
$\int_{\Omega}\mathcal{L}(x)\mathop{}\\!\mathrm{d}\mu(x)\stackrel{{\scriptstyle(\clubsuit)}}{{=}}-\int_{0}^{1}\lambda\mathop{}\\!\mathrm{d}X(\lambda)\stackrel{{\scriptstyle(\heartsuit)}}{{\approx}}\sum_{i=1}^{N}\lambda_{i}\cdot\left(X_{i-1}-X_{i}\right)\stackrel{{\scriptstyle(\spadesuit)}}{{\approx}}\sum_{i=1}^{N}\lambda_{i}\cdot
D_{i},$

where $X(\lambda)=\mu(\mathcal{L}>\lambda)$ is the survival function of $X$
and the $(X_{i})_{i=1}^{N}$ are a suitable discretization
$X_{i}=X(\lambda_{i})$ derived further below. This means that the integral in
question is

  1. $(\clubsuit)$

converted into a Riemann-Stieltjes-integral first,

  2. $(\heartsuit)$

which is approximated by a specific (randomized) quadrature rule,

  3. $(\spadesuit)$

with the probabilistic contraction ($X_{i-1}-X_{i}$) being replaced by a
deterministic estimator $D_{i}$.

In some sense, Nested Sampling can be thought of as a Lebesgue-integral
version of stochastic quadrature, in the same way that usual Monte Carlo
quadrature corresponds to Riemann integration (with probabilistically chosen
quadrature points). This point of view becomes clearer from the proof of Lemma
1.1 below. Here, we use a set of particles to approximate increasing (nested)
super-level sets of $\mathcal{L}$. The contraction of the measure of those
super-level sets is then approximated by a deterministic shrinkage estimator.

There is already an extensive body of work regarding the performance of the
second step $(\heartsuit)$, i.e. the quality of the Monte Carlo approximation
of the one-dimensional integral and a series of improvements have been
proposed [19, 12, 18, 34]. To the best of our knowledge there is a gap in the
justification of steps $(\clubsuit)$ and $(\spadesuit)$, i.e. the substitution
of the integral over $\Omega$ with an integral over the unit interval, and the
validity of the deterministic approximation. This manuscript tries to close
this gap in the literature. The crucial issue allowing Nested Sampling to be
applied for general integration tasks (especially those outlined in (i) and
(ii)), is the integration of functions with plateaus. We say, $\mathcal{L}$
has a plateau, if there is some $c\in\mathbb{R}$, for which the level set
$\\{\mathcal{L}=c\\}$ has positive Lebesgue measure.

Computational problems with using, in particular, integrands with non-
vanishing plateaus have been described by [27] and [36], who suggest a
randomization or “labelling” approach, breaking a tie between points with the
same likelihood value. We argue in this manuscript that the version proposed
by [15] is a computationally more suitable way for handling possible plateaus
in the likelihood. While this manuscript did not provide a proof of
correctness, and it is indeed not trivial to do so, we give some justification
for the correctness of this modification of Nested Sampling.

The contributions of this manuscript are, stated succinctly, as follows:

  1. Section 1.2

describes the idea of Nested Sampling and shows that the integral
transformation $(\clubsuit)$ given by
$\int_{\Omega}\mathcal{L}(x)\mathop{}\\!\mathrm{d}\mu(x)=-\int_{0}^{1}\lambda\mathop{}\\!\mathrm{d}X(\lambda)$
is valid even for very general functions $\mathcal{L}$ with plateaus. This is
achieved by viewing Nested Sampling as a numerical quadrature of a Riemann-
Stieltjes-integral, a perspective that apparently has not been explored so
far.

  2. Section 2

justifies with some examples why the modification proposed by [15] is a good
idea in order to make sure that the deterministic approximation $(\spadesuit)$
is valid. This is not entirely straightforward to see and depends on slightly
obscure properties of survival functions. Essentially, the deterministic
estimate for the contraction $X(\lambda_{i})-X(\lambda_{i+1})$ of the
available prior mass has to be justified.

  3. Section 3

proposes a way of unifying large classes of integrals of the form (1) by means
of level set surrogates, allowing efficient application of Nested Sampling to
rare event estimation, approximation of the survival function of a random
variable, and simultaneous estimation of higher-order moments.

  4. Section 4

demonstrates the performance of Nested Sampling in the tasks outlined in
Section 3.

###  1.2 The paradigm of Nested Sampling

We begin with a quick stand-alone explanation of Nested Sampling, i.e. how
exactly (2) can be derived. This derivation is partially restricted (for
purely pedagogical reasons) to the case that $\mathcal{L}$ does not have any
plateaus of non-zero measure. We will later drop this assumption.

######  Assumption 1.

We consider a probability space $(\Omega,\mathcal{A},\mu)$. Let
$\mathcal{L}:\Omega\to\mathbb{R}$ be a measurable function. We define
$\mu^{\lambda}$ as the conditional measure given $\\{\mathcal{L}>\lambda\\}$,
i.e.
$\mu^{\lambda}(A):=\frac{\mu(A\cap\\{\mathcal{L}>\lambda\\})}{\mu(\mathcal{L}>\lambda)}$
and we assume that

  1. (a)

Positivity: $\mathcal{L}>0$ almost surely, and

  2. (b)

Finite Moments: $\mathbb{E}\mathcal{L}^{p}<\infty$ for some $p>1$.

Note that this also works with any $\hat{\mathcal{L}}$ bounded from below by
any $C\in\mathbb{R}$, which can be reduced to the case $\mathcal{L}>0$ by
setting $\mathcal{L}=\hat{\mathcal{L}}-C$ and seeing that
$\int\hat{\mathcal{L}}\mathop{}\\!\mathrm{d}\mu=\int\mathcal{L}\mathop{}\\!\mathrm{d}\mu+C$.

######  Lemma 1.1 (validity of $\clubsuit$).

If Assumption 1 holds, we have

$\int_{\Omega}\mathcal{L}(x)\mathop{}\\!\mathrm{d}\mu(x)\stackrel{{\scriptstyle\clubsuit}}{{=}}\int_{0}^{\infty}\lambda\mathop{}\\!\mathrm{d}(-X(\lambda)),$

where $X(\lambda):=\mu(\mathcal{L}>\lambda)$ is the so-called _survival
function_ of $\mathcal{L}$.

######  Proof 1.2.

We define $X(\lambda)$ as stated above and see that

$\displaystyle\int_{\Omega}\mathcal{L}(x)\mathop{}\\!\mathrm{d}\mu(x)$
$\displaystyle=\mathbb{E}^{\mu}\mathcal{L}=\int_{0}^{\infty}\mu(\mathcal{L}>\lambda)\mathop{}\\!\mathrm{d}\lambda=\int_{0}^{\infty}X(\lambda)\mathop{}\\!\mathrm{d}\lambda.$
The Riemann-Stieltjes integral $\int_{0}^{r}\lambda\mathrm{d}X(\lambda)$ is
well-defined and we can rewrite
$\displaystyle\int_{0}^{\infty}X(\lambda)\mathop{}\\!\mathrm{d}\lambda$
$\displaystyle=\lim_{r\to\infty}\int_{0}^{r}X(\lambda)\mathop{}\\!\mathrm{d}\lambda=\lim_{r\to\infty}X(r)r-X(0)0+\lim_{r\to\infty}\int_{0}^{r}\lambda\mathop{}\\!\mathrm{d}(-X(\lambda)).$
By Markov’s inequality,
$0\leq\lim_{r\to\infty}X(r)r\leq\lim_{r\to\infty}\mathbb{E}\mathcal{L}^{p}r^{1-p}=0$.
Then $\displaystyle\lim_{r\to\infty}$ $\displaystyle
X(r)r-X(0)0+\lim_{r\to\infty}\int_{0}^{r}\lambda\mathop{}\\!\mathrm{d}(-X(\lambda))=\int_{0}^{\infty}\lambda\mathop{}\\!\mathrm{d}(-X(\lambda)).\qquad$

#### Approximating the one-dimensional integral ($\heartsuit$)

By definition, Riemann-Stieltjes integrals can be approximated via

$\displaystyle\int_{0}^{\infty}\lambda\mathop{}\\!\mathrm{d}(-X(\lambda))\approx\sum_{i=0}^{N-1}\lambda_{i+1}\left[X(\lambda_{i})-X(\lambda_{i+1})\right]$

for some choice of $0=\lambda_{0}<\lambda_{1}<\cdots<\lambda_{N}$ where we set
$X(\lambda_{0})=1$. Unfortunately, if we just set some fixed points
$\lambda_{i}$, we do not have any straight-forward way of calculating
$\left[X(\lambda_{i-1})-X(\lambda_{i})\right]=\mu(\mathcal{L}\in(\lambda_{i-1},\lambda_{i}])$.
The trick of Nested Sampling is now to find a clever way of choosing the
points $\lambda_{i}$ such that the quantity $X(\lambda_{i-1})-X(\lambda_{i})$
can be statistically estimated. This is the content of Lemma 1.3, which relies
on the following additional assumption to hold.

######  Assumption 2.

We assume that Assumption 1 holds, and additionally

  1. (a)

No plateaus: $\mu(\mathcal{L}=\lambda)=0$ for all $\lambda\in\mathbb{R}$
(equivalently, the cumulative distribution function of $\mathcal{L}$ is
continuous), and

  2. (b)

Super-level sampling: we can efficiently sample from $\mu^{\lambda}$.

Assumption 2(b) seems (and often is) restrictive, but is in practice fulfilled
in the context of Nested Sampling, as we will see further below.

We start with $\lambda_{0}=\inf\mathcal{L}=0$, i.e. $\mu^{\lambda_{0}}=\mu$.
We now describe how to inductively choose $\lambda_{i+1}$ from given
$\lambda_{i}$: We generate independent and identically distributed (i.i.d.)
samples $\\{x^{k}\\}_{k=1}^{J}\sim\mu^{\lambda_{i}}$, which is possible
according to Assumption 1(b). We then set

(3) $\lambda_{i+1}=\min_{k=1}^{J}\mathcal{L}(x^{k})$

and

(4) $\xi_{i+1}:=\max_{k}X(\mathcal{L}(x^{k}))=X(\lambda_{i+1}).$

Equality holds due to the fact that
$\operatorname{argmin}_{k}\mathcal{L}(x^{k})=\operatorname{argmax}_{k}X(\mathcal{L}(x^{k}))$
because $X$ is non-increasing. The next lemma characterizes what we can say
about the contraction $X(\lambda_{i})-X(\lambda_{i+1})$.

######  Lemma 1.3 (validity of $\spadesuit$, no-plateau case).

If Assumption 2 holds, then the following statements are true.

  1. 1.

If $\\{x^{k}\\}_{k}\sim\mu$, then $\\{X(\mathcal{L}(x^{k}))\\}_{k}\sim
U[0,1]$. The random variable $\xi_{1}:=\max_{k}X(\mathcal{L}(x^{k}))$ follows
a $\mathrm{Beta}(J,1)$ distribution with mean

$\mathbb{E}^{x^{k}\sim\mu}\left[\max_{k}X(\mathcal{L}(x^{k}))\right]=\frac{J}{J+1}$

Writing $\lambda_{1}:=\min_{k}\mathcal{L}(x^{k})$, i.e.,
$X(\lambda_{1})=\xi_{1}$, we have

$\mathbb{E}^{x^{k}\sim\mu}\left[1-\xi_{1}\right]=\frac{J}{J+1}.$

  2. 2.

Recursively in $i\geq 1$: Consider $\lambda_{i}\geq 0$ and $\xi_{i}$. If
$\\{x^{k}\\}_{k}\sim\mu^{\lambda_{i}}$,222i.e. $x^{k}\sim\mu$ under the
additional condition $\mathcal{L}(x^{k})>\lambda_{i}$ then
$\\{X(\mathcal{L}(x^{k}))\\}_{k}\sim U[0,\xi_{i}]$. The random variable
$\xi_{i+1}:=\max_{k}X(\mathcal{L}(x^{k}))$ follows a rescaled
$\mathrm{Beta}(J,1)$ distribution on the interval $[0,\xi_{i}]$ with mean

$\mathbb{E}^{x^{k}\sim\mu^{\lambda_{i}}}\left[\max_{k}X(\mathcal{L}(x^{k}))\right]=\frac{J}{J+1}\xi_{i}.$

Writing $\lambda_{i+1}:=\min_{k}\mathcal{L}(x^{k})$, i.e.,
$X(\lambda_{i+1})=\xi_{i+1}$,333Note that only the term $\xi_{i+1}$ is
averaged over inside the square brackets, $\xi_{i}$ being fixed from the last
iteration.

$\mathbb{E}^{x^{k}\sim\mu^{\lambda_{i}}}\left[\xi_{i}-\xi_{i+1}\right]=\frac{1}{J+1}\xi_{i}.$

  3. 3.

Defining $\lambda_{0}:=0$ and $\xi_{0}=1$ for consistency,

$\mathbb{E}\left[\xi_{i}-\xi_{i+1}\right]=\frac{J^{i}}{(J+1)^{i+1}},$

where the expectation are over recursive repeated sampling of
$\\{x^{k}\\}\sim\mu^{\lambda_{i}}$ for $i=0,\ldots,i-1$.

The last statement of this lemma motivates the form of the deterministic
estimator $D_{i}=\frac{J^{i}}{(J+1)^{i+1}}$ in (2)$(\spadesuit)$.

######  Proof 1.4.

We define $\mathcal{F}_{\mathcal{L}}$ as the cumulative distribution function
of $\mathcal{L}$. According to [10, proposition 2.(1)], or Lemma A.1,
$(\mathcal{F}_{\mathcal{L}}\circ\mathcal{L})_{\\#}\mu=U[0,1]$.444This means
that if we sample $\\{x^{k}\\}_{k}$ from $\mu$, then the transformed samples
$\\{F_{\mathcal{L}}(\mathcal{L}(x^{k}))\\}_{k}$ will be uniformly distributed
on $[0,1]$. We now compute555I.e. if we sample $\\{x^{k}\\}_{k}$ from
$\mu^{\lambda}$, what is the distribution of
$\\{F_{\mathcal{L}}(\mathcal{L}(x^{k}))\\}_{k}$?
$(\mathcal{F}_{\mathcal{L}}\circ\mathcal{L})_{\\#}\mu^{\lambda}$ by analysing
its cumulative distribution function. Note that
$\mu^{\lambda}(A)=\mu(A\cap\\{\mathcal{L}>\lambda\\})/\mu(\mathcal{L}>\lambda)$.
For brevity of exposition, we interpret
$\mathcal{F}_{\mathcal{L}}\circ\mathcal{L}$ directly as a random variable,
this is to be understood in the sense that we consider the measure
$(\mathcal{F}_{\mathcal{L}}\circ\mathcal{L})_{\\#}\mu$.

Now, using the fact that the inverse of the cumulative distribution function
$\mathcal{F}_{\mathcal{L}}$ is well-defined due to Assumption 2(a),

$\displaystyle\mu^{\lambda}$
$\displaystyle(\mathcal{F}_{\mathcal{L}}\circ\mathcal{L}\leq
r)=\frac{\mu(\mathcal{F}_{\mathcal{L}}\circ\mathcal{L}\leq r\text{ and
}\mathcal{L}>\lambda)}{\mu(\mathcal{L}>\lambda)}$
$\displaystyle=\frac{\mu(\mathcal{F}_{\mathcal{L}}\circ\mathcal{L}\leq r\text{
and
}\mathcal{F}_{\mathcal{L}}\circ\mathcal{L}>\mathcal{F}_{\mathcal{L}}(\lambda))}{1-\mathcal{F}_{\mathcal{L}}(\lambda)}$
$\displaystyle=\frac{\mu(\mathcal{F}_{\mathcal{L}}\circ\mathcal{L}\leq
r)-\mu(\mathcal{F}_{\mathcal{L}}\circ\mathcal{L}\leq\mathcal{F}_{\mathcal{L}}(\lambda))}{1-\mathcal{F}_{\mathcal{L}}(\lambda)}$
$\displaystyle=\frac{r-\mathcal{F}_{\mathcal{L}}(\lambda)}{1-\mathcal{F}_{\mathcal{L}}(\lambda)}\text{
for $r\in[\mathcal{F}_{\mathcal{L}}(\lambda),1]$ and $0$ otherwise,}$

where the last step is due to the fact that
$(\mathcal{F}_{\mathcal{L}}\circ\mathcal{L})_{\\#}\mu=U[0,1]$, i.e.
$\mu(\\{x\in\Omega:\mathcal{F}_{\mathcal{L}}(\mathcal{L}(x))\leq r\\})=r$ for
$r\in[0,1]$. This means that
$(\mathcal{F}_{\mathcal{L}}\circ\mathcal{L})_{\\#}\mu^{\lambda}=U[\mathcal{F}_{\mathcal{L}}(\lambda),1]$.

Now we use the fact that $X(\lambda)=1-\mathcal{F}_{\mathcal{L}}(\lambda)$,
and thus
$(X\circ\mathcal{L})_{\\#}\mu^{\lambda}=U[0,1-\mathcal{F}_{\mathcal{L}}(\lambda)]=U[0,X(\lambda)]$.
The remaining statements follow directly from the fact that the maximum of $J$
uniform distributions is a Beta distribution, recursive reasoning, and
elementary computation.

After this derivation of the variable substitution at the heart of Nested
Sampling (in the restricted case of Assumption 2) we want to point out a few
computational techniques relevant to Nested Sampling.

#### Concrete implementation and computational tricks

Nested Sampling performs a kind of noisy Riemann-Stieltjes integration: We
compute

$\int_{\Omega}\mathcal{L}(x)\mathop{}\\!\mathrm{d}\mu(x)=\int_{0}^{\infty}\lambda\mathop{}\\!\mathrm{d}(-X(\lambda))\approx\sum_{i=0}^{N-1}\lambda_{i+1}\left[X(\lambda_{i})-X(\lambda_{i+1})\right]$

where we can only infer $\left[X(\lambda_{i})-X(\lambda_{i+1})\right]$
statistically, e.g. by approximating it by the unbiased estimator
$\mathbb{E}\left[X(\lambda_{i})-X(\lambda_{i+1})\right]=\frac{J^{i}}{(J+1)^{i+1}}$
obtained from Lemma 1.3.

A relevant computational ingredient is the use of the logsumexp function in
order to avoid numerical underflow: The values of $\lambda_{i}$ are increasing
in $i$, while $X(\lambda_{i})-X(\lambda_{i+1})$ is set to a fixed
exponentially decreasing sequence of the form $\frac{J^{i}}{(J+1)^{i+1}}$.
This means that for most terms in this sum, either of those terms will be very
small. In addition, the value of $\int\mathcal{L}\mathop{}\\!\mathrm{d}\mu$ is
in most applications almost astronomically small, e.g. when $\mathcal{L}$ is a
likelihood in a high-dimensional Bayesian inverse problem with relatively
highly informative data, or when $\mathcal{L}$ is the characteristic function
of an event with very small probability. For this reason, we usually compute
$\log\int\mathcal{L}\mathop{}\\!\mathrm{d}\mu$ instead of
$\int\mathcal{L}\mathop{}\\!\mathrm{d}\mu$. The logsumexp function provides a
way of computing $\log(\sum_{i}\exp(a_{i}))$ iteratively without loss of
information. That this is indeed relevant can be seen with the following
simple example: Set $a_{1}=-800$ and $a_{2}=-801$. Then naive implementation
$\log(\exp(a_{1})+\exp(a_{2}))$ would lead to an error since both terms in the
sum would be evaluated to $0$ in most floating point arithmetics, and inserted
into the logarithm. On the other hand, we can write
$\log(\exp(a_{1})+\exp(a_{2}))=\log(\exp(a_{1})(1+\exp(a_{2}-a_{1})))=a_{1}+\log(1+\exp(a_{2}-a_{1}))$
which can be evaluated without any problem. This way we can iteratively
compute the logarithmic sum
$\log\sum_{i=0}^{N-1}\lambda_{i+1}\left[X(\lambda_{i})-X(\lambda_{i+1})\right]$
without incurring too much loss of accuracy.

A seemingly major stumbling block for the implementation of Nested Sampling is
the issue of Assumption 1(b): How do we sample from $\mu$ under the additional
assumption that $\mathcal{L}>\lambda_{i}$? This is not too big a challenge in
practice and we now describe how this can be done in a recursive manner over
the iterations $i$. The case $i=1$ is trivial, since $\mu^{\lambda_{0}}=\mu$,
i.e. we are sampling from the unconstrained measure
$\\{x^{k}\\}_{k=1}^{J}\sim\mu$. We pick
$\lambda_{1}=\min_{k}\mathcal{L}(x^{k})$, and estimate
$X(\lambda_{0})-X(\lambda_{1})=\frac{1}{J+1}$. This gives us our first term in
the sum, $\lambda_{1}\cdot(X(\lambda_{0})-X(\lambda_{1}))$ and we can
increment the loop counter $i$. Now we would have to sample from
$\mu^{\lambda_{1}}$, which is a priori a nontrivial task. But the samples from
the previous iteration carry information that we can use: Since
$\lambda_{1}=\min_{k}\mathcal{L}(x^{k})<\mathcal{L}(x^{l})$ for
$l\neq\operatorname{argmin}_{k}\mathcal{L}(x^{k})$, all but one of the old
samples give us a good idea of where to look since they already are samples
from $\mu^{\lambda_{0}}$ fulfilling $\mathcal{L}>\lambda_{1}$, and therefore,
are samples from $\mu^{\lambda_{1}}$. We just need to remove one sample (the
particle $x^{k}$ with $k=\operatorname{argmin}_{k}\mathcal{L}(x^{k})$) and
replace it by a sample from $\mu^{\lambda_{1}}$. But since we already have
$J-1$ samples from this measure, we can use an “explorative” sampler like
Slice Sampling or a Metropolis–Hastings variant to start there and obtain a
new sample from $\mu^{\lambda_{1}}$. In fact, we found that a good explorative
method is critical for Nested Sampling’s performance in most applications we
considered.

##  2 Nested Sampling: The plateau case

Nested Sampling was originally introduced as a methodology to compute model
evidences in Bayesian inference. In this case, the function $\mathcal{L}$ is a
likelihood. Likelihood functions $\mathcal{L}$ with plateaus, i.e.
$\mathcal{L}$ having level sets with positive prior measure, appear only
infrequently in the literature. There, we find them when considering rounded
data or general piecewise constant statistical models666In the context of
quantization in electrical engineering, signals are rounded to a grid. The
likelihood of the true signal given such a quantized version is then a uniform
distribution over the range of numbers being rounded to this value [37, 4].
They also appear in more general settings where, e.g., piecewise constant
parameter-data dependencies occur, see, e.g., [22]. that lead to a flat
plateau at the likelihood’s peak. In other scenarios, the likelihood is
supported on a compact domain, e.g. $\\{\mathcal{L}=0\\}$ has positive
measure777Some applications come with explicit knowledge about upper bounds on
the measurement error’s magnitude which corresponds to a plateau (of magnitude
0) on the range of impossible original signals (e.g. in biological
applications; in the context of image processing on grayscale images with
values in the unit interval; or when reading out an analog thermometer by
looking at the nearest labeled tick on the scale: Here, the maximum
measurement error is the distance between adjacent ticks on the scale). In
[44] the authors argue that bounded noise is more realistic in some biological
and physical context and that the choice of the correct noise model has a
large influence on the long-term behaviour of models. In econometrics (e.g.
[17]), lognormal measurement noise is sometimes used as a tractable way of
modelling noise terms guaranteeing non-negative data. See also [8] and the
references therein for more examples of bounded noise in physics, biology, and
engineering..

Figure 1: Left: Plateau as a peak. Right: Plateau by compact support.

However, we explicitly consider the plateau case in this work as we are
interested in much more general integration tasks that go beyond model
evidences. When estimating the probabilities of rare events, $\mathcal{L}$ is
an indicator function and, thus, comprised of exactly two plateaus. Integrands
with plateaus also appear in the estimation of risk measures, such as the
conditional value at risk that has recently gained significance in stochastic
control problems [24]. The statistical-physics-equivalent of estimating the
model evidence in Bayesian inference is the computation of the partition
function. Here, piecewise-constant integrands arise through piecewise-constant
potentials in, e.g., molecular dynamics [20] and statistical thermodynamics
[42]. Outside of these specific, computationally expensive problems, the
general integration of piecewise-constants functions may already be of
interest. If the positions of the discontinuities is not known, many classical
quadrature rules may fail as they usually rely on smoothness in the integrand
and its precise approximation through polynomials [38].

We showed the validity of (2), i.e. of steps $(\clubsuit)$, $(\heartsuit)$,
$(\spadesuit)$ in the no-plateau case. The first equation ($\clubsuit$) was
proven in a way which does not depend on whether the integrand has a plateau
or not. Equation ($\heartsuit$) is just the discretization of the
Riemann–Stieltjes integral, so it remains to investigate whether the
deterministic estimator in ($\spadesuit$) is a reasonable approximation even
in the plateau case.

###  2.1 Correctness of the contraction estimate for integrands with non-
vanishing plateaus

The first proof of correctness of Nested Sampling was given in [11], but does
not consider the problematic case of discontinuous survival functions
(corresponding to plateaus in the integrand, see below). Similarly, [34] gives
a very short and precise derivation via the inverse of the survival function,
but again only for integrands without plateaus of non-negligible measure. In
[6], the authors refer to [5] for justification of the integral
transformation, but in our opinion the issue is quite difficult and deserves a
more careful analysis. At the heart of the matter, Nested Sampling amounts to
integration with respect to the push-forward measure, and the basic idea
hinges on properties of survival functions and generalized inverses which are
not readily found explicitly in the existing literature.

The main difficulty arises when $\mathcal{L}$ has a non-negligible plateauwere
already hinted at (although in the context of the algorithm’s performance) in
the original publication [36] with a discussion of “cliffs” – deemed non-
problematic – and “plateaus”. The author recognizes the difficulty of plateaus
(for reasons laid out below), but rules “[…] even so, it may be possible to
generate [new active samples] efficiently.” The PhD thesis [27] also mentions
this issue and further improves on Skilling’s original suggestion.

We will show that plateaus in ${\mathcal{L}}$ are not only computationally
troublesome but also fundamentally mathematically problematic: If there is a
level $\alpha\in\mathbb{R}$ such that $\mu(\mathcal{L}=\alpha)>0$, then the
following implicit uniformity assumption about Nested Sampling is violated:
“Samples from the prior which are plotted in a $X$-$\mathcal{L}$-diagram are
uniformly distributed along the axis $[0,1]$.” This is what the authors of
[35, section 9.2] mean when they write “In terms of $\xi$, the objects are
uniformly sampled subject to the constraint $\xi<\xi^{\star}$”. This means
that the approximation $(\spadesuit)$ in (2) is invalid, because the
uniformity assumption of $X\circ\mathcal{L}$ in Lemma 1.3 no longer holds.

This uniformity assumption (which is conditional to $\mathcal{L}$ not having
plateaus of non-vanishing mass) is crucial for the statistical contraction
estimate empowering Nested Sampling, and is assumed to hold in every
exposition of Nested Sampling known to us, as with the original manuscript
[36, section 5], but also in expositions since then, [11, section 3], [21,
section 2.1], [19, section 2], to just name a few, usually derived from the
statistics of the largest of $M$ uniform deviates, where $M\in\mathbb{N}$ is
the number of active samples.

In the context of the applications considered in this manuscript, the
integrand $\mathcal{L}$ will usually have a plateau, and thus we need a theory
of what to do in this case. Fortunately, we can prove that this mathematical
problem can be circumvented. Lemma 2.3 shows that
$\mathcal{F}_{\mathcal{L}}\circ\mathcal{L}$ (and thus, $X\circ\mathcal{L}$) is
indeed not uniformly distributed, but gives the correct form of this measure,
which is, essentially, a well-behaved combination of uniform distributions and
Dirac masses.

We recall the following well-known properties of cumulative distribution
functions:

######  Lemma 2.1.

We consider a probability space $(\Omega,\mathcal{A},\mu)$ and a measurable
function $\mathcal{L}:\Omega\to\mathbb{R}$. Then the cumulative distribution
function of $\mathcal{L}$ given by
$\mathcal{F}_{\mathcal{L}}(\lambda)=\mu(\mathcal{L}\leq\lambda)$ is a non-
decreasing càdlàg function, i.e.

  * •

Nondecreasing property: For any $\lambda_{1}<\lambda_{2}$, we have
$\mathcal{F}_{\mathcal{L}}(\lambda_{1})\leq\mathcal{F}_{\mathcal{L}}(\lambda_{2})$.

  * •

Continuity from the right: For all $\lambda$, $\lim_{\varepsilon\searrow
0}\mathcal{F}_{\mathcal{L}}(\lambda+\varepsilon)=\mathcal{F}_{\mathcal{L}}(\lambda)$.

  * •

Limits from the left: For all $\lambda$, the limit
$\mathcal{F}_{\mathcal{L}}(\lambda-):=\lim_{\varepsilon\searrow
0}\mathcal{F}_{\mathcal{L}}(\lambda-\varepsilon)$ is well-defined.

Next we characterise plateaus and cliffs of a random variable: A plateau is a
level set $\lambda^{\star}\in\mathbb{R}$ such that
$\mu(\mathcal{L}=\lambda^{\star})>0$. This can be characterized as a point in
the domain of $\mathcal{F}_{\mathcal{L}}$ with a lack of left continuity:
$\mu(\mathcal{L}=\lambda^{\star})=\mathcal{F}_{\mathcal{L}}(\lambda^{\star})-\mu(\mathcal{F}_{\mathcal{L}}<\lambda^{\star})$.
Indeed, for $\varepsilon>0$,

$\displaystyle\mu(\mathcal{F}_{\mathcal{L}}<\lambda^{\star})\in\left[\mu(\mathcal{F}_{\mathcal{L}}\leq\lambda^{\star}-\varepsilon),\mu(\mathcal{F}_{\mathcal{L}}\leq\lambda^{\star})\right]$

and thus by the sandwich lemma (and letting $\varepsilon\to 0$),

$\displaystyle\mu(\mathcal{F}_{\mathcal{L}}<\lambda^{\star})\in\left[\mathcal{F}_{\mathcal{L}}(\lambda^{\star}-),\mathcal{F}_{\mathcal{L}}(\lambda^{\star})\right]$

and

$\displaystyle\mu(\mathcal{L}=\lambda^{\star})\in[0,\mathcal{F}_{\mathcal{L}}(\lambda^{\star})-\mathcal{F}_{\mathcal{L}}(\lambda^{\star}-)]$

This means that $\mu(\mathcal{L}=\lambda^{\star})>0$ implies
$\mathcal{F}_{\mathcal{L}}(\lambda^{\star})>\mathcal{F}_{\mathcal{L}}(\lambda^{\star}-)$.

A cliff of $\mathcal{L}$ is a pair of levels $\lambda_{1}<\lambda_{2}$, such
that $\mu(\mathcal{L}\in(\lambda_{1},\lambda_{2}))=0$, or equivalently,
$\mu(\mathcal{L}\leq\lambda_{1})=\mu(\mathcal{L}<\lambda_{2})$, i.e. a
connected interval $[\lambda_{1},\lambda_{2})$ such that
$\mathcal{F}_{\mathcal{L}}$ is constant on it. This motivates the following
definition.

######  Definition 2.2.

We consider a probability space $(\Omega,\mathcal{A},\mu)$, a measurable
function $\mathcal{L}:\Omega\to\mathbb{R}$ and its cdf
$\mathcal{F}_{\mathcal{L}}$.

  1. (a)

A level $\lambda^{\star}\in\mathbb{R}$ such that
$\mathcal{F}_{\mathcal{L}}(\lambda^{\star})>\mathcal{F}_{\mathcal{L}}(\lambda^{\star}-)$
is called a plateau of $\mathcal{L}$. We call
$\Delta(\lambda^{\star})=\mathcal{F}_{\mathcal{L}}(\lambda^{\star})-\mathcal{F}_{\mathcal{L}}(\lambda^{\star}-)$
the mass of the plateau and
$\alpha(\lambda^{\star})=\mathcal{F}_{\mathcal{L}}(\lambda^{\star}-)$ the
submass of the plateau.

  2. (b)

A cliff of $\mathcal{L}$ is a largest connected interval $I$ of form
$[\lambda_{1},\lambda_{2})$ such that $\mathcal{F}_{\mathcal{L}}$ is constant
on it.888Note that this does not exclude that
$F_{\mathcal{L}}(\lambda_{1})=\mathcal{F}_{\mathcal{L}}(\lambda_{2})$, i.e.
that the plateau extends to the full interval $[\lambda_{1},\lambda_{2}]$,
this just makes sure that there is no larger interval open to the right of
form $[\lambda_{1},\lambda_{2}+\varepsilon)$ on which
$\mathcal{F}_{\mathcal{L}}$ is constant.

Note that if $\lambda^{\star}$ is not a plateau, then
$\Delta(\lambda^{\star})=0$ and
$\alpha(\lambda^{\star})=\mathcal{F}_{\mathcal{L}}(\mathcal{L}(\lambda^{\star}-))$.
This further implies that $\alpha(\cdot)$ is the left-continuous version of
$\mathcal{F}_{\mathcal{L}}$. This now allows us to characterize the exact form
of $\mathcal{F}_{\mathcal{L}}\circ\mathcal{L}$, interpreted as a random
variable.

######  Lemma 2.3.

We consider a probability space $(\Omega,\mathcal{A},\mu)$ and a measurable
function $\mathcal{L}:\Omega\to\mathbb{R}$ with the following properties:

  1. (a)

Positivity: $\mathcal{L}>0$ almost surely,

  2. (b)

Plateaus: There is a countable number of plateaus $l_{1}<l_{2},\ldots$ with
non-vanishing mass $\mu(\mathcal{L}=l_{i})=\Delta_{i}$, submass
$\mu(\mathcal{L}<l_{i})=\alpha_{i}$, and supermass
$\mu(\mathcal{L}>l_{i})=\beta_{i}$. For all
$\lambda\not\in\\{l_{1},\ldots,l_{M}\\}$, we assume
$\mu(\mathcal{L}=\lambda)=0$.

Now, $(\mathcal{F}_{\mathcal{L}}\circ\mathcal{L})_{\\#}\mu$ and
$(\mathcal{F}_{\mathcal{L}}\circ\mathcal{L})_{\\#}\mu^{\lambda}$ are
probability measures on $[0,1]$ characterized by their cumulative distribution
in the form

(5)
$\displaystyle(\mathcal{F}_{\mathcal{L}}\circ\mathcal{L})_{\\#}\mu([0,\alpha])=\begin{cases}\alpha_{i}&\text{
if }\alpha\in[\alpha_{i},\alpha_{i}+\Delta_{i})\\\
\alpha&\text{else.}\end{cases}$ (6)
$\displaystyle(\mathcal{F}_{\mathcal{L}}\circ\mathcal{L})_{\\#}\mu^{\lambda}([0,\alpha])=\begin{cases}\frac{\alpha_{i}-\mathcal{F}_{\mathcal{L}}(\lambda)}{1-\mathcal{F}_{\mathcal{L}}(\lambda)}&\text{
if
}\alpha\in[\alpha_{i},\alpha_{i}+\Delta_{i}),\alpha>\mathcal{F}_{\mathcal{L}}(\lambda)\\\
\frac{\alpha-\mathcal{F}_{\mathcal{L}}(\lambda)}{1-\mathcal{F}_{\mathcal{L}}(\lambda)}&\text{
if
}\alpha\not\in\bigcup[\alpha_{i},\alpha_{i}+\Delta_{i}),\alpha>\mathcal{F}_{\mathcal{L}}(\lambda)\\\
0&\text{ if }\alpha\leq\mathcal{F}_{\mathcal{L}}(\lambda)\end{cases}$

If we define $X(\lambda):=1-\mathcal{F}_{\mathcal{L}}(\lambda)$,

(7)
$\displaystyle(X\circ\mathcal{L})_{\\#}\mu([0,\alpha])=\begin{cases}\beta_{i}+\Delta_{i}&\text{
if }\alpha\in[\beta_{i},\beta_{i}+\Delta_{i})\\\
\alpha&\text{else}\end{cases}$
$\displaystyle(X\circ\mathcal{L})_{\\#}\mu^{\lambda}([0,\beta])$ (8)
$\displaystyle=\begin{cases}\frac{\beta_{i}+\Delta_{i}}{X(\lambda)}&\text{ if
}\beta\in[\beta_{i},\beta_{i}+\Delta_{i})\text{ and }\beta<X(\lambda)\\\
\frac{\beta}{X(\lambda)}&\text{ if
}\beta\not\in\bigcup[\beta_{i},\beta_{i}+\Delta_{i})\text{ and
}\beta<X(\lambda)\\\ 1&\text{ if }\beta\geq X(\lambda)\end{cases}$

This means that

$\displaystyle(\mathcal{F}_{\mathcal{L}}\circ\mathcal{L})_{\\#}\mu$
$\displaystyle=\mathrm{Unif}\left(\bigcup_{i}[\alpha_{i},\alpha_{i}+\Delta_{i})\right)^{c}+\sum_{i}\Delta_{i}\cdot\delta_{\alpha_{i}+\Delta_{i}}$
$\displaystyle(X\circ\mathcal{L})_{\\#}\mu$
$\displaystyle=\mathrm{Unif}\left(\bigcup_{i}[\beta_{i},\beta_{i}+\Delta_{i})\right)^{c}+\sum_{i}\Delta_{i}\cdot\delta_{\beta_{i}}$

where $r\cdot\delta_{x}$ denotes a Dirac measure of strength $r$ at position
$x$, i.e. $\int f(y)\mathop{}\\!\mathrm{d}(r\cdot\delta_{x}(y))=r\cdot f(x)$.

######  Proof 2.4.

In addition to
$\mathcal{F}_{\mathcal{L}}(\lambda)=\mu(\\{x\in\Omega:\mathcal{L}(x)\leq\lambda\\})$
being the (right-continuous) cumulative distribution function of
$\mathcal{L}$, we define its (left-continuous) version
$\mathcal{F}_{\mathcal{L}}^{\circ}(\lambda)=\mu(\\{x\in\Omega:\mathcal{L}(x)<\lambda\\})$

$\displaystyle(\mathcal{F}_{\mathcal{L}}\circ\mathcal{L})_{\\#}\mu([0,\alpha))$
$\displaystyle=\mu(\\{x\in\Omega:\mathcal{F}_{\mathcal{L}}(\mathcal{L}(x))<\alpha\\})$
Since $\mathcal{F}_{\mathcal{L}}$ is right-continuous, by Lemma A.4(a), we
have equal to
$\displaystyle\mu(\\{x\in\Omega:\mathcal{F}_{\mathcal{L}}(\mathcal{L}(x))<\alpha\\})$
$\displaystyle=\mu(\\{x\in\Omega:\mathcal{L}(x)<\mathcal{F}_{\mathcal{L}}^{-}(\alpha)\\})$
$\displaystyle=\mathcal{F}_{\mathcal{L}}^{\circ}(\mathcal{F}_{\mathcal{L}}^{-}(\alpha))$
$\displaystyle=\mathcal{F}_{\mathcal{L}}^{\circ}((\mathcal{F}_{\mathcal{L}}^{\circ})^{-}(\alpha))$
$\displaystyle=\begin{cases}\alpha_{i}&\text{ if
}\alpha\in(\alpha_{i},\alpha_{i}+\Delta_{i}]\\\ \alpha&\text{else}\end{cases}$

by an application of Lemma A.5, and using the fact that the generalized
inverse does not depend on the continuity of the function considered, i.e.
$\mathcal{F}_{\mathcal{L}}^{-}=(\mathcal{F}_{\mathcal{L}}^{\circ})^{-}$, by
virtue of Lemma A.4(b).

The pushforward of $\mu^{\lambda}$ under
$\mathcal{F}_{\mathcal{L}}\circ\mathcal{L}$ is computed by proceeding with the
same computations as before after seeing that

$\displaystyle(\mathcal{F}_{\mathcal{L}}\circ\mathcal{L})_{\\#}\mu^{\lambda}([0,\alpha))$
$\displaystyle=\frac{\mu(\mathcal{F}_{\mathcal{L}}\circ\mathcal{L}<\alpha\text{
and }\mathcal{L}>\lambda)}{\mu(\mathcal{L}>\lambda)}$
$\displaystyle=\frac{\mathcal{F}_{\mathcal{L}}(\mathcal{F}_{\mathcal{L}}^{-}(\alpha)-)-\mathcal{F}_{\mathcal{L}}(\lambda)}{1-\mathcal{F}_{\mathcal{L}}(\lambda)}.$

The formula for $X$ follows from the following straightforward computation

$\displaystyle(X\circ\mathcal{L})_{\\#}\mu([0,\beta])$
$\displaystyle=1-(\mathcal{F}\circ\mathcal{L})_{\\#}\mu([0,1-\beta))$
$\displaystyle=1-\begin{cases}\alpha_{i}&\text{ if
}1-\beta\in(\alpha_{i},\alpha_{i}+\Delta_{i}]\\\
1-\beta&\text{else}\end{cases}$
$\displaystyle=\begin{cases}\beta_{i}+\Delta_{i}&\text{ if
}\beta\in(\beta_{i},\beta_{i}+\Delta_{i}]\\\ \beta&\text{else}\end{cases}$

where we used the relations $1-\alpha_{i}=\beta_{i}+\Delta_{i}$ and
$1-(\alpha_{i}+\Delta_{i})=\beta_{i}$.

We now present Algorithm 1 which is the modification of the original Nested
Sampling algorithm as proposed in [15], but where we split the “sample
generation” and the “weight contraction/integration” computations into
separate loops. This does not increase the complexity of the algorithm (also,
the weight contraction is purely deterministic and could even be precomputed
beforehand), but allows for more modularity which we will use for rare event
estimation where we swap out the integrand $\mathcal{L}$.

Data: integrand $\mathcal{L}$, probability measure $\mu$, number of live
particles $J$, number of iterations $N$

Result: estimate $Z\approx\int\mathcal{L}\mathop{}\\!\mathrm{d}\mu$, ordered
list of dead samples $R$.

Generate $J$ samples from $\mu$, call them live particles

Call $P$ set of live particles

$R\leftarrow\emptyset$ (Initialize set of dead samples)

$i\leftarrow 0$

while $i<N$ (Sample generation loop) do

$\mathcal{L}^{\star}\leftarrow\min_{j}\\{\mathcal{L}(p^{j}),\text{ live
particles }p^{j}\\}$

for all live particle $p^{j}$ such that
$\mathcal{L}(p^{j})=\mathcal{L}^{\star}$ do

$i\leftarrow i+1$

$R[i]\leftarrow p^{j}$ (Next element in list of dead samples)

$P\leftarrow P\setminus p^{j}$ (Remove particle from live set)

end for

Generate enough samples from $\mu$, conditioned to
$\mathcal{L}>\mathcal{L}^{\star}$, in order to fill $P$ to size $J$ again.

end while

$Z\leftarrow 0$ (Initialize quadrature)

$i\leftarrow 0$

while $i<N$ (Quadrature loop) do

$i\leftarrow i+1$

$x_{i}\leftarrow R[i]$ (Get $i$th dead sample)

$\xi_{i}\leftarrow\left(1-\exp(-1/J)\right)\cdot\exp(-(i-1)/J)$ (deterministic
estimate of contraction)

$Z\leftarrow Z+\xi_{i}\cdot\mathcal{L}(x_{i})$ (Increment quadrature)

end while

Algorithm 1 Modified Nested Sampling according to [15]

######  Remark 2.5.

This algorithm has a remarkable property: The selection process which
generates the list of dead samples $R$ does not depend on the actual magnitude
of the integrand function evaluations, but just on the structure of the level
sets of $\mathcal{L}$. Clearly we could replace $\mathcal{L}$ by
$\mathcal{L}+c$, or $\mathcal{L}\cdot a$, and keep the same set of dead
samples which will be valid dead samples for this modified function, too. Even
more, we can replace $\mathcal{L}$ by a completely different function under
the condition that their level sets are isomorphic. We make this more rigorous
below.

We note that [15] proposed to approximate the unknown quantities
$X_{i}=\mu(\mathcal{L}>\mathcal{L}^{\star})$, where $\mathcal{L}^{\star}$ is
the minimum $\mathcal{L}$ value in the live point ensemble at iteration $i$.
This is an unknown quantity, but it is being approximated by a deterministic
estimator, which is obtained by contracting $X_{i}$ by $e^{-1/J}$ each time we
discard an element of the set of live points. This gives rise to the following
unproven hypothesis, which we supplement with a sketch of validity in a
specific example.

We consider the setting of Lemma 2.3, i.e. we have an integrand $\mathcal{L}$
such that

$\displaystyle(X\circ\mathcal{L})_{\\#}\mu$
$\displaystyle=\mathrm{Unif}\left(\bigcup_{i}[\beta_{i},\beta_{i}+\Delta_{i})\right)^{c}+\sum_{i}\Delta_{i}\cdot\delta_{\beta_{i}}.$

is a measure on $[0,1]$.

We further set $\mathcal{L}^{\star}=\min_{j}\mathcal{L}(x^{j})$, i.e., have
that $X(\mathcal{L}^{\star})=\max_{j=1}^{J}X(\mathcal{L}(x^{j}))$. Then we can
approximate $X(\mathcal{L}^{\star})$ deterministically by
$X(\mathcal{L}^{\star})\approx e^{-n/J}$, where
$n=\operatorname{Card}(\\{x^{j}:\mathcal{L}(x^{j})=\mathcal{L}^{\star}\\})$.

######  Proof 2.6 (Sketch of the validity of Section 2.1).

Figure 2: $\mathcal{L}_{1}$ with a plateau at its minimum, and
$\mathcal{L}_{2}$ with a plateau at its maximum.

We consider the two examples in Figure 2. We denote them by $\mathcal{L}_{1}$
and $\mathcal{L}_{2}$. We assume that the measure $\mu$ is absolutely
continuous so that indeed
$\mu(\mathcal{L}_{1}=\inf\mathcal{L}_{1})=\Delta_{1}$,
$\mu(\mathcal{L}_{1}>\inf\mathcal{L}_{1})=\beta_{1}=1-\Delta_{1}$ and
$\mu(\mathcal{L}_{2}=\sup\mathcal{L}_{2})=\Delta_{2}$,
$\mu(\mathcal{L}_{2}>\sup\mathcal{L}_{2})=\beta_{2}=0$. Hence,

$\displaystyle(X\circ\mathcal{L}_{1})_{\\#}\mu$
$\displaystyle=\mathrm{Unif}[0,1-\Delta_{1})+\Delta_{1}\cdot\delta_{\beta_{1}}$
$\displaystyle(X\circ\mathcal{L}_{2})_{\\#}\mu$
$\displaystyle=\Delta_{2}\cdot\delta_{0}+\mathrm{Unif}[\Delta_{2},1]$

and the cdfs can be written down as follows:

$\displaystyle\mu\left(X\circ\mathcal{L}_{1}\leq\alpha\right)$
$\displaystyle=\begin{cases}\alpha:&\text{ if }\alpha<\beta_{1}\\\ 1&\text{ if
}\alpha\geq\beta_{1}\end{cases}$
$\displaystyle\mu\left(X\circ\mathcal{L}_{2}\leq\alpha\right)$
$\displaystyle=\begin{cases}\Delta_{2}&\text{ if }\alpha\leq\Delta_{2}\\\
\alpha:&\text{ if }\alpha>\Delta_{2}\end{cases}$

The cdf of the maximum of $J$ i.i.d. samples is just the $J$-th power of their
individual cdf, hence

$\displaystyle\mu\left(\max_{J}X(\mathcal{L}_{1}(x^{k}))\leq\alpha\right)$
$\displaystyle=\begin{cases}\alpha^{J}:&\text{ if }\alpha<\beta_{1}\\\
1&\text{ if }\alpha\geq\beta_{1}\end{cases}$
$\displaystyle\mu\left(\max_{J}X(\mathcal{L}_{2}(x^{k}))\leq\alpha\right)$
$\displaystyle=\begin{cases}\Delta_{2}^{J}&\text{ if }\alpha\leq\Delta_{2}\\\
\alpha^{J}:&\text{ if }\alpha>\Delta_{2}\end{cases}$

We start with the first example, $\mathcal{L}_{1}$, and we more specifically
demand that $J\gg 1$, $\Delta\gg 0$ (or more generally $(1-\Delta)^{J}\ll 1$).
This means that with probability $1-(1-\Delta)^{J}$, the maximum will be
attained at $1-\Delta$. Nested Sampling employs a contraction of $e^{-1/J}$
for every particle
$X(\mathcal{L}_{1}(x^{i}))=\max_{j}X(\mathcal{L}_{1}(x^{j}))=:\xi$, leading to
a total contraction of $e^{-n/J}$, where $N={\rm
card}\\{x^{j}:X(\mathcal{L}_{1}(x^{j}))=\xi\\}$. $N$ is binomially
distributed: $B\sim{\rm Bin}(J,\Delta)$. This means that the expected number
of particles all being the joint maximum is $\mathbb{E}N=J\cdot\Delta$.
Ignoring the (very improbable) possibility that the maximum is not attained at
$1-\Delta$, the following happens during the foreach loop of considering this
plateau: The accessible prior volume is contracted by $e^{-1/J}$ a number of
$n=J\cdot\Delta$ on average, i.e. on average, it will have been reduced to
$X_{n}=e^{-\Delta}$, and the evidence will have been increased (telescopic
sum!) to $Z=(1-e^{-\Delta})\mathcal{L}^{\star}$. This is exactly what we would
expect: The plateau has prior contribution $\Delta$, which is roughly equal to
$(1-e^{-\Delta})$ (at least for $\Delta$ not too large), with value
$\mathcal{L}^{\star}$. This correctly (to first order in $\Delta$) handles the
plateau, and then Nested Sampling proceeds with the rest of the integrand
which does not have any plateaus (and thus this is reduced to the known case).

Now to the second case where $(1-\Delta)^{J}\gg 0$, i.e., the maximum is with
very high probability not attained in the Dirac component. Then calculations
are even easier: $N=1$ and $\xi$ is Beta distributed with mean
$\frac{J}{J+1}(1-\Delta)$, which matches the contraction $\exp(-1/J)$ proposed
by Nested Sampling to first order (in $J^{-1}$).

Second example: The probability that the maximum is attained in the (left)
Dirac is given by $\Delta^{J}$. We can again distiguish two cases: Case I
($\Delta^{J}\approx 0$) is again straightforward and we obtain the usual
contraction rate. Case II: $\Delta^{J}\ll 0$ means that there is a nontrivial
chance that the maximum (and hence all particles) land in $x=0$. This means
that the contraction rate is
$\left(\frac{J}{J+1}\right)^{J}\simeq\exp(-\frac{J}{J+1})\simeq 1$, which is
indeed correct since the contraction is the full interval.

##  3 Application to Rare Event Estimation and beyond

In this section we mainly demonstrate how Nested Sampling can be used to
efficiently compute rare event probabilities: to this end, we introduce our
methodology of surrogate functions that we can also employ to compute moments
and to approximate cumulative distribution functions, which we come back to in
Section 3.2. The estimation of probabilities of rare events is central in many
disciplines, especially in civil engineering. Nested Sampling in rare event
estimation has previously been discussed by [16, 43]. Especially the fact that
Nested Sampling is usually implemented to work with logarithmic quantities
only allows for a very stable computation. More traditional methods for rare
event estimation include FORM [32], subset simulation [2], sequential
importance sampling [28], and the cross-entropy method [33]. We now introduce
our Nested-Sampling-based rare event estimator.

We consider now the problem of computing

(9) $\int_{\Omega}\chi_{A}\mathop{}\\!\mathrm{d}\mu$

where $A$ is a rare event, i.e. we assume that $\mu(A)$ is too small for
standard Monte-Carlo estimation. In the previous section we have seen that it
is possible to compute this integral via Nested Sampling, even though the
integrand $\chi_{A}$ consists solely of plateaus. However, this would be
inefficient: the worker is stuck in the first while-loop of Algorithm 1 until
sampling $J$ particles that are in the rare event. To solve this problem, we
now propose a variant of Nested Sampling to compute this object in practice.
The key idea will be to swap out a regularization of the characteristic
function by a level-set surrogate, which is a different function having
matching super-level sets. The level-set surrogate can be chosen such that it
allows for more efficient re-sampling than $\chi_{A}$ during Nested Sampling.
The samples and contractions obtained via the level-set surrogate can be
employed to calculate $\int_{\Omega}\chi_{A}\mathop{}\\!\mathrm{d}\mu$
afterwards.

###  3.1 A unifying computational trick: Surrogate Integrands

######  Definition 3.1 (Level-set Surrogate).

Let $f:\Omega\to\mathbb{R}$ be a continuous mapping. We call $g$ a level-set
surrogate for $f$, if for all $\lambda\in\mathbb{R}$ there exists a
$\kappa\in\mathbb{R}$ such that
$\\{x\in\Omega:f(x)>\lambda\\}=\\{x\in\Omega:g(x)>\kappa\\}$.

######  Lemma 3.2.

Let $\\{x_{i}\\}_{i=1}^{N}$, $\\{\xi_{i}\\}_{i=1}^{N}$ be the set of dead
samples and estimated contractions obtained from one run of Nested Sampling
for the computation of $\int f\mathop{}\\!\mathrm{d}\mu$. In particular,

$\int f\mathop{}\\!\mathrm{d}\mu\approx\sum_{i=1}^{N}f(x_{i})\xi_{i}.$

If $g$ is a level-set surrogate for $f$, then we can use the same set of dead
samples to approximate $\int g\mathop{}\\!\mathrm{d}\mu$:

$\int g\mathop{}\\!\mathrm{d}\mu\approx\sum_{i=1}^{N}g(x_{i})\xi_{i}$

######  Proof 3.3.

From Section 2, we already know that the estimation of the contraction
$\\{\xi_{i}\\}_{i=1}^{N}$ is independent of the integrand. We further know
that $g$ is a level-set surrogate of $f$ and, thus, it holds
$\operatorname{argmin}_{k}f(s^{k})=\operatorname{argmin}_{k}g(s^{k})$ for any
i.i.d sample set $\\{s^{k}\\}_{k=1}^{J}\sim\mu^{\lambda_{i}}$. Therefore, the
lemma follows directly as long as it is ensured that the process of dead
sample selection in Nested Sampling is unchanged by replacing $f$ with $g$.

We can exploit this observation in at least two ways:

  * •

Instead of computing an integral of form
$\int\chi_{A}\mathop{}\\!\mathrm{d}\mu$, we can construct an integrand $g$
with the property that there exist a value $\kappa$ such that
$A=\\{g>\kappa\\}$. Then $g$ is a level-set surrogate for $\chi_{A}$ and we
can construct samples $\\{x_{i}\\}$ and contractions $\\{\xi_{i}\\}$ using
Nested Sampling performed with $g$. For example, if the integration domain is
one-dimensional and $A=(3,\infty)$, then $g(x)=x$ is a valid surrogate
function (and Nested Sampling tends to work much better with this integrand
since it avoids computational problems associated with plateaus).

  * •

If there are multiple integrals of interest with the property that all
integrands are level-set surrogates for each other, then we can only need to
compute samples and contractions once, and plug them into each integrand. This
will be used for the simultaneous computation of moments as well as for the
approximation of distribution functions. A further example would be the
simultaneous computation of the probabilities of a sequence of nested rare
events.

Therefore, Lemma 3.2 tells us how to handle “demanding” integrals more
efficiently: we can approximate the calculation of
$\int_{\Omega}\chi_{A}\mathop{}\\!\mathrm{d}\mu$ by the calculation of the
integral of an appropriate regularized level-set surrogate.999Actually, it is
already sufficient to perform Nested Sampling employing the level-set
surrogate until iteration $K$ which is the first iteration in which the
removed particle is a rare event ($x_{K}\in A$). Then, the integral can be
approximated by $1-\sum\xi_{K}$. This has mainly the benefit that the
computational effort of re-sampling from $\mu^{\lambda_{i}}$ is reduced due to
the fact that an appropriate level-set surrogate allows for a more directed
search of new samples within the “exploration” method (needed for the
generation of new samples satisfying $\mathcal{L}>\lambda^{\star}$), e.g.,
gradient based methods.

Example: We consider the case of $\mu=\mathcal{N}(0,1)$ and
$\mathcal{L}=\chi_{[a,\infty)}$. So
$\int\mathcal{L}d\mu=\mathcal{N}(0,1)[a,\infty)$ for $a\gg 1$ is the
probability of the rare event of obtaining a large ($\geq a$) positive number
when sampling from a standard normal distribution. Then a suitable surrogate
integrand is $g(x)=x$ as $\mathcal{L}$ can be approximated arbitrarily well
(with $h\to 0$) by a sigmoid function of form
$f(x)=\tfrac{1}{2}\tanh((x-a)/h)+\tfrac{1}{2}$, and level sets of $f$ are
isomorphic to level sets of $g$.

###  3.2 Survival functions and higher-order moments

In addition to the rare event probability, we can also infer moments and the
cumulative distribution function from the dead samples with (almost) no
additional computational cost.

#### Calculation of moments

The $k$-th moment is given by
$m_{k}=\int_{\Omega}x^{k}\cdot\mathcal{L}(x)\mathop{}\\!\mathrm{d}\mu(x)$.
Here, we can use the particles $\\{x_{i}\\}_{i=1}^{N}$ which were removed from
the set of particles in iteration $i$ to approximate this integral:

$\displaystyle
m_{k}=\int_{\Omega}x^{k}\cdot\mathcal{L}(x)\mathop{}\\!\mathrm{d}\mu(x)\approx\frac{1}{Z}\left({x}_{i}\right)^{k}\cdot\mathcal{L}({x}_{i})\cdot\left(X(\lambda_{i-1})-X(\lambda_{i})\right).$

In fact, this reduces to a weighted sum of the removed particles with
normalized weights $\frac{1}{Z}w_{i}$ with
$w_{i}=\xi_{i}\cdot\mathcal{L}({x}_{i})$, see also [36]. The quantities
${x}_{i},\mathcal{L}({x}_{i})$ and $w_{i}$ are already computed while
iterating, see Algorithm 1. Therefore, the calculation of moment $m_{k}$
requires almost no additional computational effort. We know that the higher
the moment, the higher the contribution of samples which are “far away in the
region of the rare event”. This is relativized by the weights $w_{i}$.
Nevertheless, the precision of the moments decreases with increasing $k$.

#### Approximation of a cumulative distribution function

The cdf can be approximated by the empirical cdf derived from the $N$ weighted
particles $\tilde{x}_{i}$:

$\displaystyle
F(t)\approx\hat{F}_{N}(t)=\sum_{i=1}^{N}\left(X(\lambda_{i-1})-X(\lambda_{i})\right)\cdot\mathcal{L}(\tilde{x}_{i})\cdot\chi_{\tilde{x}_{i}\leq
t}(t),$

where again $w_{i}=(X(\lambda)-X(\lambda))\cdot\mathcal{L}(\tilde{x}_{i})$ is
the weight of sample $\tilde{x}_{i}$. As for the moments, the quantities
$\tilde{x}_{i},\mathcal{L}(\tilde{x}_{i})$ and $w_{i}$ are already computed
while iterating. Therefore, almost no additional effort is spent on the
approximation of the cdf after performing Nested Sampling.

##  4 Numerical Experiments

###  4.1 Elementary one-dimensional densities

We consider first the elementary example of computing the rare event
probability $P_{a}=\mu_{i}([a,\infty))$ in two settings: $\mu_{1}$ being a
standard Gaussian probability measure, and $\mu_{2}$ being a Cauchy
probability measure. The results of our implementation with Nested Sampling
can then be directly compared to the ground truth solution given by
$P_{a}=1-\Phi(a)$, where $\Phi$ is the Gaussian cumulative distribution
function, and similarly for the Cauchy setting. The integrand is thus given by
$\mathcal{L}(x)=\chi_{[a,\infty)}(x)$. As a surrogate function, we just use
the linear function ${\mathcal{L}_{\text{surrogate}}}(x)=x$, as motivated
earlier.

#### 1d standard Gaussian: $P_{6}$ and moments

The log probability of the rare event $X>6$ with standard Gaussian random
variable $X$ is $\log P_{6}\approx-20.74$. We approximate it via Nested
Sampling with $J\in\\{5,15,50\\}$ particles. The estimated logprobabilites of
200 independent runs are summarized in Figure 3. Additionally, we give the
approximation of the (log)-survival function in Figure 3 which is achieved
almost without additional computational effort.

#### 1d Cauchy: $P_{100}$ and cdf

The probability of the rare event $X>100$ with Cauchy random variable $X$ is
$\log P_{6}\approx-5.75$. The approximation results obtained from Nested
Sampling using $J\in\\{5,50,500\\}$ particles are shown in Figure 3 and the
approximation of the (log-)survival function is shown in Figure 3.

Figure 3: Top row: Estimation of rare event log-probability
$\log\mu_{i}([a_{i},\infty))$ for different sizes $J$ of active set. Left:
$\mu_{1}=\mathcal{N}(0,1)$, $a=6$, with $J\in\\{5,15,10\\}$. Right:
$\mu_{2}=\mathrm{Cauchy}$, $a=100$, with $J\in\\{5,50,500\\}$. True value is
shown as a blue line. Second and third row: True and estimated (via Nested
Sampling) (log-)survival function of $\mathcal{N}(0,1)$. Bottom two rows: True
and estimated (via Nested Sampling) (log-)survival function of a Cauchy
probability measure.

###  4.2 Diffusion in a double-well potential

We consider a double-well potential $V(x)=-\frac{a}{2}x^{2}+\frac{b}{4}x^{4}$,
where we set $a=2$ and $b=0.5$. Then $V$ has two wells at $\pm\sqrt{a/b}=\pm
2$. We consider the overdamped Langevin dynamics of a particle (initialized to
start in one of the wells) in this potential governed by

$\mathop{}\\!\mathrm{d}x_{t}=-\nabla
V(x_{t})+\sigma\mathop{}\\!\mathrm{d}W_{t},\quad x_{0}=2.$

Figure 4: Rare Event estimation for double-well diffusion. Top row: $J=10$.
Bottom row: $J=100$. Left column: accumulated value of $\log Z$ during 25
independent runs of Nested Sampling (boxplot shows distribution of final value
of computed rare event probability). Right column: Visualization of samples
(i.e. paths in double well potential) removed from active set along with
iteration of Nested Sampling. Color of path varies with iteration from dark
(early) to light (late).

We are interested in the probability of the path $\\{x_{t}\\}_{t\in[0,T]}$
leaving its initial well at $x_{0}=2$ within $t\in[0,10]$. This is a rare
event, and due to the attractive nature of the second well we can specify it
via $A=\\{\min_{t\in[0,T]}x_{t}\leq-1.5\\}$, i.e.
$\mu(A)=\int_{\Omega}\chi_{A}(\\{x_{t}\\}_{t})\mathop{}\\!\mathrm{d}\mu(\\{x_{t}\\}_{t})$
if $\mu$ is the distribution of the stochastic process $\\{x_{t}\\}_{t}$. As a
level-set surrogate function we can choose
$g(\\{x_{t}\\}_{t})=\min_{t\in[0,T]}x_{t}$. Figure 4 shows the evolution of
accumulated log-probability of $25$ independent runs of Nested Samplings for
ensemble sizes $J=10$ and $J=100$ each, as well as (for one of the $25$ runs)
a visualization of all discarded samples (which correspond to specific
diffusion paths) created during a specific run of Nested Sampling. Darker
paths in the right part of Figure 4 correspond to typical (prior-like)
simulation outcomes, with lighter paths corresponding to (rarer) paths closer
to the rare event of switching to the secondary potential.

It can be observed that a larger number of ensemble members $J$ decreases the
variance of the Nested Sampling estimator of the log probability, but requires
a higher number of iterations until convergence.

###  4.3 Loaded beam with random inclusions

We now study a rare event estimation problem that arises in an engineering
context. We consider a beam (1D) that consists of a generally homogeneous
material that is subject to small inclusions. These inclusions affect the
beam’s flexibility $F$. We now affix the beam horizontally on one side and
attach a small point load on the other side. Given a stochastic model for the
inclusions, we are interested in the probability of the beam vertically
deflecting beyond a certain point.

The beam has length $L$; the flexibility is given by a sum of a constant $c>0$
and a Poisson point process on $[0,\infty)$. Then, the flexibility is given as
a measure

$F(\mathrm{d}t)=c\mathrm{d}t+\sum_{i=1}^{\infty}w\delta(\mathrm{d}t-T_{i}),$

where $T_{i}=\sum_{j=1}^{i}\Delta T_{i}$, $\Delta T_{i}$ are independent
exponentially distributed random variables with rate $\lambda>0$
$(i\in\mathbb{N})$, and $w$ is a weight. We choose this model for the
flexibility to reflect that the inclusions are spatially extremely small, but
have a significant influence. The deflection response $d$ is given through an
Euler–Bernoulli equation and reads:

$d(x)=-P\int_{0}^{x}\int_{0}^{s}(L-t)F(\mathrm{d}t)\mathrm{d}s\qquad(x\in[0,L]).$

This integral can be solved analytically:

$d(x)=-P\left(\frac{cLx^{2}}{2}-\frac{cx^{3}}{6}+\sum_{i:T_{i}\leq
x}w(L-T_{i})x\right)\qquad(x\in[0,L]).$

For our experiment, we consider a beam of length $L=5$, usual flexibility
constant $c=1$, load $P:=0.01$, additional flexibility at inclusions
$w:=0.05$, and inclusion rate $\lambda=1$. We assume that the rare event
occurs when the beam deflects beyond the point $-0.55$, which is about 30%
lower than the deflection under no inclusions: $-0.417$.

To discretise the Poisson point process given above, we assume a maximum
number of inclusions of $20$; given inclusion rate and length of the rod, the
probability to have more than 20 inclusions is very small; we have
$\mathrm{Poisson}(5)([21,\infty))=8.11\cdot 10^{-8}$. This is two orders of
magnitude smaller than the probability of the rare event, which we estimated
with $10^{8}$ Monte Carlo samples and give it in the rightmost panel in Figure
5. The top panel shows a beam without inclusions as well one with rare event
inclusions leading to displacement beyond the critical threshold. The bottom
panel visualizes (similarly to our simulations for the double-well potential)
displacements corresponding to samples discarded during the iteration of
Nested Sampling. The surrogate function used was the negative displacement of
the beam.

Figure 5: Top: Loaded beam with and without inclusions. Bottom: Samples
generated by Nested Sampling during iteration (dark = early, light = late
iteration). Right: Logprobability of rare event (deflection exceeding
$-0.55$), as computed by 20 independent runs of Nested Sampling (Boxplot), in
comparison to brute force Monte Carlo simulation with $10^{9}$ runs
(bootstrapped - horizontal line).

###  4.4 Random Walk in a labyrinth

We consider a discrete state space
$\mathcal{S}:=\\{0,\ldots,N-1\\}\times\\{0,\ldots,N-1\\}$ in the form of a
two-dimensional grid with “prohibited states”(walls), forming a labyrinth as
in Figure 6. On this labyrinth we set a random walk of length $K$, starting at
$(0,0)$. In each step, the random walk considers all four possible directions
$\\{\text{North},\text{South},\text{East},\text{West}\\}$ and chooses a random
direction under the condition that it does not end up on a prohibited state or
outside of the state space. We want to calculate the probability that this
random walk of length $K$ finds the exit, located at $(N-1,N-1)$; i.e. it
visits this cell at least once within $K$ steps. Simulations below are for
$N=12$, and $K=100$. Writing $w=(w_{0},\ldots,w_{K-1})$ with
$w_{i}\in\mathcal{S}$ for the random walk, and
$d:\mathcal{S}\to\mathbb{N}_{0}$ for the step-wise distance function of a
given cell $s\in\mathcal{S}$ to the exit, the surrogate function used was
$g(w)=-\min_{i\in\\{0,\ldots,K-1\\}}\\{d(w_{i})\\}$. The method used for
finding random walks under the condition of “distance not exceeding a given
threshold” was rejection sampling (on the space of random walks).

Figure 6: A two-dimensional labyrinth with start $(0,0)$ (upper left corner)
and goal $(11,11)$ (lower right corner) and a randomly sampled random walk of
length $100$. Main plot below: Each point marks the final position of the
random walks sampled from Nested Sampling, with light colors being generated
later in Nested Sampling’s iteration. Top row: Two sample random walks, early
(left) and last (right) in Nested Sampling’s iteration. Figure 7: Estimates of
logprobability of “accidental labyrinth solve” via Nested Sampling, for $J=3$,
$J=9$, and $J=25$, with 25 separate runs each. Blue line marks a brute-force
estimate of the logprobability obtained by brute force Monte Carlo sampling
one million random walks.

##  5 Conclusions and outlook

We have studied the generalisation of Nested Sampling towards general
quadrature problems, including the plateau case, and have provided both
mathematical verification of correctness and computational improvement via the
surrogate function methodology. A quadrature problem that we have focused on
particularly is that of estimating rare event probabilities. Here, we have
shown efficiency and effectiveness of our estimator in a multitude of very
heterogeneous numerical experiments.

A rigorous proof of Section 2.1, as well as an accurate estimation of the
quadrature error produced by the Nested Sampling procedure are obvious
directions for future research.

##  Appendix A Facts about generalized inverses

The following lemma can be found, e.g., in [10, proposition 2.(1)], but is
recorded here for completeness of exposition.

######  Lemma A.1.

If $\mathcal{L}:\Omega\to\mathbb{R}$ is a measureable map without plateaus of
nonvanishing measure $\mu$, and $X$ is its survival function
$X(\lambda)=\mu(\mathcal{L}>\lambda)$, then the push-forward of $\mu$ along
$X\circ\mathcal{L}$ (or alternatively, the law of $X\circ\mathcal{L}$,
interpreted as a random variable) is the uniform measure on $[0,1]$

######  Proof A.2.

This is easiest seen by looking at the cumulative distribution function of
$\mathcal{F}_{\mathcal{L}}\circ\mathcal{L}$, where
$\mathcal{F}_{\mathcal{L}}=1-X$ is the cdf of $\mathcal{L}$. We note that
$X:\mathbb{R}\to(0,1)$ is invertible since $\mathcal{L}$ has no plateaus

$\displaystyle\mu(\mathcal{F}_{\mathcal{L}}\circ\mathcal{L}\leq\alpha)$
$\displaystyle=\mu(\mathcal{L}\leq\mathcal{F}_{\mathcal{L}}^{-1}(\alpha))=\mathcal{F}_{\mathcal{L}}(\mathcal{F}_{\mathcal{L}}^{-1}(\alpha))=\alpha$

This shows that the cdf of $\mathcal{F}_{\mathcal{L}}\circ\mathcal{L}$ is
indeed the cdf of a uniform distribution on $(0,1)$. Since
$X=1-\mathcal{F}_{\mathcal{L}}$, the same holds for $X$.

The following statements are from [41].

######  Definition A.3 (Generalized inverse).

Let $T:\mathbb{R}\to\mathbb{R}$ be a non-decreasing101010i.e. $x<y$ implies
$T(x)\leq T(y)$ function where we set $T(-\infty)=\lim_{x\to-\infty}T(x)$ and
$T(\infty)=\lim_{x\to\infty}T(x)$. Then the left-continuous generalized
inverse $T^{-}:\mathbb{R}\to\bar{\mathbb{R}}$ of $T$ is defined by

(10) $\displaystyle T^{-}(y)$ $\displaystyle=\inf\\{x\in\mathbb{R}:T(x)\geq
y\\}$

with the convention that $\inf\emptyset=\infty$.

We follow up with a list of elementary properties of $T^{+}$. This is an
adaptation of [10, Proposition 1] to our version of generalized inverse and
similar to [7, Proposition 4.2] (but with some errors fixed).

######  Lemma A.4.

Let $T:\mathbb{R}\to\mathbb{R}$ be a nondecreasing map.

  1. (a)

$y\leq T(x)$ if and only if $T^{-}(y)\leq x$.

  2. (b)

We define the left-continuous and right-continuous versions $T_{l}(x):=T(x-)$
and $T_{r}(x):=T(x+)$ of $T$. Then $T_{l}^{+}=T_{r}^{+}$ as well as
$T_{l}^{-}=T_{r}^{-}$.

######  Lemma A.5.

Let $T$ be nondecreasing and continuous from the left. We denote by
$X=\\{x_{i}\\}$ the (ordered) list of all discontinuities of $T$, i.e.
$y_{i}^{+}:=T(x_{i}+)>T(x_{i})=:y_{i}^{-}$ and $T(x+)=T(x)$ for $x\not\in X$.
We denote by $Y=\\{y_{i}\\}$ the (ordered) list of plateaus of $T$, i.e. for
each $y_{i}$ there exists a proper (maximal in the set of half-open intervals)
interval $I_{i}=(x_{i}^{-},x_{i}^{+}]$ such that $T(x)\equiv y_{i}$ for all
$x\in I_{i}$. Then

$\displaystyle T(T^{-}(y))$ $\displaystyle=\begin{cases}y_{i}^{-},&\text{ for
}y\in(y_{i}^{-},y_{i}^{+}]\\\ y,&\text{ else }\end{cases}$ $\displaystyle
T^{-}(T(x))$ $\displaystyle=\begin{cases}x_{i}^{-},&\text{ for
}x\in(x_{i}^{-},x_{i}^{+}]\\\ x,&\text{ else }\end{cases}$

## Acknowledgments

D.S and P.W. want to thank Maria Neuss-Radu for fruitful discussion.

## References

  * [1] S. Aitken and O. E. Akman, Nested sampling for parameter inference in systems biology: application to an exemplar circadian model, BMC Systems Biology, 7 (2013), p. 72, https://doi.org/10.1186/1752-0509-7-72.
  * [2] S.-K. Au and J. L. Beck, Estimation of small failure probabilities in high dimensions by subset simulation, Probabilistic Engineering Mechanics, 16 (2001), pp. 263–277, https://doi.org/https://doi.org/10.1016/S0266-8920(01)00019-4, https://www.sciencedirect.com/science/article/pii/S0266892001000194.
  * [3] R. J. Baldock, N. Bernstein, K. M. Salerno, L. B. Pártay, and G. Csányi, Constant-pressure nested sampling with atomistic dynamics, Physical Review E, 96 (2017), p. 043311.
  * [4] N. Björsell and P. Händel, Truncated gaussian noise in adc histogram tests, Measurement, 40 (2007), pp. 36–42.
  * [5] B. Burrows, A new approach to numerical integration, IMA J. Appl. Math., 26 (1980), pp. 151–173.
  * [6] N. Chopin and C. P. Robert, Properties of nested sampling, Biometrika, 97 (2010), pp. 741–755.
  * [7] A. de La Fortelle, A study on generalized inverses and increasing functions part i: generalized inverses, (2015).
  * [8] A. d’Onofrio, Bounded noises in physics, biology, and engineering, Springer, 2013.
  * [9] R. Dybowski, T. J. McKinley, P. Mastroeni, and O. Restif, Nested sampling for bayesian model comparison in the context of salmonella disease dynamics, PLOS ONE, 8 (2013), pp. 1–17, https://doi.org/10.1371/journal.pone.0082317.
  * [10] P. Embrechts and M. Hofert, A note on generalized inverses, Mathematical Methods of Operations Research, 77 (2013), pp. 423–432.
  * [11] M. Evans, Discussion of nested sampling for bayesian computations by john skilling, Bayesian Statistics, 8 (2007), pp. 491–524.
  * [12] F. Feroz, M. Hobson, E. Cameron, and A. Pettitt, Importance nested sampling and the multinest algorithm, arXiv preprint arXiv:1306.2144, (2013).
  * [13] F. Feroz and M. P. Hobson, Multimodal nested sampling: an efficient and robust alternative to markov chain monte carlo methods for astronomical data analyses, Mon. Not. R. Astron. Soc., 384 (2008), pp. 449–463, https://doi.org/10.1111/j.1365-2966.2007.12353.x.
  * [14] F. Feroz, M. P. Hobson, and M. Bridges, Multinest: an efficient and robust bayesian inference tool for cosmology and particle physics, Monthly Notices of the Royal Astronomical Society, 398 (2009), pp. 1601–1614, https://doi.org/10.1111/j.1365-2966.2009.14548.x.
  * [15] A. Fowlie, W. Handley, and L. Su, Nested sampling with plateaus, Monthly Notices of the Royal Astronomical Society, 503 (2021), pp. 1199–1205.
  * [16] A. Fowlie, S. Hoof, and W. Handley, Nested sampling for frequentist computation: fast estimation of small p-values, Physical Review Letters, 128 (2022), p. 021801.
  * [17] P. Glewwe, Measurement error bias in estimates of income and income growth among the poor: analytical results and a correction formula, Economic Development and Cultural Change, 56 (2007), pp. 163–189.
  * [18] E. Higson, W. Handley, M. Hobson, and A. Lasenby, Dynamic nested sampling: an improved algorithm for parameter estimation and evidence calculation, Stat. Comput., 29 (2019), pp. 891–913.
  * [19] E. Higson, W. Handley, M. Hobson, A. Lasenby, et al., Sampling errors in nested sampling parameter estimation, Bayesian Anal., 13 (2018), pp. 873–896.
  * [20] P. Kasperkovitz and C. Tutschka, Soft chaos in a hamiltonian system with step potential. i: Statistical properties, Chaos, Solitons & Fractals, 24 (2005), pp. 839–861, https://doi.org/https://doi.org/10.1016/j.chaos.2004.11.001, https://www.sciencedirect.com/science/article/pii/S0960077904006861.
  * [21] C. R. Keeton, On statistical uncertainty in nested sampling, Monthly Notices of the Royal Astronomical Society, 414 (2011), pp. 1418–1426.
  * [22] J. Latz, On the Well-posedness of Bayesian Inverse Problems, SIAM/ASA Journal on Uncertainty Quantification, 8 (2020), pp. 451–482, https://doi.org/10.1137/19M1247176, https://doi.org/10.1137/19M1247176.
  * [23] D. J. MacKay, Information theory, inference and learning algorithms, Cambridge university press, 2003.
  * [24] C. W. Miller and I. Yang, Optimal control of conditional value-at-risk in continuous time, SIAM Journal on Control and Optimization, 55 (2017), pp. 856–884, https://doi.org/10.1137/16M1058492, https://doi.org/10.1137/16M1058492, https://arxiv.org/abs/https://doi.org/10.1137/16M1058492.
  * [25] P. Mukherjee, D. Parkinson, and A. R. Liddle, A nested sampling algorithm for cosmological model selection, The Astrophysical Journal, 638 (2006), pp. L51–L54, https://doi.org/10.1086/501068.
  * [26] I. Murray, D. MacKay, Z. Ghahramani, and J. Skilling, Nested sampling for potts models, in Advances in Neural Information Processing Systems, 2006, pp. 947–954.
  * [27] I. A. Murray, Advances in Markov chain Monte Carlo methods, PhD thesis, University of London, 2007.
  * [28] I. Papaioannou, C. Papadimitriou, and D. Straub, Sequential importance sampling for structural reliability analysis, Structural Safety, 62 (2016), pp. 66–75, https://doi.org/https://doi.org/10.1016/j.strusafe.2016.06.002, https://www.sciencedirect.com/science/article/pii/S0167473016300169.
  * [29] D. Parkinson and A. R. Liddle, Bayesian model averaging in astrophysics: a review, Stat. Anal. Data Min., 6 (2013), pp. 3––14, https://doi.org/10.1002/sam.11179, http://dx.doi.org/10.1002/sam.11179.
  * [30] L. B. Pártay, A. P. Bartók, and G. Csányi, Nested sampling for materials: The case of hard spheres, Physical Review E, 89 (2014), p. 022302.
  * [31] N. Pullen and R. J. Morris, Bayesian model comparison and parameter inference in systems biology using nested sampling, PLOS ONE, 9 (2014), pp. 1–11, https://doi.org/10.1371/journal.pone.0088419.
  * [32] R. Rackwitz and B. Flessler, Structural reliability under combined random load sequences, Computers & Structures, 9 (1978), pp. 489–494, https://api.semanticscholar.org/CorpusID:119496005.
  * [33] R. Y. Rubinstein, Optimization of computer simulation models with rare events, European Journal of Operational Research, 99 (1997), pp. 89–112, https://doi.org/https://doi.org/10.1016/S0377-2217(96)00385-2, https://www.sciencedirect.com/science/article/pii/S0377221796003852.
  * [34] R. Salomone, L. F. South, C. C. Drovandi, and D. P. Kroese, Unbiased and consistent nested sampling via sequential monte carlo, arXiv preprint arXiv:1805.03924, (2018).
  * [35] D. S. Sivia and J. Skilling, Data Analysis: A Bayesian Tutorial, Oxford University Press, 2nd ed., 2006.
  * [36] J. Skilling et al., Nested sampling for general bayesian computation, Bayesian Anal., 1 (2006), pp. 833–859.
  * [37] A. Sripad and D. Snyder, A necessary and sufficient condition for quantization errors to be uniform and white, IEEE Transactions on Acoustics, Speech, and Signal Processing, 25 (1977), pp. 442–448.
  * [38] L. N. Trefethen, Exactness of Quadrature Formulas, SIAM Review, 64 (2022), pp. 132–150, https://doi.org/10.1137/20M1389522, https://doi.org/10.1137/20M1389522, https://arxiv.org/abs/https://doi.org/10.1137/20M1389522.
  * [39] S. Vegetti and L. V. Koopmans, Bayesian strong gravitational-lens modelling on adaptive grids: objective detection of mass substructure in galaxies, Monthly Notices of the Royal Astronomical Society, 392 (2009), pp. 945–963.
  * [40] J. Veitch, V. Raymond, B. Farr, W. Farr, P. Graff, S. Vitale, B. Aylott, K. Blackburn, N. Christensen, M. Coughlin, et al., Parameter estimation for compact binaries with ground-based gravitational-wave observations using the lalinference software library, Physical Review D, 91 (2015), p. 042003.
  * [41] P. Wacker, Please, not another note about generalized inverses, arXiv preprint arXiv:2306.06989, (2023).
  * [42] M. R. Walker and M. Vucelja, Anomalous thermal relaxation of langevin particles in a piecewise-constant potential, Journal of Statistical Mechanics: Theory and Experiment, 2021 (2021), p. 113105.
  * [43] C. Walter, Point process-based monte carlo estimation, Statistics and Computing, 27 (2017), pp. 219–236.
  * [44] Y. Yao, H. Deng, M. Yi, and J. Ma, Impact of bounded noise on the formation and instability of spiral wave in a 2d lattice of neurons, Scientific Reports, 7 (2017), p. 43151, https://doi.org/10.1038/srep43151.