Introducing Inner Nested Sampling

TL;DR

This paper presents Inner Nested Sampling, a Monte Carlo algorithm for estimating moments of functions of Dirichlet distributions within Skilling's Nested Sampling framework.

Contribution

It introduces a novel Monte Carlo algorithm specifically designed for Dirichlet distributions, expanding the applications of Nested Sampling.

Findings

Efficient estimation of moments for Dirichlet functions.

Implementation of Inner Nested Sampling within existing frameworks.

Potential improvements in Bayesian inference methods.

Abstract

In this paper we will give a Monte Carlo algorithm by which the moments of a functions of Dirichlet probability distributions can be estimated. This algorithm is called Inner Nested Sampling and is an implementation of Skilling's general Nested Sampling framework.

Tables

Table 1. Table 1. Distribution 1

	$b_1$	$b_2$
$a_1$	${\theta }_1$	${\theta }_3$	${\theta }_1+{\theta }_3$
$a_2$	${\theta }_2$	${\theta }_4$	${\theta }_2+{\theta }_4$
	${\theta }_1+{\theta }_2$	${\theta }_3+{\theta }_4$	${\theta }_1+{\theta }_2+{\theta }_3+{\theta }_4=1$

Table 2. Table 2. Count data

	$b_1$	$b_2$
$a_1$	$r_1$	$r_3$	$r_1+r_3$
$a_2$	$r_2$	$r_4$	$r_2+r_4$
	$r_1+r_2$	$r_3+r_4$	$r_1+r_2+r_3+r_4=n$

Equations

$$D=\left\{r_1,r_2,\dots ,r_n\right\},$$

$$p\left(\mathbf{\theta }\right)=\frac{\left(\sum r_i-1\right)!}{\left(r_1-1\right)!\left(r_1-2\right)!\cdots \left(r_n-1\right)!}{\theta }_1^{r_1-1}{\theta }_2^{r_2-1}\cdots {\theta }_n^{r_n-1},$$

$$\mathbf{\theta }=\left\{{\theta }_1,{\theta }_2,\dots ,{\theta }_n\right\}.$$

$$p\left(\theta ,D\right)=p\left(\theta \right)L\left(\theta ,D\right).$$

$$Z=\int p\left(\theta ,D\right)d\theta$$

$$p\left(\theta |D\right)=\frac{p\left(\theta ,D\right)}{Z}.$$

$$\text{posterior}=\frac{\text{prior}\times \text{likelihood}}{\text{evidence}}.$$

$$\text{posterior}\propto \text{prior}\times \text{likelihood},$$

$$p\left(\theta |D,M,I\right)=\frac{p\left(\theta |I\right)p\left(D|\theta ,M\right)}{\int p\left(\theta |I\right)p\left(D|\theta ,M\right)d\theta }=\frac{p\left(\theta |I\right)p\left(D|\theta ,M\right)}{p\left(D|M,I\right)},$$

$$p\left(D|M,I\right)=\int p\left(\theta ,D|M,I\right)d\theta =\int p\left(\theta |I\right)p\left(D|\theta ,M\right)d\theta$$

$$p\left(M_j|D,I\right)=\frac{p\left(M_j\right)p\left(D|M_j,I\right)}{{\sum }_{j}p\left(M_j\right).p\left(D|M_j,I\right)}$$

$$p\left(M_j|D,I\right)=\frac{p\left(D|M_j,I\right)}{{\sum }_{j}p\left(D|M_j,I\right)}.$$

$$f\left(x,y\right)=\frac{\sqrt{1-0.7^2}}{2\pi }\exp \left[-\frac{1}{2}\left(x^2+1.4xy+y^2\right)\right]$$

$${\int }_{-5}^5{\int }_{-5}^5\frac{\sqrt{1-0.7^2}}{2\pi }\exp \left[-\frac{1}{2}\left(x^2+1.4xy+y^2\right)\right]dxdy=0.9993$$

$$V_{jk}=f\left({\tilde{x}}_j,{\tilde{y}}_k\right)d{x}_{j}d{y}_k$$

$$\text{volume}\approx \sum _{j=1}^{20}\sum _{k=1}^{20}{V}_{jk}=0.9994$$

$$d{w}_i=d{x}_{j}d{y}_k,f_i=f\left({\tilde{x}}_j,{\tilde{y}}_k\right)$$

$$i\equiv \left(j-1\right)20+k$$

$$V_{jk}=f\left({\tilde{x}}_j,{\tilde{y}}_k\right)d{x}_{j}d{y}_k=f_{i}d{w}_i=A_i$$

$$\text{area}=\sum _{i=1}^{400}{A}_i=\sum _{j=1}^{20}\sum _{k=1}^{20}{V}_{jk}=\text{volume}$$

$$p\left(w\right)=\frac{1}{W},0\le w\le W.$$

$$E\left(w\right)=\frac{W}{2}$$

$$\text{std}\left(w\right)=\frac{W}{2\sqrt{3}}.$$

$$g\left(w_n\right)=f\left({\mathbf{x}}_n\right),n=1,\dots ,N$$

$$p\left(w_{\max }\right)=N{\left(\frac{w_{\max }}{W}\right)}^{N-1}\frac{1}{W}$$

$$E\left(w_{\max }\right)=\left(1-\frac{1}{N+1}\right)W$$

$$\text{std}\left(w_{\max }\right)=\sqrt{\frac{N}{{\left(N+1\right)}^2\left(N+2\right)}}W$$

$$t=\frac{w_{\max }}{W}.$$

$$p\left(t\right)=N{t}^{N-1}.$$

$$f\left({\mathbf{x}}_n^{\left(1\right)}\right),n=1,\dots ,N$$

Full text

1. Introduction

This paper is for the benefit of those readers who are interested in finding moments of functions of Dirichlet probability distributions. In a data analysis we will, typically, have some set of count data, say,

[TABLE]

which by way of a Dirichlet distribution,

[TABLE]

points to some set of underlying probabilities $\boldsymbol{\theta}$, where

[TABLE]

Stated differently, the $\boldsymbol{\theta}$ are not directly observable, they can only be inferred. So, if we wish to assign a function $u$ to $\boldsymbol{\theta}$, then we have to take our uncertainty, in regards to the actual value of the vector $\boldsymbol{\theta}$, into account. But if we do so, then this will give us highly dimensional and highly intractable integrals. In order to evaluate these integrals, we will invoke the Nested Sampling framework.

In order to implement this Nested Sampling framework, we will then proceed to introduce an algorithm, called Inner Nested Sampling. This algorithm is, for the specific problem of transforming the probability distribution $\boldsymbol{\theta}$ to a probability distribution of the function $u\!\left(\boldsymbol{\theta}\right)$, the optimal implementation, in that it takes into the specific geometry of the probability distribution (1.1).

2. Why Nested Sampling?

In order to answer question why one ought to use Nested Sampling (NS) as their Monte Carlo framework of choice, we quote Skilling himself (personal communication, 2013):

NS is designed to get the evidence value, which is equation #1 of inference.

(1)  Z = INTEGRAL dtheta JointProb(theta,data)

(2)  Posterior(theta) = JointProb(theta,data) / Z

Logically, Z precedes the posterior and should never be ignored. It’s the more important half of Bayes, not only in the algebra, but also because there’s no point in bothering to generate a posterior if the evidence falls much short of some other plausible model.

