###### Abstract

The Shannon entropy, and related quantities such as mutual information, can be
used to quantify uncertainty and relevance. However, in practice, it can be
difficult to compute these quantities for arbitrary probability distributions,
particularly if the probability mass functions or densities cannot be
evaluated. This paper introduces a computational approach, based on Nested
Sampling, to evaluate entropies of probability distributions that can only be
sampled. I demonstrate the method on three examples: a simple gaussian example
where the key quantities are available analytically; (ii) an experimental
design example about scheduling observations in order to measure the period of
an oscillating signal; and (iii) predicting the future from the past in a
heavy-tailed scenario.

###### keywords:

information theory; entropy; mutual information; monte carlo; nested sampling;
bayesian inference

To be determined Accepted 08/2017 Computing Entropies With Nested Sampling
Brendon J. Brewer<EMAIL_ADDRESS>
##  1 Introduction

If an unknown quantity $x$ has a discrete probability distribution $p(x)$, the
Shannon entropy (Shannon, 1948; Cover and Thomas, 2012) is defined as

$\displaystyle H(x)$ $\displaystyle=-\sum_{x}p(x)\log p(x)$ (1)

where the sum is over all of the possible values of $x$ under consideration.
The entropy quantifies the degree to which the issue “what is the value of
$x$, precisely?” remains unresolved (Knuth, 2005). More generally, the
‘Kullback-Leibler divergence’ quantifies the degree to which a probability
distribution $p(x)$ departs from a ‘base’ or ‘prior’ probability distribution
$q(x)$:

$\displaystyle D_{\rm KL}\big{(}p\,||\,q\big{)}$
$\displaystyle=\sum_{x}p(x)\log\left[\frac{p(x)}{q(x)}\right]$ (2)

(Knuth and Skilling, 2012; Caticha and Giffin, 2006). If $q(x)$ is taken to be
uniform, this is equivalent (up to an additive constant and a change of sign)
to the Shannon entropy $H(x)$. With the negative sign reinstated, Kullback-
Leibler divergence is sometimes called relative entropy, and is more
fundamental than the Shannon entropy (Knuth and Skilling, 2012). However, this
paper focuses on the Shannon entropy for simplicity, and because the resulting
algorithm calculates the expected value of a log-probability.

The Shannon entropy can be interpreted straightforwardly as a measure of
uncertainty. For example, if there is only one possible value which has
probability one, $H(x)$ is zero. Conventionally, $0\times\log 0$ is defined to
be equal to $\lim_{x\to 0^{+}}(x\log x)$, which is zero; i.e., ‘possibilities’
with probability zero do not contribute to the sum. If there are $N$
possibilities with equal probabilities $1/N$ each, then the entropy is
$H(x)=\log(N)$. If the space of possible $x$-values is continuous so that
$p(x)$ is a probability density function, the differential entropy

$\displaystyle H(x)$ $\displaystyle=-\int p(x)\log p(x)\,dx$ (3)

quantifies uncertainty by generalising the log-volume of the plausible region,
where volume is defined with respect to $dx$. Crucially, in continuous spaces,
the differential entropy is not invariant under changes of coordinates, so any
value given for a differential entropy must be understood with respect to the
coordinate system that was used. Even in discrete cases, using the Shannon
entropy assumes that each possibility contributes equally when counting
possibilities — if this is inappropriate, the relative entropy should be used.

Entropies, relative entropies, and related quantities, tend to be analytically
available only for a few families of probability distributions. On the
numerical side, if $\log p(x)$ can be evaluated for any $x$, then simple Monte
Carlo will suffice for approximating $H(x)$. On the other hand, if $p(x)$ can
only be sampled but not evaluated (for example, if $p(x)$ is a marginal or
conditional distribution, or the distribution of a quantity derived from $x$),
then this will not work. Kernel density estimation or the like (e.g. Szabó,
2014) may be effective in this case, but is unlikely to generalise well to
high dimensions.

This paper introduces a computational approach to evaluating the Shannon
entropy $H(x)$ of probability distributions that can only be sampled using
Markov Chain Monte Carlo (MCMC). The key idea is that Nested Sampling can
create an unbiased estimate of a log-probability, and entropies are averages
of log-probabilities.

###  1.1 Notation and conventions

Throughout this paper, I use the compact ‘overloaded’ notation for probability
distributions favoured by many Bayesian writers (Jaynes, 2003; MacKay, 2003),
writing $p(x)$ for either a probability mass function or a probability density
function, instead of the more cumbersome $P(X=x)$ or $f_{X}(x)$. In the
compact notation, there is no distinction between the ‘random variable’ itself
($X$) and a similarly-named dummy variable ($x$). Probability distributions
are implicitly conditional on some prior information, which is omitted from
the notation unless necessary. All logarithms are written $\log$ and any base
can be used, unless otherwise specified (for example by writing $\ln$ or
$\log_{10}$). Any numerical values given for the value of specific entropies
are in nats, i.e., the natural logarithm was used.

Even though the entropy is written as $H(x)$, it is imperative that we
remember it is not a property of the value of $x$ itself, but a property of
the probability distribution used to describe a state of ignorance about $x$.
Throughout this paper, $H(x)$ is used as notation for both Shannon entropies
(in discrete cases) and differential entropies (in continuous cases). Which
one it is should be clear from the context of the problem at hand. I sometimes
write general formulae in terms of sums (i.e., as they would appear in a
discrete problem), and sometimes as integrals, as they would appear in a
continuous problem.

##  2 Entropies in Bayesian inference

Bayesian inference is the use of probability theory to describe uncertainty,
often about unknown quantities (‘parameters’) $\boldsymbol{\theta}$. Some data
$\boldsymbol{d}$, initially unknown but thought to be relevant to
$\boldsymbol{\theta}$, is obtained. The prior information leads the user to
specify a prior distribution $p(\boldsymbol{\theta})$ for the unknown
parameters, along with a conditional distribution
$p(\boldsymbol{d}|\boldsymbol{\theta})$ describing prior knowledge of how the
data is related to the parameters (i.e., if the parameters were known, what
data would likely be observed?). By the product rule, this yields a joint
prior

$\displaystyle p(\boldsymbol{\theta},\boldsymbol{d})$
$\displaystyle=p(\boldsymbol{\theta})p(\boldsymbol{d}|\boldsymbol{\theta})$
(4)

which is the starting point for Bayesian inference (Caticha, 2008; Caticha and
Giffin, 2006). In practical applications, the following operations are usually
feasible and have a low computational cost:

  1. 1.

Samples can be generated from the prior $p(\boldsymbol{\theta})$;

  2. 2.

Simulated datasets can be generated from
$p(\boldsymbol{d}|\boldsymbol{\theta})$ for any given value of
$\boldsymbol{\theta}$;

  3. 3.

The likelihood, $p(\boldsymbol{d}|\boldsymbol{\theta})$, can be evaluated
cheaply for any $\boldsymbol{d}$ and $\boldsymbol{\theta}$. Usually it is the
log-likelihood that is actually implemented, for numerical reasons.

Throughout this paper I assume these operations are available and inexpensive.

###  2.1 The relevance of data

The entropy of the prior $p(\boldsymbol{\theta})$ describes the degree to
which the question “what is the value of $\boldsymbol{\theta}$, precisely?”
remains unanswered, while the entropy of the joint prior
$p(\boldsymbol{\theta},\boldsymbol{d})$ describes the degree to which the
question “what is the value of the pair
$(\boldsymbol{\theta},\boldsymbol{d})$, precisely?” remains unanswered. The
degree to which the question “what is the value of $\boldsymbol{\theta}$?”
would remain unresolved if $\boldsymbol{d}$ were resolved is given by the
conditional entropy

$\displaystyle H(\boldsymbol{\theta}|\boldsymbol{d})$
$\displaystyle=-\sum_{\boldsymbol{d}}p(\boldsymbol{d})\sum_{\boldsymbol{\theta}}p(\boldsymbol{\theta}|\boldsymbol{d})\log
p(\boldsymbol{\theta}|\boldsymbol{d})$ (5)

which is the expected value of the entropy of the posterior, averaged over all
possible datasets which might be observed.

One might wish to compute $H(\boldsymbol{\theta}|\boldsymbol{d})$, perhaps to
compare it to $H(\boldsymbol{\theta})$ and quantify how much might be learned
about $\boldsymbol{\theta}$. This would be difficult because the expression
for the posterior distribution

$\displaystyle p(\boldsymbol{\theta}|\boldsymbol{d})$
$\displaystyle=\frac{p(\boldsymbol{\theta})p(\boldsymbol{d}|\boldsymbol{\theta})}{p(\boldsymbol{d})}$
(6)

contains the marginal likelihood:

$\displaystyle p(\boldsymbol{d})$
$\displaystyle=\sum_{\boldsymbol{\theta}}p(\boldsymbol{\theta})p(\boldsymbol{d}|\boldsymbol{\theta})$
(7)

also known as the ‘evidence’, which tends to be computable but costly,
especially when the summation is replaced by an integration over a high-
dimensional parameter space.

It is important to distinguish between $H(\boldsymbol{\theta}|\boldsymbol{d})$
and the entropy of $\boldsymbol{\theta}$ given a particular value of
$\boldsymbol{d}$, which might be written
$H(\boldsymbol{\theta}|\boldsymbol{d}=\boldsymbol{d}_{\rm obs})$ (‘obs’ for
observed). The former measures the degree to which one question answers
another ex ante, and is a function of two questions. The latter measures the
degree to which a question remains unresolved after conditioning on a specific
statement (ex post), and is a function of a question and a statement.

###  2.2 Mutual information

The mutual information is another way of describing the relevance of the data
to the parameters. Its definition, and relation to other quantities, is

$\displaystyle I(\boldsymbol{\theta};\boldsymbol{d})$
$\displaystyle=\sum_{\boldsymbol{\theta}}\sum_{\boldsymbol{d}}p(\boldsymbol{\theta},\boldsymbol{d})\log\left[\frac{p(\boldsymbol{\theta},\boldsymbol{d})}{p(\boldsymbol{\theta})p(\boldsymbol{d})}\right]$
(8)
$\displaystyle=H(\boldsymbol{\theta})+H(\boldsymbol{d})-H(\boldsymbol{\theta},\boldsymbol{d})$
(9) $\displaystyle=H(\boldsymbol{d})-H(\boldsymbol{d}|\boldsymbol{\theta})$
(10)
$\displaystyle=H(\boldsymbol{\theta})-H(\boldsymbol{\theta}|\boldsymbol{d}).$
(11)

The mutual information can also be written as the expected value (with respect
to the prior over datasets $p(\boldsymbol{d})$) of the Kullback-Leibler
divergence from prior to posterior:

$\displaystyle I(\boldsymbol{\theta};\boldsymbol{d})$
$\displaystyle=\sum_{\boldsymbol{d}}p(\boldsymbol{d})D_{\rm
KL}\big{(}p(\boldsymbol{\theta}|\boldsymbol{d})\,||\,p(\boldsymbol{\theta})\big{)}.$
(12)

In terms of the prior, likelihood, and evidence, it is

$\displaystyle I(\boldsymbol{\theta};\boldsymbol{d})$
$\displaystyle=\sum_{\boldsymbol{\theta}}\sum_{\boldsymbol{d}}p(\boldsymbol{\theta})p(\boldsymbol{d}|\boldsymbol{\theta})\left[\log
p(\boldsymbol{d}|\boldsymbol{\theta})-\log p(\boldsymbol{d})\right],$ (13)

i.e., the mutual information is the prior expected value of the log likelihood
minus the log evidence. As with the conditional entropy, the computational
difficulty appears in the form of the log evidence, $\log
p(\boldsymbol{d})=\log\left(\sum_{\boldsymbol{\theta}}p(\boldsymbol{\theta})p(\boldsymbol{d}|\boldsymbol{\theta})\right)$,
which must be evaluated or estimated for many possible datasets.

For experimental design purposes, maximising the expected amount of
information obtained from the data is a sensible goal. Formally, either
maximising $I(\boldsymbol{\theta};\boldsymbol{d})$ or minimising
$H(\boldsymbol{\theta}|\boldsymbol{d})$ will produce the same result because
the prior $p(\boldsymbol{\theta})$ does not vary with the experimental design.
Reference priors (Bernardo, 2005) also maximise
$I(\boldsymbol{\theta};\boldsymbol{d})$ but vary the prior
$p(\boldsymbol{\theta})$ in the maximisation process while keeping the
experimental design fixed.

##  3 Nested Sampling

Nested Sampling (NS), introduced by Skilling (2006), is an algorithm whose aim
is to calculate the evidence

$\displaystyle p(\boldsymbol{\theta})$ $\displaystyle=\int
p(\boldsymbol{\theta})p(\boldsymbol{d}|\boldsymbol{\theta})\,d\boldsymbol{\theta},$
(14)

or, in the simplified notation common when discussing computational matters,

$\displaystyle Z$
$\displaystyle=\int\pi(\boldsymbol{\theta})L(\boldsymbol{\theta})\,d\boldsymbol{\theta},$
(15)

where $\boldsymbol{\theta}$ is the unknown parameter(s), $\pi$ is the prior
distribution, and $L$ is the likelihood function. The original NS algorithm,
and variations within the same family (Feroz et al., 2009; Brewer et al.,
2011; Handley et al., 2015), have become popular tools in Bayesian data
analysis (Knuth et al., 2015; Pullen and Morris, 2014; Brewer and Donovan,
2015) and have also found use in statistical physics (Pártay et al., 2010;
Baldock et al., 2016; Martiniani et al., 2014), which was another of its aims.
NS also estimates the Kullback-Leibler divergence from the prior
$\pi(\boldsymbol{\theta})$ to the posterior $\pi(x)L(x)/Z$, which Skilling
calls the information. This is a measure of how much was learned about
$\boldsymbol{\theta}$ from the specific dataset, and is also useful for
defining a termination criterion for NS.

While $Z$ is a simple expectation with respect to $\pi$, the distribution of
$L$-values implied by $\pi$ tends to be very heavy-tailed, which is why simple
Monte Carlo does not work. Equivalently, the integral over
$\boldsymbol{\theta}$ is dominated by very small regions where
$L(\boldsymbol{\theta})$ is high.

To overcome this, NS evolves a population of $N$ particles in the parameter
space towards higher likelihood regions. The particles are initially
initialised from the prior $\pi(\boldsymbol{\theta})$, and the particle with
the lowest likelihood, $L^{*}_{1}$, is found and its details are recorded as
output. This worst particle is then discarded and replaced by a new particle
drawn from the prior $\pi$ but subject to the constraint that its likelihood
must be above $L^{*}_{1}$. This is usually achieved by cloning a surviving
particle and evolving it with MCMC to explore the prior distribution but with
a constraint on likelihood value, $L(\boldsymbol{\theta})>L^{*}_{1}$. This
‘constrained prior’ distribution has probability density proportional to

$\displaystyle\pi(\boldsymbol{\theta})\mathds{1}\left(L(\boldsymbol{\theta})>L^{*}_{1}\right)$
(16)

where $\mathds{1}()$ is an indicator function which returns one if the
argument is true and zero if it is false.

This process is repeated, resulting in an increasing sequence of likelihoods

$\displaystyle L^{*}_{1},L^{*}_{2},L^{*}_{3},...$ (17)

Defining $X(\ell)$ as the amount of prior mass with likelihood greater than
some threshold $\ell$,

$\displaystyle X(\ell)$
$\displaystyle=\int\pi(\boldsymbol{\theta})\mathds{1}\left(L(\boldsymbol{\theta})>\ell\right)\,d\boldsymbol{\theta},$
(18)

each discarded point in the sequence can be assigned an $X$-value, which
transforms $Z$ to a one-dimensional integral. Skilling (2006)’s key insight
was that the $X$-values of the discarded points are approximately known from
the nature of the algorithm. Specifically, each iteration reduces the prior
mass by approximately a factor $e^{-1/N}$. More accurately, the conditional
distribution of $X_{i}$ given $X_{i-1}$ (and the definition of the algorithm)
is obtainable by letting $t_{i}\sim\textnormal{Beta}(N,1)$ and setting
$X_{i+1}:=t_{i}X_{i-1}$ (where $X_{1}=t_{1}$). The $t_{i}$ quantities are the
proportion of the remaining prior mass that is retained at each iteration,
after imposing the constraint of the latest $L^{*}$ value.

###  3.1 The sequence of $X$ values

The sequence of $X$ values is defined by

$\displaystyle X_{1}$ $\displaystyle=t_{1}$ (19) $\displaystyle X_{2}$
$\displaystyle=t_{2}X_{1}$ (20) $\displaystyle X_{3}$
$\displaystyle=t_{3}X_{2}$ (21)

and so on, i.e. $X_{i}$ can be written as

$\displaystyle X_{i}=\prod_{k=1}^{i}t_{k}.$ (22)

Taking the log and multiplying both sides by -1, we get

$\displaystyle-\ln X_{i}=\sum_{k=1}^{i}\left(-\ln t_{k}\right).$ (23)

Consider $t_{i}$, the fraction of the remaining prior mass that is retained in
each NS iteration. By a change of variables, the distribution for these
compression factors $t_{i}\sim$Beta$(N,1)$ corresponds to an exponential
distribution for $-\ln(t_{i})$:

$\displaystyle-\ln t_{i}$ $\displaystyle\sim\textnormal{Exponential}(N).$ (24)

From Equations 23 and 24, we see that the sequence of $-\ln X$ values produced
by a Nested Sampling run can be considered as a Poisson process with rate
equal to $N$, the number of particles. This is why separate NS runs can be
simply merged (Skilling, 2006; Henderson et al., 2017). Walter (2017) showed
how this view of the NS sequence of points can be used to construct a version
of NS that produces unbiased (in the frequentist sense) estimates of the
evidence. However, it can also be used to construct an unbiased estimator of a
log-probability, which is more relevant to information theoretic quantities
discussed in the present paper.

Consider a particular likelihood value $\ell$ whose corresponding $X$-value we
would like to know. Since the $-\ln X_{i}$ values have the same distribution
as the arrival times of a Poisson process with rate $N$, the probability
distribution for the number of points in an interval of length $w$ is Poisson
with expected value $Nw$. Therefore the expected number of NS discarded points
with likelihood below $\ell$ is $-N\ln X(\ell)$:

$\displaystyle\big{<}n(L(\boldsymbol{\theta}_{i})\leq\ell)\big{>}$
$\displaystyle=-N\ln X(\ell).$ (25)

We can therefore take $n(L(\boldsymbol{\theta}_{i})\leq\ell)/N$, the number of
points in the sequence with likelihood below $\ell$ divided by the number of
NS particles, as an unbiased estimator of $-\ln X(\ell)$. The standard
deviation of the estimator is $\sqrt{-\ln X(\ell)/N}$.

##  4 The algorithm

The above insight, that the number of NS discarded points with likelihood
below $\ell$ has expected value $-N\ln X(\ell)$, is the basis of the
algorithm. If we wish to measure a log-probability $-\ln X(\ell)$, we can use
NS to do it. $\pi(\boldsymbol{\theta})$ and $L(\boldsymbol{\theta})$ need not
be the prior and likelihood respectively, but can be any probability
distribution (over any space) and any function whose expected value is needed
for a particular application.

See Algorithm 1 for a step-by-step description of the algorithm. The algorithm
is written in terms of an unknown quantity $\boldsymbol{\theta}$ whose entropy
is required. In applications, $\boldsymbol{\theta}$ may be parameters, a
subset of the parameters, data, some function of the data, parameters and data
jointly, or whatever.

The idea is to generate a reference particle $\boldsymbol{\theta}_{\rm ref}$
from the distribution $p(\boldsymbol{\theta})$ whose entropy is required.
Then, a Nested Sampling run evolves a set of particles
$\\{\boldsymbol{\theta}\\}_{i=1}^{N}$, initially representing
$p(\boldsymbol{\theta})$, towards $\boldsymbol{\theta}_{\rm ref}$ in order to
measure the log-probability near $\boldsymbol{\theta}_{\rm ref}$ (see Figure
1). Nearness is defined using a distance function
$d(\boldsymbol{\theta};\boldsymbol{\theta}_{\rm ref})$, and the number of NS
iterations taken to reduce the distance to below some threshold $r$ (for
‘radius’) provides an unbiased estimate of

$\displaystyle 0pt=-\log\left[P(d(\boldsymbol{\theta};\boldsymbol{\theta}_{\rm
ref})<r)\right]$ (26)

which I call the ‘depth’. E.g., if $N=10$ Nested Sampling particles are used,
and it takes 100 iterations to reduce the distance to below $r$, then the
depth is 10 nats.

If we actually want the differential entropy

$\displaystyle H(\boldsymbol{\theta})$ $\displaystyle=-\int
p(\boldsymbol{\theta})\log p(\boldsymbol{\theta})\,d\boldsymbol{\theta}$ (27)

we can use the fact that density equals mass divided by volume. Assume $r$ is
small, so that

$\displaystyle P(d(\boldsymbol{\theta};\boldsymbol{\theta}_{\rm
ref})<r\,|\,\boldsymbol{\theta}_{\rm ref})$ $\displaystyle\approx
p(\boldsymbol{\theta}_{\rm
ref})\int_{d(\boldsymbol{\theta};\boldsymbol{\theta}_{\rm
ref})<r}\,d\boldsymbol{\theta}$ (28) $\displaystyle=p(\boldsymbol{\theta}_{\rm
ref})\,V$ (29)

where $V$ is the volume of the region where
$d(\boldsymbol{\theta};\boldsymbol{\theta}_{\rm ref})<r$. Then
$H(\boldsymbol{\theta})=\left<0pt\right>_{p(\boldsymbol{\theta})}+\log V(r)$.
It is important to remember that the definition of volume here is in terms of
the coordinates $\boldsymbol{\theta}$, and that a change of coordinates (or a
change of the definition of the distance function) would result in a different
value of the differential entropy. This is equivalent to counting
possibilities or volumes with a different choice of base measure
$q(\boldsymbol{\theta})$.

Reducing $r$ towards zero ought to lead to a more accurate result because the
probability density at a point will be more accurately approximated by the
probability of a small region containing the point, divided by the volume of
that region. On the other hand, smaller values of $r$ will lead to more Nested
Sampling iterations required, and hence less accuracy in the depth estimates.
Thus, there is a an awkward trade-off involved in choosing the value of $r$,
which will need to be investigated more thoroughly in the future. However, the
algorithm as it currently stands ought to be useful for calculating the
entropies of any low-to-moderate dimensional functions of an underlying high-
dimensional space.

Set the numerical parameters:

$N\in\\{1,2,3,...\\}$ $\triangleright$ the number of Nested Sampling particles
to use

$r\geq 0$ $\triangleright$ the tolerance

$\widehat{\mathbf{h}}\leftarrow\\{\\}$ $\triangleright$ Initialise an empty
list of results

while more iterations desired do

$k\leftarrow 0$ $\triangleright$ Initialise counter

Generate $\boldsymbol{\theta}_{\rm ref}$ from $p(\boldsymbol{\theta})$
$\triangleright$ Generate a reference point

Generate $\left\\{\boldsymbol{\theta}_{i}\right\\}_{i=1}^{N}$ from
$p(\boldsymbol{\theta})$ $\triangleright$ Generate initial NS particles

Calculate $d_{i}\leftarrow d(\boldsymbol{\theta}_{i};\boldsymbol{\theta}_{\rm
ref})$ for all $i$ $\triangleright$ Calculate distance of each particle from
the reference point

$i^{*}\leftarrow\textnormal{argmin}(\lambda i\rightarrow d_{i})$
$\triangleright$ Find the worst particle (greatest distance)

$d_{\rm max}\leftarrow d_{i^{*}}$ $\triangleright$ Find the greatest distance

while $d_{\rm max}>r$ do