Stated differently, let $p\!\left(\theta,D\right)$ be the product of some prior $p\!\left(\theta\right)$ of some unknown set of parameters $\theta$ and some likelihood function $L\!\left(\theta,D\right)$ of some data set $D$:

[TABLE]

Then the integral

[TABLE]

is the evidence measure which may be used to differentiate between competing models, by way of Bayesian model selection. Furthermore, the posterior of the parameters $\theta$, given the data $D$, that is, $p\!\left(\left.\theta\right|D\right)$, is given as:

[TABLE]

The Nested Sampling algorithm [1, 2] is specifically designed to evaluate the integral (2.2), giving us an estimate of the evidence $Z$. Furthermore, it also provides us with a set of representative samples from the posterior (2.3), which may function as a proxy for that posterior.

For those cases where the integral (2.2) may be evaluated analytically, one will have no need for the Nested Sampling algorithm. However, for those problems where the integral (2.2) is both intractable and highly dimensional, there one will have to take his recourse to Nested Sampling, in order to be able to evaluate the evidence (2.2) and obtain a set of representative samples from the desired posterior (2.3).

The larger the evidence, the larger the adequacy of (2.1), in terms of both parsimony and data fit. Stated differently, if we are considering some $p\!\left(\theta,D\right)$, then we only need to bother to construct a posterior for the model which has the largest evidence value. For what use are our (admittedly accurate) parameter estimates, as captured in the posterior distribution, if the model that leads us to these estimates is woefully inadequate?

Before proceeding to the Nested Sampling algorithm, we will elaborate on the evidence measure and its role in Bayesian model selection. This is done for the benefit of those of the readers who are not yet well acquainted with these concepts.

3. The Evidence in Bayesian Model Selection

Bayesian statistics has four fundamental constructs, namely, the prior, the likelihood, the posterior, and the evidence. These constructs are related in the following way:

[TABLE]

As an aside, any student of Bayesian probability theory will tend to have a firm grip on the concepts of a prior, likelihood, and posterior. However, the concept of evidence is less universally known. A possible explanation for this is that most people (used to) come to Bayesian probability theory by way of the more compact relationship

[TABLE]

which does not make any explicit mention of the evidence construct111See for example [4] throughout..

In what follows, we will employ in our analysis the correct, though notationally more cumbersome, relation (3.1), and forgo of the more compact, but incomplete, Bayesian shorthand (3.2). This is done so the reader may develop some feeling for the evidence construct, and how this construct relates to the other three Bayesian constructs of prior, likelihood, and posterior.

Let $p\left(\left.\theta\right|I\right)$ be the prior of some parameter $\theta$, where $I$ is the prior information regarding the unknown $\theta$ which we have to our disposal. Let $p\left(\left.D\right|\theta,M\right)$ be the probability of the data $D$ conditional on the value of parameter $\theta$ and some likelihood model $M$ which is used222Note that $M$ is a label that points to some likelihood model, and not so much a parameter whose value we wish to determine.; the probability of the data is also known as the likelihood of the parameter $\theta$. Let $p\left(\left.\theta\right|D,M,I\right)$ be the posterior distribution of the parameter $\theta$, conditional on the data $D$, the likelihood model $M$, and the prior information model $I$. Then

[TABLE]

where

[TABLE]

is the evidence, that is, the marginalized likelihood, of both the likelihood model $M$ and the prior information model $I$. We now will show how this evidence may be used in Bayesian model selection.

If we have a set of likelihood models $M_{j}$ we wish to choose from, and just the one prior information model $I$, then we may do so by computing the evidence values $p\left(\left.D\right|M_{j},I\right)$.

Let $p\left(M_{j}\right)$ and $p\left(\left.M_{j}\right|D,I\right)$ be, respectively, the prior and posterior probability of the likelihood model $M_{j}$. Then the posterior probability distribution of these likelihood models is given as

[TABLE]

Note that if $p\left(M_{j}\right)=p\left(M_{k}\right)$ for all $j$ and $k$, then we have that (3.5) reduces to

[TABLE]

Stated differently, if we assign equal prior probabilities to our different models, then these models may be ranked by their respective evidence values $p\left(\left.D\right|M_{j},I\right)$, [5].

4. Nested Sampling: The Idea

By reducing any $k$-variate function $f$ to a corresponding monotonic descending univariate function $g$, and by using order statistics, the integral of any $k$-variate function $f$ may be evaluated using a Monte Carlo sampling scheme called Nested Sampling, [2].

Say, we wish to numerically evaluate a bivariate distribution $f\!\left(x,y\right)$, where

[TABLE]

where $-5\leq x,y\leq 5$, Figure 1.

The total volume under the curve (4.1) is given by the integral

[TABLE]

We may evaluate the integral (4.1) through brute computational force. Say, we partition the $x,y$-plane in little squares with areas $dx_{j}dy_{k}=0.25$, for $j=1,\ldots,20$ and $k=1,\ldots,20$. Then define the center of these areas as $\left(\tilde{x}_{j},\tilde{y}_{k}\right)$, and compute the strips of volume $V_{jk}$ as

[TABLE]

In Figure 2 we give all the volume elements $V_{jk}$ together.

Summating the volume elements (4.3), the total volume under the curve (4.1) may be approximated as

[TABLE]

The 3-dimensional volume elements $V_{jk}$ may be mapped to corresponding 2-dimensional area elements $A_{i}$, without any loss of generality. In order to demonstrate this, we introduce the following notation

[TABLE]

where index $i$ is a function of the indices $j$ and $k$:

[TABLE]

and $i=1,\ldots,400$. Using (4.5), we may rewrite (4.3) as

[TABLE]

In Figure 3 we give all the 400 area elements $A_{i}$ together.

Since (4.7) is equivalent to (4.3), we have that the mapping of the 3-dimensional volume elements, $V_{jk}$, to their corresponding 2-dimensional area elements, $A_{i}$, has not led to any loss of information; that is,

[TABLE]

We may, trivially, rearrange the elements $A_{i}$ in Figure 3 in descending order, Figure 4.

All these rectangular elements have a base of $dw=dx\>dy=0.25$. So, being that there are 400 area elements, we might view Figure 4 as a representation of some hypothetical monotonic descending function $g\!\left(w\right)$, where $0\leq w\leq 100$, Figure 5.

What we have accomplished is that we have mapped 3-dimensional volume elements, Figure 2, to 2-dimensional area elements, Figure 3, after which we have rearranged the area elements, Figure 4, so as to get a monotonic descending function $g\left(w\right)$, Figure 5. If we integrate the univariate function $g\left(w\right)$, Figure 5, we get the volume we are looking for.

Any $k$-variate function may be reduced to its corresponding monotonic descending univariate function $g\!\left(w\right)$. We will see that the procedure of Nested Sampling is based upon the equivalence between any $k$-variate function and its corresponding univariate representation $g\!\left(w\right)$.

4.1. Estimating abscissa’s

Skilling’s Nested Sampling framework is unique in that it is a Monte Carlo scheme that uses probability theory to its advantage, [2].

If we have a value of $g\!\left(w\right)$, without knowing the value of $w$. Then the only thing we know about $w$ is that it must lie somewhere in the region $0\leq w\leq W$, where $W$ is the scalar value of the volume on which the $k$-variate function is defined.

So, if, based on this information alone, we were asked to assign a probability distribution to the actual value of the abscissa $w$ that corresponds with the ordinate $g\!\left(w\right)$, then we would have to assign a univariately uniformly distribution on the domain $\left[0,W\right]$, that is,

[TABLE]

Consequently, the unknown abscissa value $w$ has an expectation value of

[TABLE]

and a standard deviation of

[TABLE]

Let $\boldsymbol{x}=\left(x_{1},\ldots,x_{k}\right)$ be a point in the $k$-dimensional domain. If we have $N$ random samples of the $k$-variate domain, that is, $\left\{\boldsymbol{x}_{1},\ldots,\boldsymbol{x}_{N}\right\}$. Then, by evaluating these points, we have a set of $N$ random ordinates, that is, $\left\{f\!\left(\boldsymbol{x}_{1}\right),\ldots,f\!\left(\boldsymbol{x}_{N}\right)\right\}$. Because of the equivalence between any $k$-variate function $f$ and its univariate representation $g$, we may write

[TABLE]

So, by construction, the set of $N$ random ordinates $\left\{f\!\left(\boldsymbol{x}_{1}\right),\ldots,f\!\left(\boldsymbol{x}_{N}\right)\right\}$ corresponds with the set of random ordinates $\left\{g\!\left(w_{1}\right),\ldots,g\!\left(w_{N}\right)\right\}$.

Even though we have the ordinate $g\!\left(w_{n}\right)$, by way of (4.12), we do not know the abscissa $w_{n}$ that goes with this ordinate. As we sample the $k$-variate domain $\boldsymbol{x}$, rather than the univariate domain $w$. However, the one thing we do know is that the smallest realisation of $g\left(w_{n}\right)$ corresponds with the greatest value of $w_{n}$. This is because function $g$ is a monotonic descending function, see Figure 5. It follows that we may use an order distribution for the unknown value $\max\left(w_{n}\right)$, or, for short, $w_{\max}$.

The order probability distribution for the largest draw, $w_{\max}$, in a sample of size $N$ from the uniform distribution (4.12) is

[TABLE]

where both $N$ and $W$ are known. The order distribution (4.13) has a mean of

[TABLE]

and a standard deviation of

[TABLE]

We see in (4.15) that the standard deviation, that is, our uncertainty regarding the unknown value of $w_{\max}$ falls of with a factor $N$. We will see that (4.14) and (4.15) form the backbone of the Nested Sampling framework.

Note that Skilling himself does not arrive at Nested Sampling algorithm by way of order statistics, but by way of the closely related probability distribution of the shrinkage ratio of the admissible domain of the unknown univariate representation $g$ of the multivariate function $f$:

[TABLE]

This shrinkage ratio has a probability distribution

[TABLE]

And we may compare (4.13) with (4.17), to see that both probability distributions are equal, up to a change of variable, which scales the maximum abscissa value $w_{\max}$ to the fraction $t$, (4.16).

4.2. Nested Sampling: The Algorithm

In this initial version of the Nested Sampling framework we will not protect against under- and overflow. We will just focus here on the basic philosophy which underlies Nested Sampling.

In what follows we reserve subscripts for sample membership and superscripts for the enumeration of the algorithmic step.

Step 1

Find $N$ random values $\boldsymbol{x}$ in the $k$-variate domain and evaluate, that is, find

[TABLE]

Since we may perform the steps as shown in Figure 1 through Figure 5, it holds, by way of construction, that the $N$ values of $f\!\left(\boldsymbol{x}_{n}^{\left(1\right)}\right)$ correspond with $N$ values of $g\!\left(w_{n}\right)$, (4.12).

As stated before, we cannot explicitly map the abscissa $\boldsymbol{x}$ to the corresponding abscissa of the monotonic descending function, $w$. We could link the volume elements (4.3) to the area elements(4.10) using brute computational force. But this would amount to regular numerical integration and, thus, we would still be faced with the curse of dimensionality. So, instead, we use (4.14) to statistically approximate the $w$ which corresponds with the smallest sampled value of $f\!\left(\boldsymbol{x}_{n}^{\left(1\right)}\right)$ and, thus, get our first coordinate $\left[w^{\left(1\right)},g^{\left(1\right)}\right]$ of the unknown function $g$, Figure 5. The abscissa and ordinate of this coordinate are, respectively,

[TABLE]

and

[TABLE]

As the error of the statistically approximated $w^{\left(1\right)}$ falls of with a factor $N$, (4.15), the accuracy of the approximation (4.19) will increase as $N$, the number of random samples, increases.

Finally, we approximate the integral to the right of $w^{\left(1\right)}$ of our constructed function $g$, Figure 5, as

[TABLE]

and set the evaluated integral at the first step, $Z^{\left(1\right)}$, to

[TABLE]

Step 2

We again find $N$ random values $\boldsymbol{x}$ in the $k$-variate domain and evaluate. But now we constrain the random values $\boldsymbol{x}$ in that their mapping to the function $f$ be equal or greater than the minimum value of the mapping of the previous mapping, that is, we sample the $\boldsymbol{x}_{n}^{\left(2\right)}$ under the constraint

[TABLE]

In the first iteration our state of knowledge was that $0\leq w\leq W$. In the second iteration we sample $\boldsymbol{x}_{n}^{\left(2\right)}$ under constraint (4.23), which is equivalent to sampling $w$ under the constraint, (4.20),

[TABLE]

where $g$ is a monotonic descending function in $w$. So, the updated state of knowledge is that $0\leq w\leq w^{\left(1\right)}$.

Using (4.14) again, but replacing $W$, the initial upperbound of $w$, with $w^{\left(1\right)}$, the new upperbound of $w$, we may approximate the second coordinate $\left[w^{\left(2\right)},g^{\left(2\right)}\right]$ of the unknown function $g$:

[TABLE]

and

[TABLE]

We approximate the area of the integral between $w^{\left(2\right)}$ and $w^{\left(2\right)}$, where $w^{\left(2\right)}<w^{\left(1\right)}$, as

[TABLE]

and set the evaluated integral at the second step, $Z^{\left(2\right)}$, to

[TABLE]

Step t

In iteration step $t$ we approximate coordinate $\left[w^{\left(t\right)},g^{\left(t\right)}\right]$ of the unknown function $g$ by way of

[TABLE]

and

[TABLE]

where the $\boldsymbol{x}_{n}^{\left(t\right)}$ have been sampled under the constraint

[TABLE]

The area of the integral between $w^{\left(t\right)}$ and $w^{\left(t-1\right)}$, where $w^{\left(t\right)}<w^{\left(t-1\right)}$, is approximated as

[TABLE]

The evaluated integral at the step $t$, $Z^{\left(t\right)}$, is updated as

[TABLE]

Termination step

Because of the identity (4.29)

[TABLE]

we have that $w^{\left(t\right)}\rightarrow 0$ as $t\rightarrow\infty$. If $f$ is bounded, then by construction so is $g$. It then follows from (4.32) that bounded functions $f$ or, equivalently, bounded $g$

[TABLE]