Replace $\boldsymbol{\theta}_{i^{*}}$ with $\boldsymbol{\theta}_{\rm new}$
from $p\left(\boldsymbol{\theta}\,|\,d(\boldsymbol{\theta})<d_{\rm
max}\right)$ $\triangleright$ Replace worst particle

Calculate $d_{i^{*}}\leftarrow d(\boldsymbol{\theta}_{\rm
new};\boldsymbol{\theta}_{\rm ref})$ $\triangleright$ Calculate distance of
new particle from reference point

$i^{*}\leftarrow\textnormal{argmin}(\lambda i\rightarrow d_{i})$
$\triangleright$ Find the worst particle (greatest distance)

$d_{\rm max}\leftarrow d_{i^{*}}$ $\triangleright$ Find the greatest distance

$k\leftarrow k+1$ $\triangleright$ Increment counter $k$

$\widehat{\mathbf{h}}\leftarrow\textnormal{append}\big{(}\widehat{\mathbf{h}},k/N\big{)}$
$\triangleright$ Append latest estimate to results

$\widehat{h_{\rm final}}\leftarrow\frac{1}{\rm
num\\_iterations}\sum\widehat{\mathbf{h}}$ $\triangleright$ Average results

Algorithm 1 The algorithm which estimates the expected value of the depth:
$-\int p(\boldsymbol{\theta}_{\rm ref})\int
p(\boldsymbol{\theta})\log\left[P(d(\boldsymbol{\theta};\boldsymbol{\theta}_{\rm
ref})<r\,|\,\boldsymbol{\theta}_{\rm
ref})\right]\,d\boldsymbol{\theta}\,d\boldsymbol{\theta}_{\rm ref}$,
that is, minus the expected value of the log-probability of a small region
near $\boldsymbol{\theta}_{\rm ref}$, which can be converted to an estimate of
an entropy or differential entropy. The part highlighted in green is standard
Nested Sampling with quasi-prior $p(\boldsymbol{\theta})$ and quasi-likelihood
given by minus a distance function
$d(\boldsymbol{\theta};\boldsymbol{\theta}_{\rm ref})$. The final result,
$\widehat{h_{\rm final}}$, is an estimate of the expected depth.  Figure 1: To
evaluate the log-probability (or density) of the blue probability distribution
near the red point, Nested Sampling can be used, with the blue distribution
playing the role of the “prior” in NS, and the Euclidean distance from the red
point (illustrated with red contours) playing the role of the negative log-
likelihood. Averaging over selections of the red point gives an estimate of
the entropy of the blue distribution.

If the distance function $d(\boldsymbol{\theta}_{\rm
ref},\boldsymbol{\theta})$ is chosen to be Euclidean, the constraint
$d(\boldsymbol{\theta}_{\rm ref},\boldsymbol{\theta})<r$ corresponds to a ball
in the space of possible $\boldsymbol{\theta}_{\rm ref}$ values. The volume of
such a ball, of radius $r$ in $n$ dimensions, is

$\displaystyle V$ $\displaystyle=\frac{(\pi r^{2})^{n/2}}{(\frac{n}{2})!}.$
(30)

In one dimension, if the threshold $r$ is not small, we might think we are
calculating the entropy associated with the question ‘what is $\theta$, to
within a tolerance of $\pm r$?’. This is not quite correct. See Appendix A for
a discussion of this type of question.

In the following sections I demonstrate the algorithm on three example
problems. For each of the examples, I used 1000 repetitions of NS (i.e., 1000
reference particles), and used 10 NS particles and 10,000 MCMC steps per NS
iteration.

##  5 Example 1: Entropy of a prior for the data

This example demonstrates how the algorithm may be used to determine the
entropy for a quantity whose distribution can only be sampled, in this case, a
prior distribution over datasets.

Consider the basic statistics problem of inferring a quantity $\mu$ from 100
observations $\boldsymbol{x}=\\{x_{1},...,x_{100}\\}$ whose probability
distribution (conditional on $\mu$) is

$\displaystyle p(x_{i}|\mu)$ $\displaystyle\sim\textnormal{Normal}(\mu,1).$
(31)

I assumed a Normal$(0,10^{2})$ prior for $\mu$. Before getting the data, the
prior information is described by the joint distribution

$\displaystyle p(\mu,\boldsymbol{x})$
$\displaystyle=p(\mu)p(\boldsymbol{x}|\mu).$ (32)

The prior over datasets is given by the marginal distribution

$\displaystyle p(\boldsymbol{x})$ $\displaystyle=\int
p(\mu,\boldsymbol{x})\,d\mu.$ (33)

This is the distribution whose differential entropy I seek in this section.
Its true value is available analytically in this case, by noting that the
posterior is available in closed form for any dataset:

$\displaystyle\mu|\boldsymbol{x}$
$\displaystyle\sim\textnormal{Normal}\left(\frac{1}{100.01}\sum_{i=1}^{100}x_{i},\left[\frac{10}{\sqrt{10001}}\right]^{2}\right).$
(34)

which enables the calculation of the mutual information and hence
$H(\boldsymbol{x})$. The true value of $H(\boldsymbol{x})$ is $146.499$ nats.
However, in most situations, the entropy of a marginal distribution such as
this is not available in closed form.

To use the algorithm, a distance function $d(;)$ must be defined, to quantify
how far the NS particles are from the reference particle. To calculate the
entropy of the data, I applied the algorithm in the joint space of possible
parameters and datasets, i.e., distribution of reference particles followed
Equation 32. For the distance function, which helps the NS particles approach
the reference particle, I used a simple Euclidean metric in the space of
datasets:

$\displaystyle d\left(\boldsymbol{\theta};\boldsymbol{\theta}_{\rm
ref}\right)$ $\displaystyle=\sqrt{\sum_{i=1}^{100}\left(x_{i}-x_{{\rm
ref},i}\right)^{2}}.$ (35)

Since the distance function refers only to the data $\boldsymbol{x}$ and not
the parameter $\mu$, the sampling is effectively done in the space of possible
datasets only — the parameter $\mu$ can be thought of as merely a latent
variable allowing $p(\boldsymbol{x})$ to be explored conveniently.

I ran the algorithm and computed the average depth using a tolerance of
$r=10^{-3}\sqrt{100}$, so that the RMS difference between points in the two
datasets was about $10^{-3}$. From the algorithm, the estimate of the expected
depth was

$\displaystyle\left<0pt_{\boldsymbol{x}}\right>$ $\displaystyle=698.25\pm
0.34.$ (36)

The uncertainty here is the standard error of the mean, i.e., the standard
deviation of the depth estimates divided by square root of the number of NS
repetitions (the number of reference particles whose depth was evaluated).

The log-volume of a 100-dimensional ball of radius $r=10^{-3}\sqrt{100}$ is
$-551.76$. Therefore, the differential entropy is estimated to be

$\displaystyle H(\boldsymbol{x})$ $\displaystyle\approx 146.49\pm 0.34.$ (37)

which is very close to the true value (perhaps suspiciously close, but this
was just a fluke).

To carry out the sampling for this problem, Metropolis-Hastings moves that
explore the joint prior distribution $p(\mu,\boldsymbol{x})$ had to be
implemented. Four basic kinds of proposals were used: i) those which change
$\mu$ while leaving $\boldsymbol{x}$ fixed; ii) those which change $\mu$ and
shift $\boldsymbol{x}$ correspondingly; iii) those which resample a subset of
the $x$s from $p(\boldsymbol{x}|\mu)$; and iv) those which move a single $x$
slightly.

##  6 Example 2: Measuring the period of an oscillating signal

In physics and astronomy, it is common to measure an oscillating signal at a
set of times $\\{t_{1},...,t_{n}\\}$, in order to infer the amplitude, period,
and phase of the signal. Here, I demonstrate the algorithm on a toy version of
Bayesian experimental design: at what times should the signal be observed in
order to obtain as much information as possible about the period? To answer
this, we need to be able to calculate the mutual information between the
unknown period and the data.

As the algorithm lets us calculate the entropy of any distribution which can
be sampled by MCMC, there are several options. The mutual information can be
written in various ways, such as the following.

$\displaystyle I(\boldsymbol{\theta};\boldsymbol{d})$
$\displaystyle=H(\boldsymbol{\theta})+H(\boldsymbol{d})-H(\boldsymbol{\theta},\boldsymbol{d})$
(38) $\displaystyle I(\boldsymbol{\theta};\boldsymbol{d})$
$\displaystyle=H(\boldsymbol{d})-H(\boldsymbol{d}|\boldsymbol{\theta})$ (39)
$\displaystyle I(\boldsymbol{\theta};\boldsymbol{d})$
$\displaystyle=H(\boldsymbol{\theta})-H(\boldsymbol{\theta}|\boldsymbol{d})$
(40) $\displaystyle I(\boldsymbol{\theta};\boldsymbol{d})$ $\displaystyle=\int
p(\boldsymbol{d})\,D_{\rm
KL}\left(p(\boldsymbol{\theta}|\boldsymbol{d})\,||\,p(\boldsymbol{\theta})\right)\,d\boldsymbol{d}$
(41) $\displaystyle I(\boldsymbol{\theta};\boldsymbol{d})$
$\displaystyle=D_{\rm
KL}\left(p(\boldsymbol{\theta},\boldsymbol{d})\,||\,p(\boldsymbol{\theta})p(\boldsymbol{d})\right).$
(42)

In a Bayesian problem where $\boldsymbol{\theta}$ is an unknown parameter and
$\boldsymbol{d}$ is data, the first and second formulas would be costly to
compute, because the high dimensional probability distribution for the dataset
would require a large number of NS iterations to compute $H(\boldsymbol{d})$
and $H(\boldsymbol{\theta},\boldsymbol{d})$. The fourth involves the Kullback-
Leibler divergence from the prior to the posterior, averaged over all possible
datasets. This is straightforward to approximate with standard Nested
Sampling, but the estimate for a given dataset may be biased. It is quite
plausible that the bias is negligible compared to the variance from averaging
over many datasets. This is probably true in many situations.

However, more significantly, the above method would not work if we want
$I(\boldsymbol{\theta};\boldsymbol{d})$ for a single parameter or a subset of
them, rather than for all of the parameters. The prior-to-posterior KL
divergence measured by standard NS is the amount learned about all parameters,
including nuisance parameters. Therefore, I adopted the third formula as the
best way to compute $I(\boldsymbol{\theta};\boldsymbol{d})$. In the sinusoidal
example, the parameter of interest is the period, whereas the model will
contain other parameters as well. Since the period will have its prior
specified explicitly, $H(\boldsymbol{\theta})$ will be available analytically.
I then use the algorithm to obtain $H(\boldsymbol{\theta}|\boldsymbol{d})$ as
follows:

$\displaystyle H(\boldsymbol{\theta}|\boldsymbol{d})$ $\displaystyle=\int
p(\boldsymbol{d}_{\rm
ref})H(\boldsymbol{\theta}|\boldsymbol{d}=\boldsymbol{d}_{\rm
ref})\,d\boldsymbol{d}_{\rm ref}.$ (43)

This is an expected value over datasets, of the entropy of the posterior given
each dataset. I ran the algorithm by generating reference particles from
$p(\boldsymbol{\theta},\boldsymbol{d})$, i.e., generating parameter values
from the prior, along with corresponding simulated datasets. Given each
simulated dataset, I used NS to compute an unbiased estimate of
$H(\boldsymbol{\theta}|\boldsymbol{d}=\boldsymbol{d}_{\rm ref})$, by evolving
a set of particles (initially from the posterior) towards the reference
particle. The role of the “prior” in the NS is actually taken by the
posterior, and the “likelihood” was minus a distance function defined in terms
of the parameter(s) of interest.

To generate the initial NS particles from the posterior, I used standard MCMC
(targeting the posterior), initialised at the $\boldsymbol{\theta}$-value that
produced the simulated dataset $\boldsymbol{d}_{\rm ref}$, which is a perfect
sample from the posterior.

###  6.1 Assumptions

I consider two possible observing strategies, both involving $n=101$
measurements. The ‘even’ strategy has observations at times

$\displaystyle t_{i}=\frac{i-1}{n-1}$ (44)

for $i\in\\{1,2,...,101\\}$, that is, the observation times are
$0,0.01,0.02,...,0.99,1$ — evenly spaced from $t=0$ to $t=1$, including the
endpoints. The second, ‘uneven’ observing strategy schedules the observations
according to

$\displaystyle t_{i}=\left(\frac{i-\frac{1}{2}}{n}\right)^{3}$ (45)

which schedules observations close together initially, and further apart as
time goes on.

A purely sinusoidal signal has the form

$\displaystyle y(t)$ $\displaystyle=A\sin\left(\frac{2\pi t}{T}+\phi\right)$
(46)

where $A$ is the amplitude, $T$ is the period, and $\phi$ is the phase.
Throughout this section, I parameterise the period by its logarithm,
$\tau=\log_{10}T$. I assumed the following priors:

$\displaystyle\ln A$
$\displaystyle\sim\textnormal{Normal}\left(0,0.1^{2}\right)$ (47)
$\displaystyle\tau$ $\displaystyle\sim\textnormal{Uniform}(-1,0)$ (48)
$\displaystyle\phi$ $\displaystyle\sim\textnormal{Uniform}(0,2\pi)$ (49)

and the following conditional prior for the data:

$\displaystyle Y_{i}|A,T,\phi$
$\displaystyle\sim\textnormal{Normal}\left(y(t_{i}),0.1^{2}\right).$ (50)

That is, the data is just the signal $y(t)$ observed at particular times
$\\{t_{1},...,t_{n}\\}$, with gaussian noise of standard deviation 0.1. The
amplitude of the sinusoid is very likely to be around 10 times the noise
level, and the period is between 0.1 and 1 times the duration of the data. An
example signal observed with the even and uneven observing schedules is shown
in Figure 2.

Figure 2: A signal with true parameters $A=1$, $\tau=-0.5$, and $\phi=0$,
observed with noise standard deviation 0.1 with the even (gold points) and
uneven (green points) observing strategies.

To carry out the sampling for this problem, Metropolis-Hastings moves that
explore the posterior distribution $p(A,\tau,\phi|\boldsymbol{Y})$ had to be
implemented. I used heavy-tailed proposals (as in Section 8.3 of Brewer and
Foreman-Mackey (2016)) which modified one parameter at a time.

###  6.2 Results

Letting $\tau=\log_{10}T$, and treating $A$ and $\phi$ are nuisance
parameters, I calculated the conditional entropy $H(\tau|\boldsymbol{Y})$
using the algorithm with a tolerance of $r=10^{-5}$. The distance function was
just the absolute difference between $\tau$ and $\tau_{\rm ref}$.

The result was

$\displaystyle\left<0pt\right>$ $\displaystyle=5.379\pm 0.038$ (51)

Converting to a differential entropy by adding $\ln(2r)=-10.820$, the
conditional entropy is

$\displaystyle H(\tau|\boldsymbol{Y})$ $\displaystyle=-5.441\pm
0.038\textnormal{ nats}.$ (52)

Since the Uniform(0,1) prior for $\tau$ has a differential entropy of zero,
the mutual information is
$I(\tau;\boldsymbol{Y})=H(\tau)-H(\tau|\boldsymbol{Y})=5.441\pm 0.038$ nats.

The results for the uneven observing schedule were

$\displaystyle\left<0pt\right>$ $\displaystyle=5.422\pm 0.038$ (53)
$\displaystyle H(\tau|\boldsymbol{Y})$ $\displaystyle=-5.398\pm 0.038$ (54)
$\displaystyle I(\tau;\boldsymbol{Y})$ $\displaystyle=5.398\pm
0.038\textnormal{ nats}.$ (55)

The difference in mutual informations between the two observation schedules is
trivial. The situation may be different if we had allowed for shorter periods
to be possible, as irregular observing schedules are known to reduce aliasing
and ambiguity in the inferred period (Bretthorst, 2001). When inferring the
period of an oscillating signal, multimodal posterior distributions for the
period are common (Gregory, 2005; Brewer and Donovan, 2015). The multimodality
of the posterior here raises an interesting issue. Is the question we really
want answered “what is the value of $T$ precisely?”, to which the mutual
information relates? Most practicing scientists would not feel particularly
informed to learn that the vast majority of possibilities had been ruled out,
if the posterior still consisted of several widely separated modes! Perhaps,
in some applications, a more appropriate question is “what is the value of $T$
to within $\pm$ 10%”, or something along these lines. See Appendix A for more
about this issue.

In a serious experimental design situation, the relevance is not the only
consideration (or the answer to every experimental design problem would be to
obtain more data without limit), but it is a very important one.

##  7 Example 3: Data with Pareto distribution

Consider a situation where the probability distribution for some data values
$\boldsymbol{x}=\\{x_{1},x_{2},...,x_{n}\\}$ is a Pareto distribution,
conditional on a minimum value $x_{\rm min}=1$ and a slope $\alpha$:

$\displaystyle p(\boldsymbol{x}|\alpha)$
$\displaystyle=\prod_{i=1}^{n}\frac{\alpha}{x_{i}^{\alpha+1}}$ (56)

where all of the $x_{i}$ are greater than $x_{\rm min}=1$. If $\alpha$ is
relatively low (close to zero) then this distribution is very heavy-tailed,
and as $\alpha$ increases the probability density for the $x$s concentrates
just above $x=1$.

With heavy-tailed distributions, it can be difficult to predict future
observations from past ones, a notion popularised by Taleb (2007). This occurs
for three main reasons. First, the observations may not provide much
information about the parameters (in this case, $\alpha$). Second, even
conditional on the parameters, the predictive distribution for future data can
be heavy-tailed, implying a lot of uncertainty. Finally, one may always doubt
the probability assignment on which such a calculation is based. However, the
data clearly provides some nonzero amount of information about the future.

Consider $n=100$ data values, partitioned into two halves, so
$\boldsymbol{y}=\\{x_{1},...,x_{50}\\}=\\{y_{1},...,y_{50}\\}$ is the first
half of the data and
$\boldsymbol{z}=\\{x_{51},...,x_{100}\\}=\\{z_{1},...,z_{50}\\}$ is the second
half. Let $y_{\rm tot}$ and $z_{\rm tot}$ be the totals of the two subsets of
the data:

$\displaystyle y_{\rm tot}$ $\displaystyle=\sum_{i=1}^{50}y_{i}$ (57)
$\displaystyle z_{\rm tot}$ $\displaystyle=\sum_{i=1}^{50}z_{i}.$ (58)

In this section I use the algorithm to answer the questions: i) How uncertain
is $y_{\rm tot}$?; ii) How uncertain is $z_{\rm tot}$?; and iii) How
informative is $y_{\rm tot}$ about $z_{\rm tot}$? These are quantified by
$H(y_{\rm tot})$, $H(z_{\rm tot})$, and $I(y_{\rm tot};z_{\rm tot})$
respectively. This can be done despite the fact that $\alpha$ is unknown.

The prior for $\alpha$ was chosen to be lognormal such that the prior median
of $\alpha$ is 1 and the prior standard deviation of $\ln\alpha$ is also 1.
Due to the heavy-tailed nature of the Pareto distribution, I ran the algorithm
in terms of $\ln y_{\rm tot}$ and $\ln z_{\rm tot}$. This affects the
differential entropies but not the mutual information (which is invariant
under changes of parameterisation). The joint distribution for $\ln y_{\rm
tot}$ and $\ln z_{\rm tot}$, whose properties I aimed to obtain, is shown in
Figure 3.

Figure 3: The joint distribution for the log-sum of the first half of the
‘Pareto data’ and the second half. The goal is to calculate the entropy of the
marginal distributions, the entropy of the joint distribution, and hence the
mutual information. The distribution is quite heavy tailed despite the
logarithms, and extends beyond the domain shown in this plot.

Three types of Metropolis-Hastings proposals were used for this example: i) A
proposal which changes $\alpha$ and modifies all of the $x$s accordingly; ii)
A proposal which resamples a subset of the $x$s from the Pareto distribution;
and iii) A proposal which moves a single $x$-value slightly. A fourth proposal
which moves $\alpha$ while keeping $\boldsymbol{x}$ fixed, i.e., a posterior-
sampling move, was not implemented.

To compute the entropy of $\ln y_{\rm tot}$, the distance function was just
the absolute value of the distance between the NS particles and the reference
particle in terms of $\ln y_{\rm tot}$. The result, using a tolerance of
$r=10^{-3}$, was:

$\displaystyle\left<0pt_{\ln y_{\rm tot}}\right>$ $\displaystyle=8.851\pm
0.057.$ (59)

Correcting for the volume factor yields a differential entropy estimate of
$H(\ln y_{\rm tot})=2.366\pm 0.057$ nats. By the symmetry of the problem
specification, $H(\ln z_{\rm tot})=H(\ln y_{\rm tot})$.

The joint entropy $H(\ln y_{\rm tot},\ln z_{\rm tot})$ was also estimated by
using a Euclidean distance function and a tolerance of $10^{-3}$. The result
was

$\displaystyle\left<0pt_{(\ln y_{\rm tot},\ln z_{\rm tot})}\right>$
$\displaystyle=16.03\pm 0.11.$ (60)

Applying the log-volume correction yields a differential entropy estimate of
$H(\ln y_{\rm tot},\ln z_{\rm tot})=3.40\pm 0.11$ nats. Therefore, the mutual
information is

$\displaystyle I(y_{\rm tot};z_{\rm tot})$ $\displaystyle=H(\ln y_{\rm
tot})+H(\ln z_{\rm tot})-H(\ln y_{\rm tot},\ln z_{\rm tot})$ (61)
$\displaystyle=1.130\pm 0.082\textnormal{ nats}.$ (62)

The variance of this mutual information estimate was reduced by using a common
sequence of reference points for the marginal and joint entropy runs. This is
an instance of the ‘common random numbers’
technique111http://demonstrations.wolfram.com/TheMethodOfCommonRandomNumbersAnExample/.
In short, for those reference particles with a large depth in terms of the
$\ln y_{\rm tot}$ distance function, the depth in terms of the $(\ln y_{\rm
tot},\ln z_{\rm tot})$ distance function was also large. Since the mutual
information involves a difference of the averages of these depths, it is more
accurate to compute the average difference (i.e. compute the two depths using
the same reference particle, to compute the difference between the two depths)
than the difference of the averages.

##  8 Computational cost, limitations, and conclusion

The computational resources needed to compute these quantities was quite
large, as Nested Sampling appears in the inner loop of the algorithm. However,
it is worth reflecting on the complexity of the calculation that has been
done.

In a Bayesian inference problem with parameters $\boldsymbol{\theta}$ and data
$\boldsymbol{d}$, the mutual information is

$\displaystyle I(\boldsymbol{\theta};\boldsymbol{d})$ $\displaystyle=\iint
p(\boldsymbol{\theta},\boldsymbol{d})\ln\left[\frac{p(\boldsymbol{\theta},\boldsymbol{d})}{p(\boldsymbol{\theta})p(\boldsymbol{d})}\right]\,d\boldsymbol{\theta}\,d\boldsymbol{d}$
(63)

and this measures the dependence between $\boldsymbol{\theta}$ and
$\boldsymbol{d}$. A Monte Carlo strategy to evaluate this is to sample from
$p(\boldsymbol{\theta},\boldsymbol{d})$ and average the value of the
logarithm, which involves computing the log-evidence

$\displaystyle\ln p(\boldsymbol{d})=\ln\left[\int
p(\boldsymbol{\theta})p(\boldsymbol{d}|\boldsymbol{\theta})\,d\boldsymbol{\theta}\right]$
(64)

for each possible data set. Computing $p(\boldsymbol{d})$ for even a single
data set has historically been considered difficult.