as $t\rightarrow\infty$. The maximum of $f$, the function whose integral we wish to evaluate, is also the bound of $g$, its univariate representation. This implies the following inequality

[TABLE]

By way of (4.32) and (4.35), we then have that

[TABLE]

So, if we can determine the maximum of $f$, we can use as a possible terminating condition the point where the upperbound of $A^{\left(t+1\right)}$ does not contribute more than $N^{2}$ part to $Z^{\left(t+1\right)}$

[TABLE]

Alternatively, if we are unable determine the maximum of $f$, then we take as a stopping criterium the point where $A^{\left(t+1\right)}$ itself does not contribute more than $N^{2}$ part to $Z^{\left(t+1\right)}$

[TABLE]

However, (4.38) comes with the caveat that the function $f$, or, equivalently, $g$, might still have some unexplored regions having large enough values $w^{\left(q+1\right)}$, where $q>t$, to take

[TABLE]

which would imply a premature termination of the algorithm.

4.3. Optimal implementation

If one is to implement the steps in the previous section one would get an algorithm very differently from the pseudo-code given by Skilling, [2]. These differences are mainly differences in implementation.

The reason for us to first give the naive algorithm was to point to reader to the elegant and simple idea behind Nested Sampling, without losing ourselves too much in the technicalities of optimal implementations.

However, with the core idea behind Nested Sampling demonstrated, we now will treat in the next two paragraphs the points of optimal implementation.

4.3.1. Reduction of computational cost

In the previous treatment of the Nested Sampling framework, $N$ new samples of $\boldsymbol{x}$ were drawn at each iteration step $t$, under the constraint

[TABLE]

Now, instead of drawing $N$ new samples at each sampling step $t$, we may also realize that in iteration $\left(t-1\right)$, we already had $\left(N-1\right)$ samples of $f$ at our disposal that satisified the constraint (4.40). Namely, those samples in iteration step $\left(t-1\right)$ smaller than

[TABLE]

If we keep these $\left(N-1\right)$ samples, then we only need to sample one aditional ordinate value which satisfies the constraint (4.40), in order to obtain our needed sample of $N$ objects.

So, after each iteration $t$ we discard one object from our sample of $N$ objects. This discarded object becomes $g^{\left(t\right)}$. The $\left(N-1\right)$ surviving objects are taken to the next iteration and an additional object is sampled under constraint

[TABLE]

This implementation reduces the computational costs of Nested Sampling with an order of magnitude of $N$.

4.3.2. Guarding against under- and overflow

In many problems $Z^{\left(t\right)}$ may become so large that computational overflow may occur, that is, there is no longer a floating number representation possible for its value. To remedy this problem we will have to work with the $\log Z^{\left(t\right)}$. Furthermore, We also have that certain functions $f$ have certain values in their domain so small that computational underflow may occur. So, in what follows we will evaluate $\log f$ instead of $f$. Likewise, because of the fact that $w_{t}\rightarrow 0$ as $t\rightarrow\infty$, we have that for sufficiently large $t$ computational underflow may occur, that is, $w_{t}$ may become so small that there is no longer a floating number representation possible for its value. To remedy the latter situation we will go from the $w$ scale to the $\log w$ scale.

To go to the $u=\log w$ scale, we will have to make a proper change of variable for the order distribution (4.13). We have that

[TABLE]

Substituting (4.41) into (4.13), we may obtain

[TABLE]

for $-\infty<u_{\max}\leq\log W$. with mean standard deviation of

[TABLE]

and

[TABLE]

With a repeated application of (4.43) we may find the limits of the $u$ scale after the $t$th iteration to be

[TABLE]

From (4.45), it follows that that the width of the $t$th interval on the original $w$ scale may be written as

[TABLE]

Consequently, we have

[TABLE]

We then have that that the log of the area element at the $t$th iteration is the sum

[TABLE]

With (4.3.2) the computation of $A^{\left(t\right)}$ is sufficiently protected from underflow of both $w^{\left(t\right)}$ and $f$, or, equivalently, $g$.

We may update $\log Z^{\left(t\right)}$ with $\log A^{\left(t\right)}$ using the formula for logarithmic addition, (Skilling, 2006),

[TABLE]

If we set $x=\log Z^{\left(t-1\right)}$ and $y=\log A^{\left(t\right)}$, then (4.3.2) gives us

[TABLE]

With (4.50) the computation of $\log Z^{\left(t\right)}$ is sufficiently protected from the overflow of $\log Z^{\left(t-1\right)}$.

To summarize, in order to protect the algorithm from under- and overflow, which may easily occur in actual problems, the original algorithm remains unchanged, except that each iteration $\log A^{\left(t\right)}$ is computed, instead of $A^{\left(t\right)}$, by way of (4.3.2), and $\log Z^{\left(t\right)}$ is updated, instead of $Z^{\left(t\right)}$, by (4.50). Note that the termination conditions (4.37) and (4.38) transform, respectively, to

[TABLE]

if we know the maximum of the function $f$, and

[TABLE]

if we do not.

5. Applying Nested Sampling to Unscaled Posteriors

Nested Sampling is a general purpose algorithm for the evaluation of the integrals of multivariate functions. A special class of multivariate functions are unscaled posterior distributions.

Let $p\!\left(x,y\right)$ be some prior distribution of the unknown parameters $x$ and $y$. Let $p\!\left(\left.D\right|x,y\right)$ be the probability of some observed data set $D$ given the parameters $x$ and $y$, or, equivalently, the likelihood $L\!\left(x,y\right)$. Then the unscaled posterior distribution is

[TABLE]

If we integrate $p\!\left(x,y,D\right)$ over the domain of $x$ and $y$ we get the marginal probability of the observed data set $D$

[TABLE]

The marginal probility $p\!\left(D\right)$ is also called the evidence. Scaling $p\!\left(x,y,D\right)$ with $p\!\left(D\right)$ we get, by way of the product rule, the scaled posterior distribution:

[TABLE]

Furthermore, the evidence is an important quantity in its own right. It may be used for model selection. Say, we have two competing likelihood models, $L_{1}$ and $L_{2}$. Combining (5.1) and (5.2), we may compete the respective evidences of the competing likelihood models:

[TABLE]

Then the likelihood model that fits the data best will be the model that gives us the greatest marginal probability of obtaining the observed data set $D$. The evaluation of the integral (5.1) is a non-trivial matter, especially for highly variate likelihood functions. It was with this application in mind that Skilling developped his Nested Sampling framework.

We now proceed to apply the Nested Sampling framework to (5.1). Let the function $f$ be

[TABLE]

Let the prior $p\!\left(x,y\right)$ be the uniform distribution over the domain of $x$ and $y$. If this domain has a total area of $W$, then we have that

[TABLE]

Substituting (5.6) into (5.5), and substituting the resulting function $f$ into (4.3.2), we get

[TABLE]

where $\min L_{n}^{\left(t\right)}$ is the smallest likelihood from the sample of $N$ likelihoods at the the $t$th step, which were sampled under the constraint that

[TABLE]

The Nested Sampling framework is most often presented in the form of (5), which is the special case of $f$ being an unscaled posterior wich takes as its prior the uniform distribution of the parameter space. However, the framework is much more general in that is applicable to any function $f$, not just unscaled posterior distributions.