However, if we want the mutual information between the data and a subset of
the parameters (i.e., there are nuisance parameters we don’t care about),
things become even more tricky. Suppose we partition the parameters
$\boldsymbol{\theta}$ into important parameters $\boldsymbol{\phi}$ and
nuisance parameters $\boldsymbol{\eta}$, such that we want to calculate the
mutual information $I(\boldsymbol{\phi};\boldsymbol{d})$. This is given by

$\displaystyle I(\boldsymbol{\phi};\boldsymbol{d})$ $\displaystyle=\iint
p(\boldsymbol{\phi},\boldsymbol{d})\ln\left[\frac{p(\boldsymbol{\phi},\boldsymbol{d})}{p(\boldsymbol{\phi})p(\boldsymbol{d})}\right]\,d\boldsymbol{\phi}\,d\boldsymbol{d}$
(65)

which is equivalent to Equation 63 but with $\boldsymbol{\theta}$ replaced
throughout by $\boldsymbol{\phi}$. However, this is an even more difficult
calculation than before, as $p(\boldsymbol{d}|\boldsymbol{\phi})$, the
marginal likelihood function with the nuisance parameters integrated out, is
typically unavailable in closed form. If we were to marginalise out the
nuisance parameters $\boldsymbol{\eta}$ explicitly, this would give us
Equation 65 with every probability distribution written as an explicit
integral over $\boldsymbol{\eta}$:

$\displaystyle I(\boldsymbol{\phi};\boldsymbol{d})$
$\displaystyle=\iint\left(\int
p(\boldsymbol{\phi},\boldsymbol{\eta})p(\boldsymbol{d}|\boldsymbol{\phi},\boldsymbol{\eta})\,d\boldsymbol{\eta}\right)\ln\left[\frac{\int
p(\boldsymbol{\phi},\boldsymbol{\eta})p(\boldsymbol{d}|\boldsymbol{\phi},\boldsymbol{\eta})\,d\boldsymbol{\eta}}{p(\boldsymbol{\phi})\int
p(\boldsymbol{\phi},\boldsymbol{\eta})p(\boldsymbol{d}|\boldsymbol{\phi},\boldsymbol{\eta})\,d\boldsymbol{\phi}\,\boldsymbol{d}\boldsymbol{\eta}}\right]\,d\boldsymbol{\phi}\,d\boldsymbol{d}$
(66)

It should not be surprising that this is costly. Indeed, it seems possible to
me that any general algorithm for computing entropies will have the property
that NS (or a similar algorithm capable of measuring small volumes) appears in
an inner loop, since entropies are sums of terms, each of which relates to the
downset of a given statement — the downset of a statement $S$ being the set of
all statements that imply $S$ (Knuth, 2005).

The current implementation of the algorithm (see Appendix B) is in C++ and has
not been parallelised. Given the high computational costs, parallelisation
should be attempted. This could follow the method of (Henderson et al., 2017),
or alternatively, different threads could work on different reference
particles simultaneously.

The algorithm as presented will tend to become less accurate on high
dimensional problems, as the required depth of the runs will increase
correspondingly, since many NS iterations will be required to get a high
dimensional vector of quantities to be close to the reference vector. The
accuracy of the $\ln 0pt$ estimates should degrade in proportion to
$\sqrt{0pt}$, due to the Poisson-process nature of the $-\ln X$ sequence.

Throughout this paper, the question of how to sample from the constrained
distributions required by NS has not been discussed in detail. For the
examples, I used the original Skilling (2006) suggestion of using a short MCMC
run initialised at a surviving point. In some problems, this will not be
effective, and this will cause the resulting estimates to be incorrect, as in
standard Nested Sampling applications. Testing the sensitivity of the results
to changes in the number of MCMC steps used is a simple way to increase one’s
confidence in the results, but is not a guarantee (no general guarantee can
exist, because any algorithm based on sampling a subset of a space might miss
a crucial but hard-to-locate dominant spike). Despite this weakness, I hope
the algorithm will be useful in a wide range of applications.

###### Acknowledgements.

Acknowledgements I would like to thank John Skilling (Maximum Entropy Data
Consultants) and the anonymous referee for their comments. It is a pleasure to
thank the following people for interesting and helpful conversations about
this topic, and/or comments on a draft: Ruth Angus (Flatiron Institute), Ewan
Cameron (Oxford), James Curran (Auckland), Tom Elliott (Auckland), David Hogg
(NYU), Kevin Knuth (SUNY Albany), Thomas Lumley (Auckland), Iain Murray
(Edinburgh), Jared Tobin (jtobin.io). This work was supported by Centre for
eResearch at the University of Auckland. Conflicts of Interest The authors
declare no conflicts of interest.

## References

  * Shannon (1948) Shannon, C.E.  A mathematical theory of communication.  The Bell System Technical Journal 1948, 27, 379–423.
  * Cover and Thomas (2012) Cover, T.M.; Thomas, J.A.  Elements of information theory; John Wiley & Sons, 2012.
  * Knuth (2005) Knuth, K.H.  Toward question-asking machines: the logic of questions and the inquiry calculus 2005.
  * Knuth and Skilling (2012) Knuth, K.H.; Skilling, J.  Foundations of inference.  Axioms 2012, 1, 38–73.
  * Caticha and Giffin (2006) Caticha, A.; Giffin, A.  Updating Probabilities.  Bayesian Inference and Maximum Entropy Methods In Science and Engineering. AIP Publishing, 2006, Vol. 872, pp. 31–42.
  * Szabó (2014) Szabó, Z.  Information Theoretical Estimators Toolbox.  Journal of Machine Learning Research 2014, 15, 283–287.
  * Jaynes (2003) Jaynes, E.T.  Probability theory: The logic of science; Cambridge university press, 2003.
  * MacKay (2003) MacKay, D.J.  Information theory, inference and learning algorithms; Cambridge university press, 2003.
  * Caticha (2008) Caticha, A.  Lectures on probability, entropy, and statistical physics.  arXiv preprint arXiv:0808.0012 2008.
  * Bernardo (2005) Bernardo, J.M.  Reference analysis.  Handbook of statistics 2005, 25, 17–90.
  * Skilling (2006) Skilling, J.  Nested sampling for general Bayesian computation.  Bayesian analysis 2006, 1, 833–859.
  * Feroz et al. (2009) Feroz, F.; Hobson, M.; Bridges, M.  MultiNest: an efficient and robust Bayesian inference tool for cosmology and particle physics.  MNRAS 2009, 398, 1601–1614.
  * Brewer et al. (2011) Brewer, B.J.; Pártay, L.B.; Csányi, G.  Diffusive nested sampling.  Statistics and Computing 2011, 21, 649–656.
  * Handley et al. (2015) Handley, W.; Hobson, M.; Lasenby, A.  POLYCHORD: next-generation nested sampling.  MNRAS 2015, 453, 4384–4398.
  * Knuth et al. (2015) Knuth, K.H.; Habeck, M.; Malakar, N.K.; Mubeen, A.M.; Placek, B.  Bayesian evidence and model selection.  Digital Signal Processing 2015, 47, 50–67.
  * Pullen and Morris (2014) Pullen, N.; Morris, R.J.  Bayesian model comparison and parameter inference in systems biology using nested sampling.  PloS one 2014, 9, e88419.
  * Brewer and Donovan (2015) Brewer, B.J.; Donovan, C.P.  Fast Bayesian inference for exoplanet discovery in radial velocity data.  MNRAS 2015, 448, 3206–3214.
  * Pártay et al. (2010) Pártay, L.B.; Bartók, A.P.; Csányi, G.  Efficient sampling of atomic configurational spaces.  The Journal of Physical Chemistry B 2010, 114, 10502–10512.
  * Baldock et al. (2016) Baldock, R.J.; Pártay, L.B.; Bartók, A.P.; Payne, M.C.; Csányi, G.  Determining pressure-temperature phase diagrams of materials.  Physical Review B 2016, 93, 174108.
  * Martiniani et al. (2014) Martiniani, S.; Stevenson, J.D.; Wales, D.J.; Frenkel, D.  Superposition enhanced nested sampling.  Physical Review X 2014, 4, 031034.
  * Henderson et al. (2017) Henderson, R.W.; Goggans, P.M.; Cao, L.  Combined-chain nested sampling for efficient Bayesian model comparison.  Digital Signal Processing 2017, 70, 84 – 93.
  * Walter (2017) Walter, C.  Point process-based Monte Carlo estimation.  Statistics and Computing 2017, 27, 219–236.
  * Brewer and Foreman-Mackey (2016) Brewer, B.J.; Foreman-Mackey, D.  DNest4: Diffusive Nested Sampling in C++ and Python.  Journal of Statistical Software, accepted. arxiv: 1606.03757 2016, [arXiv:stat.CO/1606.03757].
  * Bretthorst (2001) Bretthorst, G.L.  Nonuniform sampling: Bandwidth and aliasing.  AIP conference proceedings. AIP, 2001, Vol. 567, pp. 1–28.
  * Gregory (2005) Gregory, P.C.  A Bayesian Analysis of Extrasolar Planet Data for HD 73526.  ApJ 2005, 631, 1198–1214.
  * Taleb (2007) Taleb, N.N.  The Black Swan: The Impact of the Highly Improbable.  ISBN: 978-1400063512 2007.

##  Appendix A Precisional questions

The central issue about the value of a parameter $\boldsymbol{\theta}$ asks
“what is the value of $\boldsymbol{\theta}$, precisely?”. However, in practice
we often don’t need or want to know $\boldsymbol{\theta}$ to arbitrary
precision. Using the central issue can lead to counterintuitive results if you
don’t keep its specific definition in mind. For example, suppose $x$ could
take any integer value from 1 to 1 billion. If we learned the last digit of
$x$, we will have ruled out nine tenths of the possibilities, and therefore
obtained a lot of information about the central issue. However, this
information might be useless for practical purposes. If we learned the final
digit was a 9, $x$ could still be 9, or 19, or 29, any number ending in 9 up
to 999,999,999.

In practice, we may want to ask a question that is different from the central
issue. For example, suppose $x\in\\{1,...,10\\}$ and we want to know the value
of $x$ to within a tolerance of $\pm 1$. Any of the following statements would
resolve the issue:

  * •

$x\in\\{1,2,3\\}$ and anything that implies it

  * •

$x\in\\{2,3,4\\}$ and anything that implies it

  * •

$x\in\\{3,4,5\\}$ and anything that implies it

  * •

and so on.

According to Knuth (2005), the entropy of a question is computed by applying
the sum rule over all statements that would answer the question.

We can write the question of interest, $Q$, as a union of ideal questions (the
ideal questions are downsets of statements: the set of all statements that
imply $S$ is the downset of $S$ and is denoted $\downarrow S$):

$\displaystyle Q$
$\displaystyle=\left[\downarrow(x\in\\{1,2,3\\})\right]\cup\left[\downarrow(x\in\\{2,3,4\\})\right]\cup\left[\downarrow(x\in\\{3,4,5\\})\right]\cup...$
(67)

The entropy of this question is

$\displaystyle H(Q)$ $\displaystyle=h_{1,2,3}$ (68)
$\displaystyle\quad\quad+\left(h_{2,3,4}-h_{2,3}\right)$ (69)
$\displaystyle\quad\quad+\left(h_{3,4,5}-h_{3,4}\right)$ (70)
$\displaystyle\quad\quad+...$ (71)
$\displaystyle\quad\quad+\left(h_{8,9,10}-h_{8,9}\right).$ (72)

where $h_{x}=-P(x)\log P(x)$.

### Continuous case

Consider a probability density function $f(x)$, defined on the real line. The
probability contained in an interval $[x_{0},x_{0}+r]$, which has length $r$,
is

$\displaystyle\int_{x_{0}}^{x_{0}+r}f(x)\,dx$
$\displaystyle=F(x_{0}+r)-F(x_{0})$ (73)

where $F(x)$ is the cumulative distribution function (CDF). This probability
can be considered as a function of two variables, $x_{0}$ and $r$, which I
will denote by $P(,)$:

$\displaystyle P(x_{0},r)$ $\displaystyle=F(x_{0}+r)-F(x_{0}).$ (74)

The contribution of such an interval to an entropy expression is $-P\log P$,
i.e.,

$\displaystyle Q(x_{0},r)$ $\displaystyle=-P(x_{0},r)\log P(x_{0},r).$ (75)

Consider the rate of change of $Q$ as the interval is shifted to the right,
but with its width $r$ held constant:

$\displaystyle\frac{\partial Q}{\partial x_{0}}$
$\displaystyle=\frac{\partial}{\partial x_{0}}\left[-P\log P\right]$ (76)
$\displaystyle=-\left[1+\log P(x_{0},r)\right]\frac{\partial P}{\partial
x_{0}}$ (77) $\displaystyle=-\left[1+\log
P(x_{0},r)\right]\left[f(x_{0}+r)-f(x_{0})\right].$ (78)

Consider also the rate of change of $Q$ as the interval is expanded to the
right, while keeping its left edge fixed:

$\displaystyle\frac{\partial Q}{\partial r}$
$\displaystyle=\frac{\partial}{\partial r}\left[-P(x_{0},r)\log
P(x_{0},r)\right]$ (79) $\displaystyle=-\left[1+\log
P(x_{0},r)\right]\frac{\partial P}{\partial r}$ (80)
$\displaystyle=-\left[1+\log P(x_{0},r)\right]f(x_{0}+r).$ (81)

The entropy of the precisional question is built up from $Q$ terms. The extra
entropy from adding the interval $[x_{0},x_{0}+r]$ and removing the overlap
$[x_{0},x_{0}+r-h]$, for small $h$, is

$\displaystyle\delta H$ $\displaystyle=Q(x_{0},r)-Q(x_{0},r-h)$ (82)
$\displaystyle=h\frac{\partial Q}{\partial r}$ (83)
$\displaystyle=-h\left[1+\log P(x_{0},r)\right]f(x_{0}+r).$ (84)

Therefore, the overall entropy is

$\displaystyle H$ $\displaystyle=-\int_{-\infty}^{\infty}\left[1+\log
P(x-r,r)\right]f(x)\,dx$ (85)
$\displaystyle=-\int_{-\infty}^{\infty}\left(1+\log\left[F(x)-F(x-r)\right]\right)f(x)\,dx.$
(86)

Interestingly, calculating the log-probability of an interval to the right of
$x$ gives an equivalent result:

$\displaystyle H^{\prime}$
$\displaystyle=-\int_{-\infty}^{\infty}\left(1+\log\left[F(x+r)-F(x)\right]\right)f(x)\,dx$
(87) $\displaystyle=H.$ (88)

This is not quite the same as the expected log-probability calculated by the
version of the algorithm proposed in this paper (when the tolerance is not
small). However, the algorithm can be made to estimate the entropy of the
precisional question, by redefining the distance function. For a one-
dimensional quantity of interest, the distance function can be defined such
that $d(x;x_{\rm ref})$ can only ever be below $r$ if $x<x_{\rm ref}$ (or
alternatively, if $x>x_{\rm ref}$).

##  Appendix B Software

A C++ implementation of the algorithm is available in a git repository located
at

    https://github.com/eggplantbren/InfoNest

and can be obtained using the following git command, executed in a terminal:

    git clone https://github.com/eggplantbren/InfoNest

The following will compile the code and execute the first example from the
paper:

    cd InfoNest/cpp
    make
    ./main

The algorithm will run for 1000 ‘reps’, i.e., 1000 samples of
$\boldsymbol{\theta}_{\rm ref}$, which is time consuming. Output is saved to
output.txt. At any time, you can execute the Python script postprocess.py to
get an estimate of the depth:

    python postprocess.py

By default, postprocess.py estimates the depth with a tolerance of
$r=10^{-3}$. This value can be changed by calling the postprocess function
with a different value of its argument tol. E.g., postprocess.py can be edited
so that its last line is postprocess(tol=0.01) instead of postprocess().

The numerical parameters, and the choice of which problem is being solved, are
specified in main.cpp. The default problem is the first example from the
paper. For Example 2, since it is a conditional entropy which is being
estimated (which requires a slight modification to the algorithm), an
additional argument InfoNest::Mode::conditional_entropy must be passed to the
InfoNest::execute function.


t(\vec{n},\vec{m}):=(t_{1}(\vec{n}),\dots,t_{\mathtt{A}}(\vec{n}),t_{\mathtt{A}+1}(\vec{m}),\dots,t_{\mathtt{A}+\mathtt{B}}(\vec{m}))\,,$
(B.32)

defined by (B.26), it results

$c(\mathtt{A}){\Big{(}\sum_{a=1}^{\mathtt{A}}n_{a}+\sum_{b=1}^{\mathtt{B}}m_{b}\Big{)}^{-\tau_{2}}}\leq|{\rho(x(\vec{n},\vec{m}),t(\vec{n},\vec{m}))}|\leq
C(\mathtt{A})\Big{(}\sum_{a=1}^{\mathtt{A}}n_{a}+\sum_{b=1}^{\mathtt{B}}m_{b}\Big{)}^{-1}$
(B.33)

with $\tau_{2}:=\mathtt{A}+1+12\binom{\mathtt{A}}{2}$.

###### Proof.

The upper bound (B.33) directly follows by (B.27). The lower bound (B.33)
follows by (B.27) and the fact that, since $1\leq n_{1}<\ldots<n_{\mathtt{A}}$
are all distinct,

$\displaystyle\Big{|}(gx_{a}^{2}x_{0}^{4}+\frac{\gamma^{2}}{4}t_{a}^{2}x_{0}^{6})x_{b}^{6}-(gx_{b}^{2}x_{0}^{4}+\frac{\gamma^{2}}{4}t_{b}^{2}x_{0}^{6}))x_{a}^{6}\Big{|}$
$\displaystyle=x_{0}^{12}(\vec{n},\vec{m})(gn_{a}+\frac{\gamma^{2}}{4}\tanh(\mathtt{h}n_{a}))(gn_{b}+\frac{\gamma^{2}}{4}\tanh(\mathtt{h}n_{b}))\Big{|}\frac{n_{a}^{3}}{gn_{a}+\frac{\gamma^{2}}{4}\tanh(\mathtt{h}n_{a})}-\frac{n_{b}^{3}}{gn_{b}+\frac{\gamma^{2}}{4}\tanh(\mathtt{h}n_{b})}\Big{|}$
$\displaystyle\geq x_{0}(\vec{n},\vec{m})^{12}\,g^{2}\,\min_{y\geq
1}\Big{(}\frac{{\rm d}}{{\rm
d}y}\frac{y^{3}}{gy+\frac{\gamma^{2}}{4}\tanh(\mathtt{h}y)}\Big{)}|n_{a}-n_{b}|\geq
c\,x_{0}(\vec{n},\vec{m})^{12}$

for some constant $c>0$, having used in the last passage that, for any $y\geq
1$, $g>0$, $\gamma\in\mathbb{R}$ and $\mathtt{h}>0$,

$\displaystyle\frac{{\rm d}}{{\rm
d}y}\frac{y^{3}}{gy+\frac{\gamma^{2}}{4}\tanh(\mathtt{h}y)}$
$\displaystyle=\frac{4y^{2}(8gy+3\gamma^{2}\tanh(\mathtt{h}y)-\gamma^{2}\mathtt{h}y\,{\rm
sech}^{2}(\mathtt{h}y))}{(4gy+\gamma^{2}\tanh(\mathtt{h}y))^{2}}$
$\displaystyle\geq$
$\displaystyle\frac{2y^{2}}{16g^{2}y^{2}+\gamma^{4}}\big{(}8g+\gamma^{2}(3\tanh(\mathtt{h}y)-\mathtt{h}y\,{\rm
sech}^{2}(\mathtt{h}y))\big{)}\geq\frac{2}{16g^{2}+\gamma^{4}}8g$

using that the function $3\tanh(y)-y\,{\rm sech}^{2}(y)$ is positive for
$y>0$. ∎

We show now that the assumptions (H1) and (H2) of Theorem B.3 hold true.

Verification of (H1). By (B.31), if $\rho(x,t)\neq 0$ then, by (B.31),

$x_{a}\neq 0\,,\ \forall a=0,\ldots,\mathtt{A}\,,\quad t_{a}\not=0\,,\ \forall
a=1,\ldots,\mathtt{A}\,.$ (B.34)

In particular on the set $\\{(x,t)\in X\colon\rho(x,t)\neq 0\\}\times Y$ the
function $f_{\vec{c},\vec{d}}$ in (B.30) is real analytic.
Verification of (H2). For any $(x,t)$ such that $\rho(x,t)\neq 0$, the
analytic function $\kappa\mapsto f_{\vec{c},\vec{d}}\,(x,t,\kappa)$ possesses
only a finite number of zeros as a consequence of the next lemma.

######  Lemma B.8.

For any $(x,t)$ such that $\rho(x,t)\neq 0$, the analytic function
$\kappa\mapsto f_{\vec{c},\vec{d}}\,(x,t,\kappa)$ is not identically zero in
$\mathcal{I}$.

###### Proof.

Assume by contradiction that there exists

$(x,t)\in X\,,\ x=(x_{a})_{0\leq a\leq\mathtt{A}}\,,\ t=(t_{a})_{0\leq
a\leq\mathtt{A}+\mathtt{B}}\quad\text{with}\quad\rho(x,t)\neq 0\quad\text{such
\ that}\quad f_{\vec{c},\vec{d}}\,(x,t,\kappa)=0$

for any $\kappa$ in the interval $\mathcal{I}$. Then, by analyticity, the
function $\kappa\mapsto f_{\vec{c},\vec{d}}\,(x,t,\kappa)$ is identically zero
also on the larger interval $(-\delta,+\infty)$ where $\delta:=\min_{1\leq
a\leq\mathtt{A}}t_{a}^{2}\frac{\gamma^{2}}{4}x_{0}^{6}$. Note that $\delta>0$
as $t_{a}^{2}$, $x_{0}^{6}>0$ by (B.34). In particular, for any
$l\in\mathbb{N}$, all the derivatives
$\partial_{\kappa}^{l}f_{\vec{c},\vec{d}}\,(x,t,\kappa)\equiv 0$ are zero in
the interval $\kappa\in(-\delta,+\infty)$.

We now compute such derivatives at $\kappa=0$ differentiating (B.30). The
derivatives of the function $\lambda(y,s,x_{0},\kappa)$ defined in (B.29) are
given by, for suitable constants $C_{l}\neq 0$,

$\partial_{\kappa}^{l}\lambda(y,s,x_{0},\kappa)=C_{l}y^{6l}\,\lambda(y,s,x_{0},\kappa)^{1-2l}\,,\quad\forall
l\in\mathbb{N}\,.$

Thus we obtain

$\partial_{\kappa}^{l}\lambda(y,s,x_{0},\kappa)|_{\kappa=0}=C_{l}\,\mu(y,s,x_{0})^{l}\,\lambda(y,s,x_{0},0)\quad\text{where}\quad\mu(y,s,x_{0}):=\frac{y^{6}}{gy^{2}x_{0}^{4}+\frac{\gamma^{2}}{4}s^{2}x_{0}^{6}}\,,$
(B.35)

and, recalling (B.30),

$\partial_{\kappa}^{l}f_{\vec{c},\vec{d}}\,(x,t,\kappa)|_{\kappa=0}=C_{l}\sum_{a=1}^{\mathtt{A}}c_{a}t_{a}\mu(x_{a},t_{a},x_{0})^{l}\,\lambda(x_{a},t_{a},x_{0},0)\,,\quad\forall
l\in\mathbb{N}\,.$

As a consequence, the conditions
$\partial_{\kappa}^{l}f_{\vec{c},\vec{d}}\,(x,t,\kappa)|_{\kappa=0}=0$ for any
$l=1,\ldots,\mathtt{A}$ imply that

$A(x,t)\vec{r}=0$ (B.36)

where $A(x,t)$ is the $\mathtt{A}\times\mathtt{A}$ matrix