Last but not least, if $f$ is an unscaled posterior, then the Nested Sampling framework not only evaluates the evidence. It also gives us a Monte Carlo proxy of the posterior of interest. Let $T$ be the termination step of the Nested Sampling run. Let $\boldsymbol{x}^{\left(t\right)}$ be the points in the $k$-variate domain that corresponded with the $\boldsymbol{x}_{n}^{\left(t\right)}$ for which

[TABLE]

where $t=1,\ldots,T$ and $n=1,\ldots,N$. Then we may assign to these $\boldsymbol{x}^{\left(t\right)}$ the probability weights, (4.32) and (4.33),

[TABLE]

This leaves us with a set of weighted random Monte Carlo samples. The weighted samples are a proxy for the posterior distribution of interest, [2] .

Let the posterior $p\!\left(\left.\boldsymbol{x}\right|D\right)$ be given as

[TABLE]

Let $h$ be some function defined on the $k$-variate domain of the $\boldsymbol{x}$. Then we may approximate the weighted $h$ by way of Nested Sampling as, (5.10),

[TABLE]

It is (5.12) that will us enable us to compute the first and second moments of the relevances defined on a Dirichlet distribution.

6. The Issue of Confidence Bounds

Say we wish to assign the function $u$ to a bivariate probability distributions of the form:

Then, in practice, we only have indirectly access, by way of our data, to the probability distribution of interest; that is, rather then the known probabilities, Table 1, we only have some observed count data, Table 2:

where $n$ is the total number of observations.

Let $D=\left(r_{1},r_{2},r_{3},r_{4}\right)$ be the observed count data in Table 2. Let $\boldsymbol{\theta}=\left(\theta_{1},\theta_{2},\theta_{3},\theta_{4}\right)$ be the vector of the unknown probabilities in Table 1. Then the likelihood function of the unknown probabilities is assumed to follow multinomial distribution:

[TABLE]

As a prior for the unknown probabilities we assign the uninformative Dirichlet prior

[TABLE]

which, if marginalized, collapses to the uninformative Beta prior. Combining the likelihood (6.1) and prior(6.2), by way of the product rule, and normalizing, by way of the sum rule, we obtain the multivariate Dirichlet posterior distribution of the theta’s given the observed count data:

[TABLE]

Each realization of the function $u\!\left(\boldsymbol{\theta}\right)$ maps onto a corresponding probability $p\!\left(\left.\boldsymbol{\theta}\right|D\right)$. By arranging the values $u\!\left(\boldsymbol{\theta}\right)$ on the $x$-axis and the corresponding $p\!\left(\left.\boldsymbol{\theta}\right|D\right)d\boldsymbol{\theta}$ on the $y$-axis, we may obtain the univariate probability distribution of the function $u\!\left(\boldsymbol{\theta}\right)$. This probability distribution of $u\!\left(\boldsymbol{\theta}\right)$ takes the uncertainty into account we have in regards to the unknown $\boldsymbol{\theta}$ and, consequently, lets us put confidence bounds on this function.

If we only have four unknown theta’s we may use brute computational force to partition the domain of $\boldsymbol{\theta}$ and compute of for each partitioning the corresponding pair $\left[u\!\left(\boldsymbol{\theta}\right),p\!\left(\left.\boldsymbol{\theta}\right|D\right)d\boldsymbol{\theta}\right]$, after which we then order the $u\left(\boldsymbol{\theta}\right)$ and plot them together with their corresponding probabilities $p\left(\left.\boldsymbol{\theta}\right|D\right)d\boldsymbol{\theta}$. However, for large distributions, having many unknown $\boldsymbol{\theta}$’s, this quickly becomes unpractical because of the curse of dimensionality.

In the next section we present an implementation of the Nested Sampling framework by Skilling that allows us to evaluate the probability distribution of $u\left(\boldsymbol{\theta}\right)$ for highly variate distributions of $\boldsymbol{\theta}$, by way of Monte Carlo sampling. The implementation is accomplished by way of the Inner Nested Sampling algorithm. The Inner Nested Sampling algorithm allows us to sample uniformly from the constrained likelihood space of the $\boldsymbol{\theta}$’s, a necessary prerequisite of the Nested Sampling framework.

7. Inner Nested Sampling

Let $L$ be a likelihood function defined on a highly multivariate parameter space $\boldsymbol{\theta}$. Then Nested Sampling is a Monte Carlo framework with which this multivariate likelihood function $L\!\left(\boldsymbol{\theta}\right)$ may be evaluated. The Nested Sampling framework needs uniform samples within the multivariate geometry of some likelihood constraint $L^{*}$ in order to work. However, this framework does not tell us how to obtain these samples, that is, its optimal implementation is an open-ended research question. In this paper we give an algorithm, called Inner Nested Sampling, that obtains such uniform samples.

The idea behind Inner Nested Sampling is that we obtain a set of differentials of the multivariate geometry of the initial likelihood constraint $L^{*}$ at iteration step $t=0$ of Nested Sampling proper. These differentials are defined by a direction $\mathbf{e}$ and a radius $R\!\left(\mathbf{e}\right)$ and serve as a proxy for the actual geometry and have the nice property that they may be uniformly sampled. Furthermore as with each iteration step $t$ the geometry defined by likelihood constraint will shrinks, the radii $R\!\left(\mathbf{e}\right)$ may be updated so as to reflect this shrinkage. This then allows us to continue the uniform sampling of these differentials and, by proxy, the likelihood geometry of interest.

7.1. A change of variable

Say, we have a $m$-variate parameter vector $\boldsymbol{\theta}=\left(\theta_{1},\ldots,\theta_{m}\right)$. Then the likelihood constraint $L\!\left(\boldsymbol{\theta}\right)\leq L^{*}$ defines some sub-domain $V^{*}$ of the total paramer space $V$ of $\boldsymbol{\theta}$.

Now, if we have an (approximate) modus $\boldsymbol{\hat{\theta}}$ of the likelihood function $L\!\left(\boldsymbol{\theta}\right)$ we may translate the origin of $\boldsymbol{\theta}$ to the location of this modus and proceed to make a change of variables from a Cartesian coordinate system to unit-vector coordinates,

[TABLE]

where $\mathbf{e}=\left(e_{1},\ldots,e_{m}\right)$ is a point on the unit-sphere and $r$ is the distance from the $\boldsymbol{\hat{\theta}}$ to the likelihood constraint $L^{*}$ in the direction of $\mathbf{e}$. Unit-vector coordinates map the parameter vector $\boldsymbol{\theta}$ to a radius $r$ and $\left(m-1\right)$ non-redundant coordinates:

[TABLE]

where the $m$th redundant coordinate of the unit vector $\mathbf{e}$ may be found through the identity

[TABLE]

The differential of the unit-vector transformation is

[TABLE]

where $dS$ are equi-volume ‘patches’ of the surface $S$ of the $m$-dimensional unit-sphere. Note that the surface $S$ itself has dimensionality $\left(m-1\right)$. The differentials (7.3) are the volumes of pyramids with base $dS$ and height $R\!\left(e_{1}\cdots e_{m}\right)$. So, the integral

[TABLE]

For example, say, we have the likelihood

[TABLE]

which is given in Figure 6. The likelihood constraint $L^{*}=0.041$ corresponds with an ellips having ‘volume’ $V^{*}=8.162$, Figure 7.