$A(x,t):=\begin{pmatrix}\mu(x_{1},t_{1},x_{0})\lambda(x_{1},t_{1},x_{0},0)&\cdots&\mu(x_{\mathtt{A}},t_{\mathtt{A}},x_{0})\lambda(x_{\mathtt{A}},t_{\mathtt{A}},x_{0},0)\\\
\mu(x_{1},t_{1},x_{0})^{2}\lambda(x_{1},t_{1},x_{0},0)&\cdots&\mu(x_{\mathtt{A}},t_{\mathtt{A}},x_{0})^{2}\lambda(x_{\mathtt{A}},t_{\mathtt{A}},x_{0},0)\\\
\vdots&\ddots&\vdots\\\
\mu(x_{1},t_{1},x_{0})^{\mathtt{A}}\lambda(x_{1},t_{1},x_{0},0)&\ldots&\mu(x_{\mathtt{A}},t_{\mathtt{A}},x_{0})^{\mathtt{A}}\lambda(x_{\mathtt{A}},t_{\mathtt{A}},x_{0},0)\end{pmatrix}$

and $\vec{r}$ is the vector

$\vec{r}:=\begin{bmatrix}r_{1}\\\ \vdots\\\
r_{\mathtt{A}}\end{bmatrix}\,,\quad
r_{a}:=c_{a}t_{a}\in\mathbb{Z}\setminus\\{0\\}\,,\
\forall\,a=1,\dots,\mathtt{A}\,,$

because by assumption each $c_{a}\neq 0$ and (B.34) holds. Since $\vec{r}\neq
0$, we deduce by (B.36) that the matrix $A(x,t)$ has zero determinant. On the
other hand, by the multilinearity of the determinant,

$\displaystyle\det A(x,t)$
$\displaystyle=\prod_{a=1}^{\mathtt{A}}\mu(x_{a},t_{a},x_{0})\lambda(x_{a},t_{a},x_{0},0)\det\begin{pmatrix}1&\cdots&1\\\
\mu(x_{1},t_{1},x_{0})&\cdots&\mu(x_{\mathtt{A}},t_{\mathtt{A}},x_{0})\\\
\vdots&\ddots&\vdots\\\
\mu(x_{1},t_{1},x_{0})^{\mathtt{A}-1}&\ldots&\mu(x_{\mathtt{A}},t_{\mathtt{A}},x_{0})^{\mathtt{A}-1}\end{pmatrix}$
$\displaystyle=\prod_{a=1}^{\mathtt{A}}\mu(x_{a},t_{a},x_{0})\lambda(x_{a},t_{a},x_{0},0)\,\prod_{1\leq
a<b\leq\mathtt{A}}(\mu(x_{a},t_{a},x_{0})-\mu(x_{b},t_{b},x_{0}))\,.$ (B.37)

The condition $\rho(x,t)\neq 0$ implies, by (B.34), (B.35) and (B.29), that

$\prod_{a=1}^{\mathtt{A}}\mu(x_{a},t_{a},x_{0})\lambda(x_{a},t_{a},x_{0},0)\neq
0\,,$

and, in view of (B.37), the determinant $\det A(x,t)=0$ if only if
$\mu(x_{a},t_{a},x_{0})=\mu(x_{b},t_{b},x_{0})$ for some $1\leq
a<b\leq\mathtt{A}$. By the definition of the function
$(y,s,x_{0})\to\mu(y,s,x_{0})$ in (B.35) it follows that

$\mu(x_{a},t_{a},x_{0})=\mu(x_{b},t_{b},x_{0})\quad\Rightarrow\quad\big{(}gx_{a}^{2}x_{0}^{4}+\frac{\gamma^{2}}{4}t_{a}^{2}x_{0}^{6}\big{)}x_{b}^{6}-\big{(}gx_{b}^{2}x_{0}^{4}+\frac{\gamma^{2}}{4}t_{b}^{2}x_{0}^{6}\big{)}x_{a}^{6}=0\,.$

In view of (B.31) this contradicts $\rho(x,t)\neq 0$. ∎

We have verified assumptions (H1) and (H2) of Theorem B.3. We thus conclude
that there are $N_{0}\in\mathbb{N}$, $\alpha_{0},\delta,C>0$, such that for
any $\alpha\in(0,\alpha_{0}]$, any $N\in\mathbb{N}$, $N\geq N_{0}$, any
$(x,t)\in X$ with $\rho(x,t)\neq 0$,

${\rm meas}\big{\\{}\kappa\in\mathcal{I}\colon\ \ \
|f_{\vec{c},\vec{d}}\,(x,t,\kappa)|\leq\alpha|\rho(x,t)|^{N}\big{\\}}\leq
C\alpha^{\delta}\,|\rho(x,t)|^{N\delta}\,.$ (B.38)

For $\vec{n}=(n_{1},\ldots,n_{\mathtt{A}})\in\mathbb{N}^{\mathtt{A}}$ with
$1\leq n_{1}<\ldots<n_{\mathtt{A}}$ and
$\vec{m}=(m_{1},\ldots,m_{\mathtt{B}})\in\mathbb{N}^{\mathtt{B}}$ we consider
the set

$\mathcal{B}_{\vec{c},\vec{d},\vec{m},\vec{n}}(\alpha,N):=\left\\{\kappa\in\mathcal{I}\colon\
\ \
|f_{\vec{c},\vec{d}}\,(x(\vec{n},\vec{m}),t(\vec{n},\vec{m}),\kappa)|\leq\alpha|\rho(x(\vec{n},\vec{m}),t(\vec{n},\vec{m}))|^{N}\right\\}$
(B.39)

where $x(\vec{n},\vec{m})$ and $t(\vec{n},\vec{m})$ are defined in (B.32). By
(B.33) we deduce that $\rho(x(\vec{n},\vec{m}),t(\vec{n},\vec{m}))\neq 0$, and
(B.38) implies that

${\rm meas}\,\mathcal{B}_{\vec{c},\vec{d},\vec{m},\vec{n}}(\alpha,N)\leq
C\alpha^{\delta}|\rho(x(\vec{n},\vec{m}),t(\vec{n},\vec{m}))|^{N\delta}\,.$
(B.40)

Consider the set

$\mathcal{K}(\alpha,N):=\bigcup_{\begin{subarray}{c}\vec{n}=(n_{1},\ldots,n_{\mathtt{A}})\in\mathbb{N}^{\mathtt{A}},\,1\leq
n_{1}<\ldots<n_{\mathtt{A}}\\\ \vec{m}=(m_{1},\dots
m_{\mathtt{B}})\in\mathbb{N}^{\mathtt{B}}\\\
\vec{c}\in(\mathbb{Z}\setminus\\{0\\})^{\mathtt{A}},\,|\vec{c}|_{\infty}\leq\mathtt{M}\\\
\vec{d}\in(\mathbb{Z}\setminus\\{0\\})^{\mathtt{B}},|\vec{d}|_{\infty}\leq\mathtt{M}\end{subarray}}\mathcal{B}_{\vec{c},\vec{d},\vec{m},\vec{n}}(\alpha,N)\subset\mathcal{I}\,.$
(B.41)

By (B.40) and (B.33) one gets

${\rm meas}\,\mathcal{K}(\alpha,N)\leq
C(\mathtt{A},\mathtt{M})\alpha^{\delta}\sum_{\begin{subarray}{c}n_{1},\ldots,n_{\mathtt{A}}\in\mathbb{N}\\\
m_{1},\dots,m_{\mathtt{B}}\in\mathbb{N}\end{subarray}}\frac{1}{(\sum_{a=1}^{\mathtt{A}}n_{a}+\sum_{b=1}^{\mathtt{B}}m_{b})^{N\delta}}\leq
C^{\prime}(\mathtt{A},\mathtt{M})\alpha^{\delta}$ (B.42)

for some finite constant $C^{\prime}(\mathtt{A},\mathtt{M})<+\infty$, provided
$N\delta>\mathtt{A}+\mathtt{B}$. We fix

$\underline{N}:=[(\mathtt{A}+\mathtt{B})\delta^{-1}]+1\quad\text{and
define}\quad{\mathcal{K}}_{\alpha}:=\mathcal{K}(\alpha,\underline{N})$

whose measure satisfies $|{\mathcal{K}}_{\alpha}|=O(\alpha^{\delta})$, in view
of (B.42).

In conclusion, for any $\kappa\in\mathcal{I}\setminus\mathcal{K}_{\alpha}$,
for any $\vec{n}=(n_{1},\ldots,n_{\mathtt{A}})$ with $1\leq
n_{1}<\ldots<n_{\mathtt{A}}$ and $\vec{m}\in\mathbb{N}^{\mathtt{B}}$, any
$\vec{c}\in(\mathbb{Z}\setminus\\{0\\})^{\mathtt{A}}$ with
$|\vec{c}\,|_{\infty}\leq\mathtt{M}$ and
$\vec{d}\,\in(\mathbb{Z}\setminus\\{0\\})^{\mathtt{B}}$ with
$|\vec{d}\,|_{\infty}\leq\mathtt{M}$, it results, by (B.41) and (B.39), that

$|{f_{\vec{c},\vec{d}}\,(x(\vec{n},\vec{m}),t(\vec{n},\vec{m}),\kappa)}|>\alpha|\rho(x(\vec{n},\vec{m}),t(\vec{n},\vec{m}))|^{\underline{N}}\stackrel{{\scriptstyle\eqref{rho.estgh}}}{{\geq}}\frac{c(\mathtt{A})\alpha}{\left(\sum_{a=1}^{\mathtt{A}}n_{a}+\sum_{b=1}^{\mathtt{B}}m_{b}\right)^{\tau_{2}\underline{N}}}\,.$
(B.43)

Recalling the definition of $f_{\vec{c},\vec{d}}$ in (B.30) and
$x_{0}(\vec{n},\vec{m})$ in (B.26), the lower bound (B.43) implies (B.24) with
$\tau:=\tau_{2}\underline{N}-3$ (cfr. (B.28)) and re-denoting $\alpha
c(\mathtt{A})\leadsto\alpha$. ∎