Since we have a two-dimensional likelihood space, the unit sphere is actually an unit circle and it follows that the surface $S$ is actually a circumference, where $S=2\pi$. The $dS$ is obtained by partitioning the circumference $S$ in $n$ equi-distant line elements, that is $dS=2\pi/n$. Now, if we let $n=16$, then

[TABLE]

and the radii $R_{i}$, for $i=1,\ldots,n$ at the centers of these line elements $dS$ are shown in Figure 8.

By summating the approximate differentials, (7.3) and (7.6),

[TABLE]

we may approximate the integral (7.4):

[TABLE]

As we let $n\rightarrow\infty$, this approximation will become evermore accurate.

We summarize, the constrained likelihood space $V^{*}$ may be represented as a collection of differentials $\left\{dV^{*}_{1},\ldots,dV^{*}_{n}\right\}$. As it will turn out, it is trivially simple to sample uniformly from the set of the differentials $dV^{*}_{i}$, and by doing so we actually sample uniformly the space $V^{*}$ itself as we let $n\rightarrow\infty$, (7.4).

7.2. Obtaining differentials

The differentials depicted in Figure 8 were obtained through brute computational force. For highly multiariate likelihood spaces such methods are bound to fail due to the curse of dimensionality. However, Skilling’s Nested Sampling framework comes here to the rescue, as we shall now demonstrate.

In analogy to the treatment of Nested Sampling proper in Chapter 6, we partion $dS$ in $n=400$ equi-distant line elements, compute the area of the corresponding differentials $\left(dV_{1}^{*},\dots,dV_{400}^{*}\right)$, and plot them, Figure 9.

Again we are free to reorder these area elements as we like, Figure 10.

Now all these rectangular elements have a base of $dw=dS=2\pi/400$. Being that there are 400 area elements we again might view Figure 11 as a representation of some monotonic descending function $g\left(w\right)$, where $0\leq w\leq 2\pi$, Figure 11.

Were we to integrate this function $g\left(w\right)$, we obtain an approximation of the area $V^{*}$, (7.8). However, this approximation is not our primary interest in Inner Nested Sampling. Though it may serve a check for those rare cases were we actually know the value of $V^{*}$. Rather, the collection of differentials is what we are after, since these may sampled uniformly, thus, allowing us, by proxy, to uniformly sample $V^{*}$ itself.

The Nested Sampling framework, previously discussed, translates directly to Inner Nested Sampling, were we set $W$ to be the surface $S$ of the $m$-dimensional unit sphere

[TABLE]

and take as the function to be evaluated

[TABLE]

In doing so, we end up with a collection of directions, $\mathbf{e}^{\left(q\right)}$, for $q=1,\ldots,Q$, with associated radii $R\!\left(\mathbf{e}_{q}\right)$, which are the distances from the modus $\boldsymbol{\hat{\theta}}$ to the constraint $L^{*}$, and weights $A^{\left(q\right)}$, (4.29), (4.32), and (7.10):

[TABLE]

where $N$ is the number of objects used in the Inner Nested Sampling run and the $N$ directions $\mathbf{e}$ have been sampled under the constraint

[TABLE]

The weights $A^{\left(q\right)}$ correspond with the volumes $dV_{q}^{*}$. By substituting (7.9) in (7.2), we find, (7.3):

[TABLE]

where

[TABLE]

Note that for the differentials $dV^{*}$ and $dS$ we have let the subscripts ennumerate the iteration step.

8. Obtaining uniform samples for Inner Nested Sampling

The vector $\boldsymbol{\theta}=\left(\theta_{1},\dots,\theta_{m}\right)$ is a point in the parameter space and $L\left(\boldsymbol{\theta}\right)$ is the likelihood function that is defined on this parameter space. In highly dimensional problems the region of interest of $L\left(\boldsymbol{\theta}\right)$ becomes exponentially small relative to the total sample space. Therefore, the change that random sampling will yield a non-negligible likelihood value becomes exponentially small as the dimensionality of problem grows.

In Inner Nested Sampling not the likelihood function $L\left(\boldsymbol{\theta}\right)$ is evaluated but the radius function $R\left(\mathbf{e}\right)$, where $\mathbf{e}$ is a point on the $m$-dimensional unit sphere. The radius function $R\left(\mathbf{e}\right)$ is more uniformly spread out over the surface of the unit-sphere than $L\left(\boldsymbol{\theta}\right)$ is over the parameter space. This makes random sampling of the radius function much more feasible. Nonetheless, as the bound $R^{*}$ on $R\left(\mathbf{e}\right)$ becomes tighter and tighter, the number of rejections of the random draws will grow so much as to make a random walk necessary.

8.0.1. Random sampling of the unit sphere

Random points on the unit-sphere are easily obtained by letting

[TABLE]

A rejection occurs whenever $R\left(\mathbf{e}\right)<R^{*}$. If the number of rejections of the random draws of the $\mathbf{e}$ exceeds some number $M$, where $M>>1$, a random walk on the unit sphere is performed.

8.0.2. Random walk on the unit sphere

The random walk on the unit sphere takes at its starting point some $\mathbf{e}_{k}$ that has a radius $R_{k}$ which is known to satisfy the radius constraint $R_{k}>R^{*}$. and we set

[TABLE]

If the parameter space is 3-dimensional, then $\mathbf{e}^{\left(0\right)}$ will be a point on the surface of a 3-dimensional sphere. Then, for a given length $s$ of the random step to be taken, there is a 2-dimensional circle lying on the surface of this sphere with equi-probable possible candidates $\mathbf{e}^{\left(1\right)}$. To obtain a realisation $\mathbf{e}^{\left(1\right)}$, we sample a random point $\mathbf{u}$ on this circle, where

[TABLE]

The circumference of the 2-dimensional circle whith the candidates $\mathbf{e}^{\left(1\right)}$ will change as the length $s$ of the random step changes. To be more more precise, let $\mathbf{e}^{\left(1\right)}$ be the proposal point of the first random step and let $\alpha$ be the angle between the starting point $\mathbf{e}^{\left(0\right)}$ and $\mathbf{e}^{\left(1\right)}$. Then the radius of the 2-dimensional circle will equal $\sin\alpha$, where $0\leq\alpha\leq\pi$.

If we translate the origin of our original coordinate system to $\left(\cos\alpha\right)\mathbf{e}^{\left(0\right)}$ and then rotate the axes of this 3-dimensional system so that the 2-dimensional circle, on which all the $\mathbf{e}^{\left(1\right)}$ are to be found, lies in the plane spanned by axes 2 and 3. Then the coordinates of $\mathbf{e}^{\left(1\right)}$ in this new coordinate system become, (8.3):

[TABLE]

In order to transform the coordinates of this alternative coordinate system, $\tilde{\mathbf{e}}^{\left(1\right)}$, back into the coordinates of the original coordinate system, $\tilde{\mathbf{e}}^{\left(1\right)}$, we must realize that in our original coordinate system the 2-dimensional circle of the $\tilde{\mathbf{e}}^{\left(1\right)}$ is orthogonal to $\left(\cos\alpha\right)\mathbf{e}^{\left(0\right)}$. This implies that the rotation matrix $B$, which accomplishes this re-transformation, has as its first column:

[TABLE]

The other two orthonormal columns of $B$, $\mathbf{b}_{2}$ and $\mathbf{b}_{3}$, then can be easily constructed by applying a Gram-Schmidt process. Let

[TABLE]

Then

[TABLE]

and

[TABLE]

Likewise, we have

[TABLE]

and

[TABLE]

This results in the rotation matrix

[TABLE]

We then have that

[TABLE]

In the random walk the proposed step $\mathbf{e}^{\left(t\right)}$ is rejected whenever either, (7.1),

[TABLE]

or

[TABLE]

The random walk step must be modulated in order that the angle $\alpha$ becomes smaller as the number of rejections increases and, visa versa, larger as the number acceptances increases. Now, since a step length of $s$ corresponds with some angle $\alpha$, we are free to modulate $\alpha$ instead of $s$. So, let

[TABLE]

where $\tau\geq t$ and $u^{\left(\tau\right)}$ should not equal 0 or 1. Let

[TABLE]

where $step^{\left(0\right)}=0$. Then, constraining $0<\alpha^{\left(\tau\right)}<\frac{\pi}{2}$, we let

[TABLE]

where it is understood that $u^{\left(\tau\right)}$, (8.13), should not equal 0 or 1.

For the case where the dimensionality of the parameter space is $m=3$, this procedure generalizes to

[TABLE]

where $t$ is the number of random steps accepted in the random walk, $\tau$ is the total number of iterations performed by the random walk, and, (8.4),

[TABLE]

where, (8.3),

[TABLE]

In our preliminary empirical studies the number of random walk acceptances fluctuated around the $67\%$, which is as it should be. Note that this modulation of a random walk on the unit sphere is an adaptation of the modulating algorithm on the Cartesian parameter space as given by Skilling, [1].

8.1. Drawing Uniform samples from Inner Nested Sampling realizations

If for some $m$-variate parameter space we have a collection of sampled unit-vectors $\left\{\mathbf{e}^{\left(1\right)},\ldots,\mathbf{e}^{\left(Q\right)}\right\}$ with corresponding radii $\left\{R^{\left(1\right)},\ldots,R^{\left(Q\right)}\right\}$. Then this allows us to draw a uniform sample from the corresponding differentials $\left\{dV_{1},\ldots,dV_{Q}\right\}$, and by doing so we actually sample uniformly the space $V^{*}$ itself as we let $Q\rightarrow\infty$, (7.4).

Step 1

First we uniformly draw a differential $dV_{q}$. We do this by drawing from the uniform distribution:

[TABLE]

where Z is understood to be the approximation of the $V^{*}$, (4.33) and (7.2):

[TABLE]

We then find the smallest index value $q$ for which the following inequality holds, (7.8):

[TABLE]

This $q$ then is the index number of the uniformly drawn $dV_{q}$.

Step 2

Let

[TABLE]

where $0\leq v\leq A^{\left(q\right)}$. Since the differentials $dV_{q}$ are pyramids in the limit $n\rightarrow\infty$, Figure 8, the realisation $v$, $0\leq\nu\leq dV_{q}$, geometrically corresponds with the volume of a sub-pyramid having height $\rho$, where $0\leq\rho\leq R_{q}$, that is

[TABLE]

where (7.14):

[TABLE]

From (8.25), we then have

[TABLE]

Step 3 The proposal of the uniformly sampled constrained likelihood space of $\boldsymbol{\theta}$, then simply becomes, (7.1):

[TABLE]

where $\boldsymbol{\hat{\theta}}$ is the modus of the likelihood function which is evaluated, $\rho$ is (8.26), and $\mathbf{e}^{\left(q\right)}$ is the direction for which the radius $R^{\left(q\right)}$ was determined.

Integrals are limit cases of simple summation. In calculus differentials go to zero in order to obtain infinite precision. So, integrals in general have an infinite amount of differentials. Nested Sampling steps away from the limit case of calculus by approximating an infinite amount of differentials that go to zero with a finite amount of differentials greater than zero, which it then summates. Now, where we to sample to uniformly from $V^{*}$ then, because of identity (7.4), we may just as well sample uniformly over the infinite amount of differentials $dV^{*}$.

Inner Nested Sampling, just like Nested Sampling, also steps away this limit case of infinite precision. By analogy, the finite amount of differentials that were obtained through random sampling are uniformly sampled, rather than the actual infinite amount of differentials. The finite set of differentials, obtained through random sampling, serves as a proxy for the infinite set of differential. So, by sampling uniformly over the finite set of differentials one approximately samples over $V^{*}$. The larger the set of differentials, the more accurate this approximation.

8.2. Implementing Nested Sampling

We now give the algorithm that integrates the Inner Nested Sampling within the Nested Sampling framework.

Inner Nested Sampling

Let the Inner Nested Sampling do a run for some large likelihood constraint that is found by shooting a large number of points in some preset parameter space of $\boldsymbol{\theta}$. The smallest likelihood value is taken as the initial likelihood constraint $L^{*}_{0}$. Let the Inner Nested Sampling algorithm explore the geometry defined by this $L^{*}_{0}$, using $M$ objects. Then we have an approximation of the corresponding domain $V^{*}$, in the form of the weights $A^{\left(q\right)}$, (7.13), and the associated directions $\mathbf{e}^{\left(q\right)}$ and corresponding radii $R^{\left(q\right)}=R\!\left(\mathbf{e}^{\left(q\right)}\right)$.

Nested Sampling proper

Using the procedure in Section 8.1, draw $\left(N-1\right)$ random objects with likelihoods greater then $L^{*}_{0}$ from the collection of Inner Nested Sampling differentials $\left\{dV_{1}^{*},\ldots,dV_{Q}^{*}\right\}$ and add these objects with the object corresponding with the initial likeklihood constraint $L^{*}_{0}$. We then have $N$ objects with which to perform Nested Sampling proper, as described in Chapter 6.

As the likelihood contour shrinks with every further Nested Sampling iteration, the probability grows that the proposed $\boldsymbol{\theta}$, (8.27), will be rejected. When the ratio of rejections becomes prohibitively large we update the monotonic function $g$ of the Inner Nested Sampling run, Figure 11. This can be done by either determining the new $R^{\left(q\right)}$ corresponding with each stored direction vector $\mathbf{e}^{\left(q\right)}$. Or, by taking advantage of the fact that the likelihood contour generally will retain its shape as it shrinks, one alternatively may, if pressed for time, take with intervals values of $R^{\left(q\right)}$ on the $w$-‘axis’ and use some interpolating scheme for the other $R^{\left(q\right)}$ values. And sampling for Nested Sampling proper then continues for this updated function of $g$ of the Inner Nested Sampling run.

Note that the number of possible directions in the Nested Sampling algorithm is constrained to the number of vectors stored in the Inner Nested Sampling run $\mathbf{e}$. In order to have many possible directions $\mathbf{e}$, a large number of objects may be taken be taken in the Inner Nested Sampling run. In our preliminary study we found that for a 20-variate normal likelihood $M=1000$ initial Inner Nested Sampling objects resulted in $Q=12.000$ stored direction vectors $\mathbf{e}$, from which $T=2000$ random points were actually sampled in the Nested Sampling proper, which was initialized with $N=100$ objects.

8.3. A caveat

Before we finish the outline of Inner Nested Sampling algorithm we have one more point to make. In Section 8.1, we sampled uniformly from the differential $dV_{q}^{*}$ by sampling the corresponding radius $R^{\left(q\right)}$. This correspond geometrically with sampling the center line of the differential $dV_{q}^{*}$. Being that this center line itself is sampled uniformly from the unit sphere, either by way of randam draw or a random walk, this will still constitute a random draw in case this differential $dV_{q}^{*}$ is only sampled once during Nested Sampling proper. However, in the case were we draw the same differential twice, as we sample uniformly over the collection of differentials, uniformity is compromised, in that given a previous draw over that differential, we have pertinent knowledge of the possible realisations of the samples points.

This may be remedied as follows. At the preliminary phase of Inner Nested Sampling algorithm run, as the radius constraint is not yet so strong that a random walk has to be employed, a very great amount of rejected directions $e$ are generated through pure random sampling. Now, instead of throwing these rejections away, we can store a predesignated number $q$ of the rejections which have radii closest to the radii $R^{\left(q\right)}$ already in place in a ‘safety matrix’. This safety matrix is updated as the number of rejections, thus, getting radii which are evermore closer the $R^{\left(q\right)}$ already in place.

Now, it is expected that this procedure will go a long way to remedy the possible distortion of uniformity in the sampling within the differentials themselves, since it is especially these differentials $dV_{q}^{*}$ in the preliminary phase of Inner Nested Sampling which are most likely to be sampled more than once. As can be deduced from (4.34):

[TABLE]

as $q\rightarrow\infty$. However, should it happen that one of the differentials $dV_{q}^{*}$ obtained through a random walk is sampled more than once, however unlikely, two choices are left to obtain extra samples for this differentials. Either try to obtain a $\tilde{R}^{\left(q\right)}$ by performing a random walk under the contraint

[TABLE]

whenever in the uniform sampling of the differentials a $R^{\left(q\right)}$ is drawn which already has been drawn once. Or, alternatively, ignore possible distortion in uniformity, especially since the $dV_{q}^{*}$ more and more become like true differentials, that is, infinitely thin slices of volume, as $q\rightarrow\infty$, (8.28). And the ever mounting computational cost in exploring the radius constraint (8.29) will eventually not outweigh the added benefits that the object $\tilde{R}^{\left(q\right)}$ has to offer in terms of the restoration of uniformity, since the distortions in uniformity are so minute.

9. Resolving the Issue of Confidence Bounds

Let $\theta_{ij}=\left(\theta_{11},\theta_{12},\ldots,\theta_{IJ}\right)$ for $i=1,\dots,I$, $j=1,\dots,J$, where $I\times J=M$. Also let $\theta_{i+}=\left(\theta_{1+},\ldots,\theta_{I+}\right)$ and $\theta_{+j}=\left(\theta_{+1},\ldots,\theta_{+J}\right)$. The issue of interest were will apply the Inner Nested Sampling is to evaluate the value $u\left(\theta_{ij}\right)$. This evaluation is over all plausible values of the $\theta_{ij}$, weighted according to the multivariate posterior distribution of these theta’s, (6.3):

[TABLE]

where $D=\left(r_{11},r_{12},\dots,r_{IJ}\right)$.

In the case of a large $M$ the brute computational evaluation of the integrals which give us, respectively, the first and second momemt of the function $u\left(\theta_{ij}\right)$:

[TABLE]

and

[TABLE]

quickly becomes unpractical because of the curse of dimensionality. So, the Nested Sampling framework, as implemented by the Inner Nested Sampling algorithm, will be used to evaluate (9.2). This is done by evaluating, respectively, the sums, (8.19):

[TABLE]

and

[TABLE]

By way (9.4) and (9.5), we may put the following approximate confidence bounds on the relevance $u$:

[TABLE]

9.1. Inner Nested Sampling for the Dirichlet Distribution

As it turns out, the dirichlet distribution, that is, the multinomial beta distribution, (9.1) is particullary amenable to the Inner Nested Sampling algorithm, since it has known modus as well as analytical derivable $R\!\left(\mathbf{e}\right)$ for the first run, when no constraints radius constraint is yet in place, other than the simplex form of the domain on which (9.1) is defined. This is particullary helpful, since in general we have that the finding of the radii $R\!\left(\mathbf{e}\right)$ is the most computational intensive part of the Inner Nested Sampling, as it requires some kind of search algorithm, such as slice sampling, for example.

We now revert back to the notation of Section 6, and let

[TABLE]

and

[TABLE]

The modus of (9.8), using some mathematical package, may be derived as:

[TABLE]

where $n=\sum_{i}r_{i}$. The radius constraint, for any modus $\hat{\boldsymbol{\theta}}$, is the distance from this modus to the boundary of the $\left(M-1\right)$-simplex defined by the vertices

[TABLE]

where the $\left(M-1\right)$-simplex is the $M$-dimensional domain of the Dirichlet parameters $\boldsymbol{\theta}$.

The $\left(M-1\right)$-simplex is spanned by the $\left(M-1\right)$ vectors, (9.1),

[TABLE]

The coordinate system defined by (9.1) is not orthonormal. And in order to make it so, we first orthogonalize the system by way of a Gramm-Schmidt process. First, arbitrily, let

[TABLE]

Then

[TABLE]

and

[TABLE]

etc., until $\mathbf{b}_{M-1}$. Then, normalizing the orthogonal vectors (9.14),

[TABLE]

we obtain the $\left(M-1\right)$-dimensional orthonormal coordinate system $W$ we are looking for:

[TABLE]

Let $\mathbf{e}$ be a random point on the $M-1$-dimensional unit-sphere, obtained by letting, (8.1),

[TABLE]

And let

[TABLE]

Then any point

[TABLE]

will lie on the same ‘plane’ as the $\left(M-1\right)$-simplex. In order to get the maximum radius $R$ for which $\mathbf{p}$ is a point on the the bondary of the $\left(M-1\right)$-simplex, we must realize that on this bondary one of the elements of the $\mathbf{p}$ must necessarily be zero. So we solve for every element of $\mathbf{p}$:

[TABLE]

and obtain

[TABLE]

As only positive radii $R_{i}$ are allowed, the maximum radius $R$ for which $\mathbf{p}$ still lies on the boundary of the $\left(M-1\right)$-simplex is

[TABLE]

As any radius greater than $R$ will automtically result in negative $p_{i}$ values, which then would imply that $\mathbf{p}$ lies outside of the $\left(M-1\right)$-simplex.

We summarize, in order to use Inner Nested Sampling on the $M$-dimensional Dirichlet distribution (9.8), we draw a random point on the $\left(M-1\right)$-dimensional unit-sphere, (9.17), which we then transform to the direction (9.18). And the corresponding distance $R$ from the modus $\hat{\boldsymbol{\theta}}$ to the boundary of the domain of the Dirichlet parameters, $\boldsymbol{\theta}$, which is the $\left(M-1\right)$-simplex defined by the vertices (9.1), is found by way of (9.21) and (9.22). Furthermore, the volume of this $\left(M-1\right)$-simplex is known to be

[TABLE]

and (9.23) then may serve as check on the accuracy of the estimation, (8.19),

[TABLE]