### References

  * [1] Alazard T., Baldi P., _Gravity capillary standing water waves_. Arch. Rat. Mech. Anal. 217(3), 741-830, 2015.
  * [2] Alazard T., Burq N., Zuily C.,  On the water-wave equations with surface tension.  Duke Math. J., 158, 413-499, 2011.
  * [3] Alazard T., Burq N., Zuily C., Cauchy theory for the gravity water waves system with non-localized initial data. Ann. Inst. H. Poincaré Anal. Non Linéaire 33.2, pp. 337-395, 2016.
  * [4] Alazard T., Burq N., Zuily C.,  On the Cauchy problem for gravity water waves.  Invent. Math., 198, 71–163, 2014.
  * [5] Alazard T., Delort J-M.,  Global solutions and asymptotic behavior for two dimensional gravity water waves.  Ann. Sci. Éc. Norm. Supér., 48(5): 1149-1238, 2015.
  * [6] Alazard T., Métivier G., _Paralinearization of the Dirichlet to Neumann operator, and regularity of the three dimensional water waves_. Comm. Partial Differential Equations 34(10-12), 1632-1704, 2009.
  * [7] Ambrose D., Masmoudi N., The zero surface tension limit two-dimensional water waves. Comm. on Pure and Applied Math., 58(10): 1287-1315, 2005.
  * [8] Baldi P., Berti M., Haus E., Montalto R., _Time quasi-periodic gravity water waves in finite depth_. Invent. Math. 214(2): 739-911, 2018.
  * [9] Bambusi D., Birkhoff normal form for some nonlinear PDEs. Comm. Math. Phys. 234(2): 253-285, 2003.
  * [10] Bambusi D., Asymptotic Stability of Ground States in Some Hamiltonian PDEs with Symmetry. Comm. Math. Phys. 320(2): 499-542, 2013.
  * [11] Bambusi D., Delort J.-M., Grébert B., Szeftel J., Almost global existence for Hamiltonian semilinear Klein-Gordon equations with small Cauchy data on Zoll manifolds. Comm. Pure Appl. Math. 60, 2007.
  * [12] Bambusi D., Feola R., Montalto R., Almost global existence for some Hamiltonian PDEs with small Cauchy data on general tori. Arxiv.org/abs/2208.00413
  * [13] Bambusi D., Grébert B., Birkhoff normal form for Pdes with tame modulus. Duke Math. J., 135(3): 507-567, 2006.
  * [14] Bambusi D., Maspero A., Birkhoff coordinates for the Toda Lattice in the limit of infinitely many particles with an application to FPU. J. Funct. Anal., 270(5): 1818-1887, 2016.
  * [15] Bambusi D., Maspero A., Freezing of Energy of a Soliton in an External Potential. Comm. Math. Phys., 344(1): 155-191, 2016.
  * [16] Berti M., Delort J.-M.,  Almost Global Solutions of Capillary-gravity Water Waves Equations on the Circle.  UMI Lecture Notes 2018, ISBN 978-3-319-99486-4.
  * [17] Berti M., Feola R., Franzoi L., Quadratic life span of periodic gravity-capillary water waves. Water Waves 3(1): 85-115, 2021.
  * [18] Berti M., Feola R., Pusateri F., Birkhoff Normal Form and Long Time Existence for Periodic Gravity Water Waves. Comm. Pure Applied Math., doi.org/10.1002/cpa.22041, 2022.
  * [19] Berti M., Franzoi L., Maspero A., Traveling quasi-periodic water waves with constant vorticity. Arch. Ration. Mech. Anal., 240, 99–202, 2021.
  * [20] Berti M., Franzoi L., Maspero A., Pure gravity traveling quasi-periodic water waves with constant vorticity. Arxiv. 2101.12006, to appear in Comm. Pure Applied Math.
  * [21] Berti M., Montalto R., _Quasi-periodic standing wave solutions of gravity-capillary water waves_. Memoires American Mathematical Society. Volume 263, 1273, ISSN 0065-9266, 2020.
  * [22] Berti M., Maspero A., Murgante F.,  Local well posedness of the Euler-Korteweg equations on $\mathbb{T}^{d}$.  J. Dyn. Diff. Equat., 33: 1475-1513, 2021.
  * [23] Beyer K., Günther M., On the Cauchy problem for a capillary drop. I. Irrotational motion. Math. Methods Appl. Sci. 21(12): 1149-1183, 1998.
  * [24] Bernier J., Faou E. and Grébert B., Rational normal forms and stability of small solutions to nonlinear Schrödinger equations. Annals of PDE 6, 14, 2020.
  * [25] Bernier J., Grébert B. Birkhoff normal forms for Hamiltonian PDEs in their energy space. Journal de l’Ecole polytechnique Mathematiques, Tome 9, pp. 681-745, 2022.
  * [26] Biasco L., Massetti J., Procesi M., An abstract Birkhoff normal form theorem and exponential type stability of the 1d NLS. Comm. Math. Phys. 375, no. 3, 2089–2153, 2020.
  * [27] Castro A., Córdoba D., Fefferman C., Gancedo F., Gómez-Serrano J., Finite time singularities for the free boundary incompressible Euler equations. Annals of Math., 178, 1061-1134, 2013.
  * [28] Constantin A., Ivanov R.I., Prodanov E.M., Nearly-Hamiltonian structure for water waves with constant vorticity. J. Math. Fluid Mech. 10: 224-237, 2008.
  * [29] Coutand D., Shkoller S., Well-posedness of the free-surface incompressible Euler equations with or without surface tension. J. Amer. Math. Soc. 20(3): 829-930, 2007.
  * [30] Coutand D., Shkoller S., On the finite-time splash and splat singularities for the 3-D free-surface Euler equations. Comm. Math. Phys., 325, 143-183, 2014.
  * [31] Craik A. D., The origins of water wave theory. Annual Review of Fluid Mechanics, Vol. 36, pp. 128. Annual Reviews, Palo Alto, CA, 2004.
  * [32] Craig W., An existence theory for water waves and the Boussinesq and Korteweg-de Vries scaling limits. Comm. Partial Differential Equations, 10(8): 787-1003, 1985.
  * [33] Craig W., Sulem C.,  Numerical simulation of gravity waves.  J. Comput. Phys., 108(1): 73-83, 1993.
  * [34] Cuccagna S.,  The Hamiltonian Structure of the Nonlinear Schrodinger Equation and the Asymptotic Stability of its Ground States.  Comm. Math. Phys., 305(2): 279-331, 2011.
  * [35] Cuccagna S.,  On asymptotic stability of moving ground states of the nonlinear Schrödinger equation.  Trans. Amer. Math. Society, 366(6): 2827-2888, 2014.
  * [36] Delort J.-M., A quasi-linear Birkhoff normal forms method. Application to the quasi-linear Klein-Gordon equation on S1. Astérisque No. 341, vi+113 pp., 2012.
  * [37] Delort J.-M., Quasi-linear perturbations of Hamiltonian Klein-Gordon equations on spheres. Mem. Amer. Math. Soc. 234, no. 1103, vi+80 pp, 2015.
  * [38] Delort J.-M., Szeftel J.,  Long-time existence for small data nonlinear Klein-Gordon equations on tori and spheres.  Int. Math. Res. Not., 37: 1897-1966, 2004.
  * [39] Deng Y., Ionescu A. D., Pausader B., Pusateri F., Global solutions for the 3D gravity-capillary water waves system. Acta Math. 219(2): 213-402, 2017.
  * [40] Deng Y., Ionescu A. D., Pusateri F., On the wave turbulence theory of 2D gravity waves, I: deterministic energy estimates, preprint, https://arxiv.org/abs/2211.10826.
  * [41] Feola R., Grébert B., and Iandoli F., Long time solutions for quasilinear Hamiltonian perturbations of Schrödinger and Klein-Gordon equations on tori, to appear on Analysis and PDEs.
  * [42] Feola R., Giuliani F., Quasi-periodic traveling waves on an infinitely deep fluid under gravity. ArXiv:2005.08280, to appear on Memoires American Mathematical Society.
  * [43] Feola R., Iandoli F., Local well-posedness for quasi-linear NLS with large Cauchy data on the circle. Ann. Inst. H. Poincaré Anal. Non Linéaire, 36(1): 119-164, 2019.
  * [44] Feola R., Iandoli F., Long time existence for fully nonlinear NLS with small Cauchy data on the circle. Annali della Scuola Normale Superiore di Pisa, XXII(1): 109-182, 2021.
  * [45] Feola R., Iandoli F., A non-linear Egorov theorem and Poincaré-Birkhoff normal forms for quasi-linear pdes on the circle.  ArXiv:2002.12448, 2020.
  * [46] Feola R., Iandoli F., Murgante F., Long-time stability of the quantum hydrodynamic system on irrational tori. Math. in Engineering, 4(3), 2022.
  * [47] Germain P., Masmoudi N., Shatah J., Global solutions for the gravity water waves equation in dimension 3. Ann. of Math. (2), 175, 691–754, 2012.
  * [48] Germain P., Masmoudi N., Shatah J., Global solutions for capillary waves equation in dimension 3. Comm. Pure Appl. Math., 68(4): 625-687, 2015.
  * [49] Harrop-Griffiths B., Ifrim M., Tataru D., Finite depth gravity water waves in holomorphic coordinates. Ann. PDE 3, 4, 2017.
  * [50] Hunter J., Ifrim M., Tataru D., Two dimensional water waves in holomorphic coordinates. Comm. Math. Phys. 346, 483-552, 2016.
  * [51] Ifrim M., Tataru D., Two dimensional water waves in holomorphic coordinates II: global solutions. Bull. Soc. Math. France 144, 369-394, 2016.
  * [52] Ifrim M., Tataru D., The lifespan of small data solutions in two dimensional capillary water waves. Arch. Ration. Mech. Anal. 225(3): 1279-1346, 2017.
  * [53] Ifrim M., Tataru D., Two-dimensional gravity water waves with constant vorticity I: Cubic lifespan. Anal. PDE 12(4): 903-967, 2019.
  * [54] Ionescu A., Pusateri F., Global solutions for the gravity water waves system in 2d. Invent. Math., 199(3): 653-804, 2015.
  * [55] Ionescu A., Pusateri F., Global regularity for 2d water waves with surface tension. Mem. Amer. Math. Soc. 256, no. 1227, v+123 pages, 2018.
  * [56] Ionescu A., Pusateri F., Long-time existence for multi-dimensional periodic water waves. Geom. Funct. Anal. 29(3): 811-870, 2019.
  * [57] Iooss G., Plotnikov P., Toland J.,  Standing waves on an infinitely deep perfect fluid under gravity.  Arch. Ration. Mech. Anal., 177(3): 367-478, 2005.
  * [58] John F., Blow-up for quasilinear wave equations in three space dimensions. Comm. Pure Appl. Math., 34, 29-51, 1981.
  * [59] Kuksin S., Perelman G.,  Vey theorem in infinite dimensions and its application to KdV.  Discrete Cont. Dyn. Systems, 27(1): 1-24, 2010.
  * [60] Lannes D., Well-posedness of the water-waves equations. J. Amer. Math. Soc., 18(3): 605-654, 2005.
  * [61] Lindblad H., Well-posedness for the motion of an incompressible liquid with free surface boundary. Annals of Math., 162(1): 109-194, 2005.
  * [62] Ming M., Zhang Z., Well-posedness of the water-wave problem with surface tension. J. Math. Pures Appl., 92(5), 429-455, 2009.
  * [63] Nalimov V. I., The Cauchy-Poisson problem. Dinamika Splosn. Sredy Vyp. 18 Dinamika Zidkost. so Svobod. Granicami, 10-210, 254, 1974.
  * [64] Plotnikov P., Toland J.,  Nash-Moser theory for standing water waves.  Arch. Rat. Mech. Anal., 159, 1–83, 2001.
  * [65] Schweizer B., On the three-dimensional Euler equations with a free boundary subject to surface tension. Ann. Inst. H. Poincaré C Anal. Non Linéaire 22(6): 753-781, 2005.
  * [66] Shatah J., Zeng C., Geometry and a priori estimates for free boundary problems of the Euler equation. Comm. Pure Appl. Math., 61, 698-744, 2008.
  * [67] Shatah J., Zeng C., Local well-posedness for fluid interface problems. Arch. Ration. Mech. Anal., 199, 653-705, 2011.
  * [68] Sideris T., Formation of singularities in three-dimensional compressible fluids. Comm. Math. Phys., 101, 475-485, 1985.
  * [69] Wahlén E., A Hamiltonian formulation of water waves with constant vorticity. Letters in Math. Physics, 79, 303-315, 2007.
  * [70] Wahlén E., Steady periodic capillary-gravity waves with vorticity. SIAM J. Math. Anal. 38, 921-943, 2006.
  * [71] Yosihara H., Gravity waves on the free surface of an incompressible perfect fluid of finite depth. Publ. Res. Inst. Math. Sci. 18, 49-96, 1982.
  * [72] Zakharov V.E., Stability of periodic waves of finite amplitude on the surface of a deep fluid. Zhurnal Prikladnoi Mekhaniki i Teckhnicheskoi Fiziki 9(2): 86-94, 1969.
  * [73] Zakharov V.E., Dyachenko A.I., Is free-surface hydrodynamics an integrable system? Physics Letters A 190 144-148, 1994.
  * [74] Wu S., Well-posedness in Sobolev spaces of the full water wave problem in 2-D. Invent. Math., 130(1):39-72, 1997.
  * [75] Wu S., Well-posedness in Sobolev spaces of the full water wave problem in 3-D. J. Amer. Math. Soc., 12(2):445-495, 1999.
  * [76] Wu S., Almost global wellposedness of the 2-D full water wave problem. Invent. Math., 177, 45–135, 2009.
  * [77] Wu S., Global well-posedness of the 3-D full water wave problem. Invent. Math., 184, 125–220, 2011.
  * [78] Wu S., The quartic integrability and long time existence of steep water waves in 2d, arxiv.org/abs/2010.09117.

Supported by PRIN 2020XB3EFL, Hamiltonian and Dispersive PDEs.