# Nested Sampling with Normalising Flows for Gravitational-Wave Inference

Michael J. Williams  John Veitch  Chris Messenger  SUPA, School of Physics and
Astronomy
University of Glasgow
Glasgow G12 8QQ, United Kingdom

###### Abstract

We present a novel method for sampling iso-likelihood contours in nested
sampling using a type of machine learning algorithm known as normalising flows
and incorporate it into our sampler Nessai. Nessai is designed for problems
where computing the likelihood is computationally expensive and therefore the
cost of training a normalising flow is offset by the overall reduction in the
number of likelihood evaluations. We validate our sampler on 128 simulated
gravitational wave signals from compact binary coalescence and show that it
produces unbiased estimates of the system parameters. Subsequently, we compare
our results to those obtained with dynesty and find good agreement between the
computed log-evidences whilst requiring 2.32 times fewer likelihood
evaluations. We also highlight how the likelihood evaluation can be
parallelised in Nessai without any modifications to the algorithm. Finally, we
outline diagnostics included in Nessai and how these can be used to tune the
sampler’s settings.

††preprint: APS/123-QED

##  I Introduction

Gravitational-wave astronomy has contributed to our understanding of physics
and astrophysics from the atomic scale up to the scale of the Universe [1, 2,
3]. This is set to continue as the detectors of the LIGO-Virgo-Kagra (LVK)
Collaboration continue to improve in sensitivity [4] and new detectors come
online [5] increasing the potential of detecting previously unseen types of
sources.

Prior to the recent third observing run there were 11 confirmed detections of
gravitational waves from compact binary coalescence (CBC) [6], including the
first multi-messenger detection [7]. The first half of the last year-long
observing run has resulted in a further 39 candidate detections [8]. Each of
the candidate events requires extensive analysis to determine whether they are
astrophysical in origin and consequently allow us to understand the
astrophysics that produce such phenomena. The nature of these analyses and
ever increasing number of candidates necessitates more efficient analysis
techniques.

Our understanding of the sources that produce the detected GWs hinges on the
ability to infer the parameters that describe them. This inference is carried
out in a Bayesian framework [9, 10] centred around Bayes’ theorem which allows
for prior knowledge to be updated with observations to obtain posterior
distributions, or more explicitly:

$p(\bm{\theta}|d,H)=\frac{p(d|\bm{\theta},H)p(\bm{\theta}|H)}{p(d|H)},$ (1)

where $\bm{\theta}$ are the parameters that describe the model $H$, $d$ is
observed data, $p(\bm{\theta}|d,H)$ is the posterior, $(d|\bm{\theta},H)$ is
the likelihood, $p(\bm{\theta}|H)$ is the prior and $p(d|H)$ is the evidence.
Whilst computing the posterior is trivial in lower dimensions, the
gravitational-wave signals from compact binary coalescence are described by a
minimum of 15 parameters and the resulting parameter space can be multi-modal
and highly correlated [11]. Typically, this requires applying stochastic
sampling techniques such as Markov Chain Monte Carlo (MCMC) [12] and Nested
Sampling [13]. These techniques are computationally expensive and their cost
is directly related to the cost of evaluating the likelihood and amount of
data being analysed. There have been various efforts to reduce this cost by
means such as re-parameterisations [14], reduced order methods [15, 16],
parallel sampling methods [17, 18, 19] and problem-specific sampling
algorithms [11, 20, 21, 22].

Machine learning algorithms have successfully been applied to various aspects
of gravitational-wave data analysis [23] including data-quality improvement,
waveform modelling, searches and, most relevant to this work, parameter
estimation [24, 25, 26, 27]. These algorithms, once trained, greatly reduce
the cost of computing posterior distributions. However they are still in their
infancy and there are various challenges that have yet to be solved, for
example, they are constrained by the distribution of training data and cannot
reliably be applied to data that is outside that distribution, e.g. with
different detector noise curves or longer duration signals.

Normalising flows are a type of generative machine learning algorithm that
have recently been applied to a wide variety of problems including generation,
inference and representation learning [28, 29]. The fundamental idea behind
normalising flows is to model a complex probability distribution as the
transformation of simple distribution. This transformation is constructed to
be highly flexible but maintain an explicit mathematical description with a
tractable Jacobian. This makes them well suited to applications in the
physical sciences where the apparent opacity of other machine learning
algorithms can hinder their widespread adoption.

In this paper we propose a modified nested sampling algorithm that
incorporates a novel proposal method using normalising flows. During sampling
the normalising flow is trained such that the resulting distribution can be
re-sampled according to the prior to produce new live points within a given
iso-likelihood contour which are independent of the current worst point. This
is akin to other nested sampling algorithms which sample from the constrained
prior thus avoiding the need to evolve points with a random walk. These
algorithms typically require fewer evaluations of the likelihood compared to
those that use random walks which reduces the computational cost. However, the
challenge when implementing these algorithms is determining the constrained
prior which often requires defining multiple bounding distributions. Using a
normalising flow allows us to use a single bounding distribution but comes at
the cost of training during sampling. This additional cost is easily offset in
problems where there is significant computational cost associated with
computing the likelihood. As such we choose to apply our algorithm to
gravitational-wave inference.

This paper is structured as follows: in section II we outline the background
theory for gravitational-wave inference, nested sampling and normalising
flows. We then introduce our method in section III and discuss related work in
section IV. In section V we present results for gravitational-wave parameter
estimation, compare to a commonly used nested sampler, and discuss the
implication of our results. Finally in section VI we summarise our results and
draw conclusions.

##  II Background

###  II.1 Gravitational-wave likelihood for compact binary coalescence

The strain data from gravitational-wave interferometers is modelled as a
signal with additive noise in the time domain [11]. The standard definition of
the likelihood for compact binary coalescence then assumes that the noise is
Gaussian and stationary over short observation times and that the power
spectral density $S_{n}$ is independent of the model parameters. This allows
us to define the likelihood as an independent product in the frequency domain

$p(\bm{d}|H,\bm{\theta})\sim\exp\sum_{i}\left[\frac{-2|\tilde{h}_{i}(\bm{\theta})-\tilde{d}_{i}|^{2}}{TS_{n}(f_{i})}\right],$
(2)

where $\tilde{h}_{i}(\bm{\theta})$ and $\tilde{d}_{i}$ are the waveform model
and data in the frequency domain. There are various different waveform models
$\bm{h}(\bm{\theta})$ which are characterised by sets of parameters. The exact
number of parameters depends on the physics that the model describes, if both
compact objects are assumed to be black holes, then typically 15 parameters
are required [11]. Different parameterisations have proven to improve sampling
efficiency and convergence, most notably reparameterising the component masses
in terms for chirp mass $\mathcal{M}=(m_{1}m_{2})^{3/5}(m_{1}+m_{2})^{-1/5}$
and asymmetric mass ratio $q=m_{2}/m_{1}$ [11] and using the system-frame when
describing the orientation of the binary [14].

###  II.2 Nested sampling

Nested sampling is a stochastic algorithm for Bayesian inference that was
proposed by Skilling in [13] and computes the Bayesian evidence
$Z=p(\bm{d}|H)=\int p(\bm{d}|\bm{\theta},H)\textrm{d}\bm{\theta}$ in eq. 1. In
nested sampling the evidence integral is simplified by considering the total
prior volume $X$ contained within a given iso-likelihood contour
$L^{*}=p(\bm{d}|\bm{\theta},H)$:

$X(L^{*})=\int_{p(\bm{d}|\bm{\theta},H)>L^{*}}\textrm{d}\bm{\theta}p(\bm{\theta}|H).$
(3)

where $\textrm{d}\bm{\theta}p(\bm{\theta}|H)=dX$. The evidence integral can
then be re-written using the likelihood $L$ by inverting eq. 3

$Z=\int_{0}^{1}L(X)\textrm{d}X.$ (4)

Since the integrand is positive and decreasing, the function is well behaved
and the integral can therefore be approximated by considering an ordered
sequence of decreasing of $M$ points in the prior volume $X_{i}$, evaluating
the likelihood at each point $L_{i}=L(X_{i})$ and finally, for example, using
the trapezoid rule

$Z=\sum_{i}^{M}\frac{1}{2}(X_{i-1}-X_{i+1}L_{i}).$ (5)

The complete algorithm is detailed in [13] but in short it requires first
sampling a set of points, known as live points, from the prior distribution.
The point with the lowest likelihood $L^{*}$ is then removed and replaced by
another sample drawn from within the likelihood contour defined by the
original point and according to the prior distribution. This replacement
process is then continued until a stopping criteria is met [13]. In
gravitational-wave inference this is typically chosen to be when the
fractional change in log evidence at a given iteration $\Delta\ln Z_{i}$
caused by replacing another point is less than $0.1$ [11], where $\Delta\ln
Z_{i}\equiv\ln\left(Z_{i}+\Delta Z_{i}\right)-\ln Z_{i}$ and $\Delta
Z_{i}\approx L^{*}X_{i}$. Once the algorithm has terminated the final evidence
is computed including the final set of live points and posterior samples can
be produced from the set of discarded and final live points [13].

The challenge when implementing this algorithm is efficiently proposing the
replacement samples within the current iso-likelihood contour since these
samples must be independently and identically distributed (i.i.d.) according
to the prior. One commonly used method draws new points based on the location
of existing points which are then evolved using a random walk to ensure they
are i.i.d. [11, 18]. The efficiency of this method depends on the length of
the MCMC chain required since each step in the chain requires evaluating the
likelihood, for GW inference proposing a single new sample often requires of
order $10^{3}$ likelihood evaluations. Other approaches avoid using a random
walk by sampling from the constrained prior at a given iteration. This
requires determining the distribution of the current live points and
constructing one or more bounding distributions, for example, using
multidimensional ellipsoids as described in [30]. In this case the challenge
is efficiently constructing the bounding distributions and sampling them
without over-constraining the prior and consequently under-sampling a region
in parameter space. What is more, these two approaches are sometimes used in
conjunction [31] to tackle particularly challenging problems, such as
gravitational-wave inference [32], where the individual methods can be
inefficient. In section III we propose a method for directly sampling from a
constrained prior using normalising flows to learn the iso-likelihood
contours.

###  II.3 Normalising flows

Normalising flows [33] are a type of generative machine learning algorithm
that map samples $x$ in the physical space111We use the term physical space to
distinguish it from the strain data however in the literature $\mathcal{X}$
often referred to as the data space. $\mathcal{X}$ to samples $z$ in a latent
space $\mathcal{Z}$ such that the samples $z$ are distributed according to a
prior probability distribution $p_{\mathcal{Z}}$, known as the latent prior.
If we assume this mapping to be bijective and denoted
$f:\mathcal{Z}\to\mathcal{X}$, then using the change of variable formula we
can define a probability distribution in the physical space $p_{\mathcal{X}}$

$p_{\mathcal{X}}(x)=p_{\mathcal{Z}}(f(x))\left|\textrm{det}\left(\frac{\partial
f(x)}{\partial x^{T}}\right)\right|,$ (6)

where $\partial f(x)/\partial x^{T}$ is the Jacobian of $f$ at $x$. Samples
can then be drawn from $p_{\mathcal{X}}$ by sampling the latent space $z\sim
p_{\mathcal{Z}}$ and applying the inverse mapping $f^{-1}$. This requires the
bijective mapping $f$ to have a tractable Jacobian, which for functions with
multidimensional domains and codomains is often computationally expensive.
These mappings are therefore carefully constructed to meet these criteria and
usually parameterised by a neural network which allows for the mapping to be
learned.

Normalising flows generally fall in two categories depending on how the
mappings are defined: autoregressive flows [34, 35] and flows based on
coupling transforms [36, 37, 38]. Each have distinct advantages and
disadvantages, most notably autoregressive flows are often more flexible but
have a greater computational cost and are not always analytically invertible
whereas coupling based flows tend to be less flexible but computationally
cheaper to evaluate and analytically invertible. In this work we use the
latter since we aim to minimise computational cost and later rely on the
invertablility of normalising flows.

Figure 1: Diagram of a normalising flow $f(x)$ composed of four coupling
transforms which maps an $n$-dimensional input vector x to an $n$-dimensional
latent vector z. Each transform splits $x$ in two $[x_{1:m},x_{m+1:n}]$ and
updates one part conditioned on the other. In the first and third transforms
$x_{1:m}$ is used as the input to a neural network (NN) which then produces
the scale $s$ and translation $t$ vectors of length $m$. The element-wise
product ($\odot$) is then computed between $x_{1:m}$ and $\exp(s)$ followed by
the sum of the output and $t$. This is shown in the left transform. In the
second and fourth transforms $x_{1:m}$ is updated conditioned on $x_{m+1:n}$
as shown in the right transform.

Flows based on coupling transforms map an $n$-dimensional input vector $x$ to
an $n$-dimensional output vector $y$ by splitting the input vector in two
parts $x=[x_{1:m},x_{m+1:n}]$ where $m<n$ and then transforming one part
conditioned on the other unchanged part. The choice of splitting depends on
the specific coupling transform though typically an alternating pattern is
used. For example the coupling transform proposed in [36]:

$\displaystyle y_{1:m}$ $\displaystyle=x_{1:m},$ (7a) $\displaystyle
y_{m+1:n}$
$\displaystyle=x_{m+1:n}\odot\exp\left[s(x_{1:m})\right]+t(x_{1:m}),$ (7b)

where $s(x_{1:m})$ and $t(x_{1:m})$ are the output of a neural network and
$\odot$ denotes the element-wise product. We show a schematic of eq. 7 in fig.
1. Since $y_{1:m}=x_{1:m}$ the Jacobian matrix of these transforms is lower
triangular and the Jacobian determinant is simply the product of the diagonal
entries. This transform can easily be inverted by inverting eq. 7 and using
$y_{1:m}$ as the input to $s$ and $t$:

$\displaystyle x_{1:m}$ $\displaystyle=y_{1:n},$ (8a) $\displaystyle
x_{m+1:n}$
$\displaystyle=\left[y_{m+1:n}-t(y_{m:1})\right]\odot\exp\left[-s(y_{1:m})\right].$
(8b)

As with eq. 7 the Jacobian determinant is trivial to compute. Coupling
transforms can also be stacked by alternating which part of $x$ is updated,
$x_{1:m}$ or $x_{m+1:n}$, to produce more flexible transforms. It is also
common practice to include permutations between coupling transforms, these
permute the inputs to the transforms making the mapping more expressive [36].
More recently this has been generalised to linear transforms in which the
permutation is learnt during training [37]. In fig. 1 we show a schematic of
how coupling transforms can be stacked with permutations. For a more detailed
description of normalising flows, the different types and their application we
point the reader to [28, 29].

##  III Method

We present a novel method for sampling within a given iso-likelihood using a
normalising flow. The normalising flow learns the distribution of a set of
live points and is constructed such that the learnt distribution can be
sampled from analytically. We introduce additional steps to ensure that the
samples are distributed according to the sampling prior222We use sampling
prior to denote the prior used for nested sampling and to distinguish it from
the latent prior use in normalising flows. and bounded by the iso-likelihood
contour. This eliminates the need to evolve new samples and allows us to
efficiently draw new samples within a complex iso-likelihood contour. A more
efficient proposal equates to few rejected points and, since the likelihood
must be computed before a point can be rejected, this is also equivalent to a
reduction in the number of likelihood evaluations. We first describe our
method in the isolated case of a single of set of live points at a given
iteration and then present a nested sampling algorithm which incorporates it.

Figure 2: Example of how a normalising flow trained on a set of live points
can produce samples within current iso-likelihood contour for simple two-
dimensional parameter space. Top: example of training samples in the physical
space $\mathcal{X}$ and learned mapping to the latent space $\mathcal{Z}$ with
the iso-likelihood contour for the current worst point shown in orange.
Middle: samples drawn from a truncated Guassian within the iso-likelihood
contour in $\mathcal{Z}$ and mapped to $\mathcal{X}$ using the inverse
mapping. Bottom: pool of accepted samples after applying rejection sampling
until 1000 points are obtained shown in both $\mathcal{Z}$ and $\mathcal{X}$.

###  III.1 Sampling within a iso-likelihood contour

At any given point in the nested sampling algorithm there is a current set of
live points which by definition is contained within the iso-likelihood contour
defined by the worst point with likelihood $L^{*}$. We describe the
implementation in our algorithm in terms of four steps:

#### How to define an iso-likelihood contour using a normalising flow.

If we treat the sampling space as the physical space $\mathcal{X}$ and the
live points as the data $x$ we can then train a normalising flow to
approximate the distribution of the live points to within some error and use
it to draw new samples in $\mathcal{X}$. This requires sampling from the
latent prior $p_{z}$, which we choose to be an $n$-dimensional Gaussian, and
then applying the inverse mapping learned by flow $f^{-1}$. However, since our
choice of latent prior $p_{\mathcal{Z}}$ has an infinite domain the resulting
distribution of samples $p_{\mathcal{X}}$ will not be bounded by the iso-
likelihood contour or distributed according to the sampling prior. We
therefore examine the notion of an iso-likelihood contour in the context of
the normalising flow with an $n$-dimensional Gaussian latent prior.

#### How to determine the contour given the current set of live points.

Once trained, the normalising flow can be used to map the current worst point
in the physical space $x^{*}$ to the latent space $\mathcal{Z}$ as $z^{*}$.
This point has a likelihood in the latent space $L_{\mathcal{Z}}^{*}$ given by
$p_{\mathcal{Z}}(z^{*})$ and, since $p_{\mathcal{Z}}$ is an $n$-dimensional
Gaussian, points of equal likelihood lie on the $(n-1)$-sphere with radius
$r^{*}$ given by $z^{*}$. If we assume a perfect mapping, then this iso-
likelihood contour in the latent space can be mapped to an iso-likelihood
contour in the physical space. We can therefore sample within the contour in
latent space and use the inverse mapping to produce samples within the contour
in the data space.

#### How to sample within the contour.

We use two approaches for drawing $K$ new samples $\bm{z}_{i}$ in the latent
space given a radius $r^{*}$, these produce normally and uniformly distributed
samples respectively. Both start by drawing $K$ samples on the $(n-1)$-sphere
using the algorithm proposed in [39, 40] where $K$-dimensional vectors are
drawn from an $n$-dimensional unit Gaussian and then normalised using the
Euclidean norm. These samples $\bm{y}_{i}$ can then be rescaled to obtain
samples within in the $n$-ball

$\bm{z}_{i}=\rho_{i}\frac{\bm{y}_{i}}{||\bm{y}_{i}||_{2}},$ (9)

where the choice of distribution for $\rho_{i}$ determines how the resulting
samples $\bm{z}_{i}$ are distributed in $\mathcal{Z}$. For uniformly
distributed samples $\rho_{i}=u^{1/n}$ where $u\sim\mathcal{U}(0,r^{*})$ and
for normally distributed samples $\rho_{i}\sim\chi(n)$ where $\chi(n)$ is a
chi-distribution with $n$ degrees of freedom truncated at $r$. The inverse
mapping of the normalising flow can then be applied to $\bm{z}_{i}$ to obtain
samples in the physical space. We consider two approaches because sampling
from a truncated Gaussian in high-dimensional space can become inefficient for
large values of $r^{*}$.

#### How to ensure new samples are drawn according to the prior.

The samples obtained in the previous step must be re-sampled such that they
are distributed according to the sampling prior. We use rejection sampling and
compute weights $\alpha_{i}$ for each sample

$\alpha_{i}=\frac{p(\bm{x}_{i})}{q(\bm{x}_{i})},$ (10)

where $x_{i}=f^{-1}(z)$, $p(x)$ is the sampling prior and $q(x)$ the proposal
probability which is computed using the inverse of eq. 6

$q(x)=q(f^{-1}(z))=p_{\mathcal{Z}}(z)\left|\textrm{det}\left(\frac{\partial
f^{-1}(z)}{\partial x^{T}}\right)\right|^{-1}.$ (11)

The choice of latent prior $p_{\mathcal{Z}}$ will depend on which method was
used to draw the samples in the latent space. The weights eq. 10 are then
rescaled such that their maximum value is one. We then draw $N$ samples
$u\sim\mathcal{U}[0,1]$ and accept samples for which $\alpha_{i}/u_{i}>1$. In
fig. 2 we show an example of this process for a simple two-dimensional case.

###  III.2 Algorithm details

We identify three key stages of the sampling approach detailed in section
III.1 that we then incorporate into the nested sampling algorithm:

  * •

Training: a normalising flow is trained on the current set of $K$ live points
by minimising a Monte Carlo approximation of the of the KL divergence between
the target distribution and $p_{\mathcal{X}}(x)$ as defined in eq. 6. The
explicit loss function is derived in appendix A.

  * •

Population: once the normalising flow is trained, samples are drawn within the
$n$-ball of radius $r$ defined by the worst point. These samples are then
mapped to the data space $\mathcal{X}$, re-sampled according to the prior and
finally stored in the pool of new samples.

  * •

Proposal: once the pool of new samples has been populated, new live points are
drawn at random from the pool and then removed until the pool is empty or the
normalising flow is retrained.

The standard nested sampling algorithm is modified to include these stages, we
call the algorithm Nessai.

The start of the algorithm remains unchanged: K live points are drawn from the
prior distribution and their log-likelihoods computed. We then start the
iterative process of determining the worst live point with log-likelihood
$L^{*}$ and drawing a replacement live point that lies within the iso-
likelihood contour. In our modified algorithm we use standard rejection
sampling from the prior for the first $M$ points (typically $2K$) or until it
becomes inefficient. The normalising flow is then trained on the current $K$
live points allowing us to map the worst point to the latent space to obtain
the worst latent point $z_{w}$ and the radius of the corresponding $n$-ball
$r$. This radius is then used for the population stage when drawing samples in
the latent space. Once populated, a replacement point is drawn from the pool,
its log-likelihood is computed and if it is greater than $L^{*}$, the point is
accepted; if not, more points are drawn until one is accepted. The proposal
stage is then repeated for subsequent worst live points until one of four
criteria is met:

  * •

the proposal pool is depleted: the normalising flow is retrained using the
current live points and worst live point is used to compute a new radius and
the population stage is repeated,

  * •

the acceptance rate falls below a user-defined criteria: the current proposal
pool is discarded and the normalising flow is retrained. This threshold is
defined by the user,

  * •

the criteria for retraining the normalising flow is met: the normalising flow
is retrained with the current $K$ live points, this happens by default every
$K$ live points,

  * •

the nested sampling convergence criterion is met: the algorithm terminates.

Since this algorithm relies on the normalising flows’ ability to approximate
the distribution of live points at various stages throughout the sampling, we
include a series of reparameterisations of $\mathcal{X}$ to reduce the
complexity of the data space and removing certain features. We denote this
reparameterised space $\mathcal{X}^{\prime}$ and include the Jacobian for each
reparameterisation in eq. 11. These reparamerisations are:

  * •

Rescaling: we add the option to rescale the input data according to either the
sampling priors or the current minimum and maximum values such that all of the
parameters in $\mathcal{X}^{\prime}$ are defined over the same domain. As as
default we use $[-1,1]^{n}$.

  * •

Boundary inversion: we observe that asymmetric distributions with high density
regions near the prior bounds are often under-sampled. To mitigate this effect
we add the option to mirror the live points around such bounds and train on
the resulting symmetric distribution. Further details are provided in appendix
B.

We also introduce additional settings which help with convergence and sampling
efficiency, some of these and discussed in section V.5 and a comprehensive
list can be found in the online documentation for our sampler [41].

###  III.3 Gravitational-wave reparameterisations

The gravitational-wave parameter space is typically $15$-dimensional and
contains various degeneracies between parameters such as the masses,
inclination and luminosity distance which can make sampling inefficient.
Previous work has shown that certain reparameterisations can improve sampling
efficiency [11]. We use two of these: chirp mass $\mathcal{M}$ and asymmetric
mass ratio $q$ replace the component masses and we use the system-frame
parameterisation in place of the radiation-frame to describe the orientation
of the binary [14].

More than half of the parameters to sample are angles and we note that the
periodicity of these angles is not encoded in the mapping learned by the
normalising flow since the latent space $\mathcal{Z}$ is continuous and
unbounded. We therefore include a further reparameterisation specifically for
the angular parameters $\theta_{i}$. We assume that each angle has a
corresponding radial component $\rho_{\theta_{i}}$ and together they describe
a position in a two-dimensional plane. We can therefore use standard
transformations to express this position in Cartesian coordinates
$(x_{\theta_{i}},y_{\theta_{i}})$:

$\begin{split}x_{\theta_{i}}=\rho_{\theta_{i}}\cos{\theta_{i}},\\\
y_{\theta_{i}}=\rho_{\theta_{i}}\sin{\theta_{i}}.\end{split}$ (12)

If we choose the distribution of radial components such that
$\rho_{\theta_{i}}\in[0,\infty)$ then
$x_{\theta_{i}},y_{\theta_{i}}\in(-\infty,\infty)$. Since we are using a
Gaussian latent prior $p_{\mathcal{Z}}$, we sample $\rho_{\theta_{i}}$ from a
chi-distribution with two degrees of freedom such that, if the angle is
uniformly distributed on $[0,2\pi]$, the resulting distribution of
$(x_{\theta_{i}},y_{\theta_{i}})$ is Gaussian. We use this treatment for the
phase, inclination, polarisation and all four spin angles, for polarisation we
rescale the angles to $[0,2\pi]$ before applying the transformation to
Cartesian coordinates. This reparameterisation also naturally includes
periodic boundary conditions for the angles with uniform priors.

The sky location is described by a further two angles, right ascension
$\alpha$ and declination $\delta$. For these angles we extend the previous
treatment from two-dimensional to three-dimensional Cartesian coordinates
$(x,y,z)$ and draw the radial component $\rho$ from a chi-distribution with
three degrees of freedom. For the standard priors,
$p(\alpha)\sim\mathcal{U}[0,2\pi]$ and $p(\delta)\sim\cos{\delta}$, the
resulting distribution of $(x,y,z)$ is again Gaussian.

The spin magnitudes $\chi_{1}$ and $\chi_{2}$ also require a specific
treatment. They are typically defined on $[0,0.99]$ with uniform priors and,
importantly, the posterior distributions are often broad and span the entire
prior range. We consider applying the boundary inversion to both bounds but in
practice find this ineffective. We instead opt to map $\chi_{i}$ into a two-
dimensional plane with positions described using Cartesian coordinates
$x_{\chi_{i}}$ and $y_{\chi_{i}}$. We achieve this by first defining a
rescaled magnitude $\hat{\chi}_{i}\in[0,1]$ which is obtained using the
corresponding priors. Then, we consider the angle defined by
$\hat{\chi_{i}}\pi$ and, again, introduce a radial component
$\rho_{\chi_{i}}\sim\chi(2)$. The corresponding Cartesian coordinates
$(x_{\chi_{i}},y_{\chi_{i}})$ are defined on $[0,\infty)$ and
$(-\infty,\infty)$ respectively. However we know that the coupling transforms
we have chosen to use are better suited to unbounded domains. To avoid this,
we introduce a random variable $k$ which is drawn from a Rademacher
distribution and include it in the Cartesian coordinate transform

$\begin{split}x_{\chi_{i}}&=\rho_{\chi_{i}}\cos{\hat{\chi_{i}}\pi},\\\
y_{\chi_{i}}&=k\rho_{\chi_{i}}\sin{\hat{\chi_{i}}\pi}.\end{split}$ (13)

As a result of the including $k$,
$(x_{\chi_{i}},y_{\chi_{i}})\in(-\infty,\infty)$ and the hard boundary at
$x_{\chi_{i}}=0$ has been avoided.

We choose to reparameterise the luminosity distance $d_{\text{L}}$ such that
the prior for the resulting parameter $d_{\text{U}}$ is uniform. The exact
reparameterisation therefore depends on the prior used for $d_{\text{L}}$. In
this work we choose to use a prior on $d_{\text{L}}$ that is uniform in co-
moving volume so we first convert the luminosity distance to a co-moving
distance $d_{\text{C}}$ and then the uniform parameter is simply
$d_{\text{U}}=d_{\text{C}}^{3}$. Similar reparametersitions can be determined
for other commonly used distances priors such as a power law. Additionally, we
allow boundary inversion as described in section III.2 but limit it to only
the upper bound since in practice the luminosity distance posterior will not
rail against the lower bound.

The remaining parameters are the chirp mass $\mathcal{M}$ and mass ratio $q$
for which we use the reparameterisations from $\mathcal{X}$ to
$\mathcal{X}^{\prime}$ mentioned in section III.2, allowing boundary inversion
for $q$.

###  III.4 Implementation

We use the implementation of normalising flows in PyTorch [42] available in
nflows [43] which allows for a wide variety of normalising flows to be used.
However we choose to use coupling transforms [36] because of the tractable
Jacobian and ease of computing the inverse mapping. As suggested in [44, 37],
we include invertible linear transforms that randomly permute the parameters
before each coupling transform allowing all of the parameters to interact with
each other. We also include batch normalisation [45] after each coupling
transform as described in [36]. We use a residual neural network [46, 47] for
computing the parameters for each transform. We train the normalising flows
with the Adam optimiser [48]. In appendix D we detail the specific parameters
used for the results presented in section V.

Our sampler, Nessai (Nested Sampling with Artificial Intelligence), is
available as an open source package [49] and documentation is also available
online [41].

##  IV Related Work

Different frameworks and samplers have been developed for gravitational-wave
inference. LALInference [11] implements nested sampling and MCMC with specific
proposal methods for the gravitational-wave parameter space and has been used
extensively for analyses of the first gravitational wave detections [50] and
GWTC-1 [6]. More recently, the Python package Bilby [22] has been developed to
use off-the-shelf samplers, such as dynesty [31], and been shown to achieve
comparable results to LALInference on GWTC-1 [32].

Machine learning has previously been incorporated into stochastic sampling
algorithms; in [51] the likelihood function is approximated with a neural
network, and in [52] neural networks are used to generalise Hamiltonian Monte
Carlo. More closely related to our work, normalising flows have been used to
improve the efficiency of MCMC methods by reparameterising the sampling space
[53] and a similar approach has also been extended to MCMC sampling in nested
sampling in [54].

Recent work has shown that likelihood-free inference using conditional
variational autoencoders [24, 25] and normalising flows [26, 27] can produce
posterior distributions for compact binary coalescence from binary black
holes. These approaches promise to drastically reduce the cost of producing
posterior samples when compared to traditional stochastic sampling methods.
However, they require large amounts of training data and they currently lack
the flexibility to deal with, for example, different PSDs, high sampling
frequencies and long duration signals.

(a) With phase marginalisation (b) With phase and distance marginalisation

Figure 3: Probability-probability (P-P) plot showing the confidence interval
versus the fraction of the events within that confidence interval for the
posterior distributions obtained using our analysis Nessai for 128 simulated
compact binary coalescence signals produced with Bilby and bilby_pipe. The 1-,
2- and 3-$\sigma$ confidence intervals are indicated by the shaded regions and
$p$-values are shown for each of the parameters and the combined $p$-value is
also shown.

##  V Results

We chose to evaluate the performance of our sampling algorithm, Nessai, with
simulated gravitational-wave signals from CBC. The parameter is multi-
dimensional and various parameters are correlated which can prove challenging
when sampling. The likelihood, as defined in section II.1, is also typically
computationally costly to evaluate, making it well suited to our sampler,
however, this depends on the waveform approximant used and length of the
observation. We first use probability-probability plots to check the
consistency of our sampler and then compare our results to those obtained
using dynesty. We then highlight how the likelihood computation can be
parallelised in Nessai before finally discussing various diagnostics that can
ne used to identify problems during sampling and tune the sampler settings.

We use bilby_pipe and Bilby [22] to simulate 128 injections using IMRPhenomPv2
[55, 56] sampled at 2048 Hz with 4-second observing time in a three detector
network with AdLIGO Hanford, AdLIGO Livingston and AdVirgo at design
sensitivity [4]. We choose uniform priors on chirp mass
$\mathcal{M}\sim\mathcal{U}[25,35\,\mathrm{]}{\text{M}_{\odot}}$ and
asymmetric mass ratio $q\sim\mathcal{U}[0.125,1.0]$, a prior on luminosity
distance that is uniform in co-moving volume on
$[100,2000\,\mathrm{]}\,{\text{Mpc}}$, a uniform prior for the reference time
at the geocentre with width 0.2 and the remaining priors are set to the
defaults for precessing binary black holes in Bilby [32], see appendix C for a
complete list.

We analyse the injections with our sampling algorithm, Nessai, outlined in
section III.2 and include the specific reparameterisations for gravitational
wave analyses described in section III.3. We choose to analyse each injection
twice: once with just phase marginalisation and once with both phase and
distance marginalisation. Further details of the exact settings used for
Nessai are provided in appendix D.

###  V.1 Result validation

Probability-probability (P-P) plots are a standard method of verifying the
performance of sampling algorithms [57, 58]. They test whether the correct
proportion of injected values are recovered at a given confidence interval for
a specific prior distribution. These tests are particularly useful when using
a Gaussian likelihood, such as eq. 2, since the fraction of events within a
given confidence interval should be uniformly distributed and we can therefore
compute $p$-values for each parameter and a combined $p$-value for all of the
parameters. We produce P-P plots for both of our analyses using Bilby and
present the results in fig. 3. For an idealised sampler for the $p\%$
confidence interval, $p\%$ of the events should be recovered, this would
correspond to a diagonal line. In practice we expect to see deviation from the
diagonal, as such the 1-, 2- and 3-$\sigma$ confidence intervals are also
shown in fig. 3. These results show that Nessai consistently recovers for the
posteriors for the 128 injections but also indicate that luminosity distance
is consistently harder to sample. The combined $p$-values of 0.3394 and 0.6818
for our analyses without and with distance marginalisation serve as further
verification.

###  V.2 Comparison to dynesty

To further validate our results we compare them to those obtained with dynesty
[31], another nested sampling algorithm commonly used in gravitational-wave
inference [22, 32]. We use the configuration described in [32] but increase
the number of live points to 2000 and run on a single thread to ensure as
direct of a comparison with Nessai as possible. With these settings dynesty
passes the P-P test (see appendix E) but we note that these settings are the
minimum required to produce reliable results and in practice more conservative
settings are often used. Additionally, several injections required a second
analysis with a different sampling-seed in order reach convergence. The
results obtained with dynesty allow us to verify the log-evidences returned by
Nessai since these cannot be computed analytically and provides a point of
reference when considering the number of likelihood evaluations and total
computational time.

Figure 4: Difference between the log evidences $\Delta\ln Z$ obtained using
dynesty and Nessai for all 128 injections with distance marginalisation
(dashed line) and without distance marginalisation (solid line).

In fig. 4 we compare the log-evidences returned by dynesty and Nessai. If
Nessai was consistently over or under-estimating the log-evidence when
compared to dynesty, this would indicate a potential problem during sampling,
such as over- or under-constraining, which would lead to biased results. The
results in fig. 4 show no such bias. However, since sampling is a stochastic
process there is an error associated with the computed log-evidence. The
theoretical error can be approximated using the information $H$ and the number
of live points $K$, $\delta\log Z\approx\sqrt{H/K}$. To quantify this error we
repeat the analysis on a single injection with 50 different sampling seeds and
compute an approximate error $\delta\log Z\approx 0.092$. In practice we
observe a wider spread of log-evidences of $0.11$, this is consistent with
previous analyses which determined that there are additional sources of
uncertainty [59].

Figure 5: Distribution of the total number of likelihood evaluations required
to reach convergence for and total run time for dynesty (blue) and Nessai
(orange) when applied 128 simulated signals from compact binary coalescence
with the priors and sampler settings described in sections V and D. The
results distance marginalisation disable are shown with solid lines and those
with distance marginalisation enabled are shown with dashed lines.

Nessai is designed with the aim of improving the efficiency of drawing
replacement live points at the cost of repeatedly training a normalising flow
and populating a pool of live points. An improvement in the efficiency
translates to a reduction in the total number of likelihood evaluations since
the likelihood must be computed for each rejected point. We therefore compare
the total number of likelihood evaluations required to reach convergence for
each sampler in fig. 5 with and without distance marginalisation. Nessai
requires a median of $5.04\times 10^{6}$ and $7.22\times 10^{6}$ likelihood
evaluations to converge with and without distance marginalisation respectively
and dynesty requires $10.44\times 10^{6}$ and $9.67\times 10^{6}$. In contrast
to dynesty, sampling with Nessai is more efficient without distance
marginalisation, we attribute this to a combination of the reparameterisation
used for luminosity distance and the sampler settings we converged on for
Nessai.

This, however, does not directly translate into a decrease in the total run-
time since it does factor in the cost of training and population. In fig. 5 we
also show the total run-time for each sampler and when comparing the median
run-times we observe that Nessai is 2.07 times faster than dynesty without
distance marginalisation and 1.34 times faster with it. We also examine the
proportion of the run-time spent on training and population and find that on
average population-time accounts for approximately 39.6% of the total run-time
and training-time accounts for a further 8.1%. This explains why, despite a
2.32 times reduction in likelihood evaluations, we do not observe a 2.32 times
reduction in sampling time compared to dynesty. We also note that the cost of
training and population does not depend on the cost of evaluating the
likelihood, as such, the fraction of the total run-time will decrease as the
cost of evaluating the likelihood cost increases.

Figure 6: Comparison of the total time (in hours) spent on each stage of the
algorithm for increasing number of threads for a single injection with a fixed
noise seed. The time spent evaluating the likelihood decreases as the number
of threads increases, the theoretical reduction is shown in black. The
training and population stages remain approximately constant, as such act a
lower bound on the minimum run-time. The sum of the time spent on likelihood
evaluation, training and population is approximately equal to the total time
spent sampling, indicating minimal overhead.

For each injection we can also compare the posterior distributions produced by
each sampler. These allow us to quickly identify discrepancies between
samplers for specific injections or regions of the parameter space. However
these differences are not easily quantified and the correlations between more
than two parameters are not clearly represented. We show an example of such a
comparison is appendix F.

###  V.3 Parallelisation of the likelihood computation

Our sampler is designed such that candidate live points are drawn
simultaneously in the population stage. This allows for simple parallelisation
of the likelihood computation since the pool of candidate live points can be
distributed over a number of threads and likelihood values computed and stored
until needed for the proposal stage. In fig. 6 we compare the run-time and
time spent evaluating the likelihood for the same injection using increasing
number of threads for the likelihood computation. We use an additional thread
for the main sampling process. The time spent evaluating the likelihood is
inversely proportional to the number of threads allocated although the overall
run-time is not. With a single thread it accounts for 53.9% of the total run-
time and this decreases to 9.0% when using 16 threads. As mentioned
previously, there is a cost associated with populating stage and further
smaller cost associated with training, for this injection these are 36.4% and
9.3% respectively. These remain approximately constant when increasing the
number of threads available and act as a lower limit on the theoretical
minimum run-time, this is shown in fig. 6. There remaining $<1\%$ of the run-
time is general overhead associated with running the sampler.

###  V.4 Diagnostics

As mentioned previously, there are various challenges when implementing a
sampling algorithm. Nessai is designed to sample from within the constrained
prior, in this case care must be taken to ensure that the prior is not over-
constrained since this will lead to regions of parameter being under-sampled
which in turn will bias the results. There are also specific problems that
arise from the nature of the parameter space, such as multi-modality and
correlations. We use a series of diagnostics to identify possible problems
during sampling. In section V.5 we also discuss how some of these diagnostics
can be used to tune the sampler settings described in sections III.2 and D.

Figure 7: Example of the distribution of insertion indices for two nested
sampling runs with 2000 live points. Uniformly distributed indices indicate no
under- or over-constraining and deviations from uniformity indicate the
opposite. The result with an orange dashed line shows over-constraining and
the result with a solid blue line shows the correctly converged run. The
shaded region indicates the 2-$\sigma$ errors on the expected distribution.

We use the cross-checks proposed in [60] as a heuristic for determining if the
nested sampling algorithm has converged without over or under-constraining the
posterior distributions. These checks rely on order statistics and the
assumption that new live points should be inserted uniformly into the existing
live points which allows for a $p$-value to be computed using the Kolomorgov-
Smirnov statistic [61] with the additional consideration that underlying
distribution is discrete [62]. In fig. 7 we show an example of the
distribution of the indices of newly inserted live points and in fig. 8 we
show the $p$-values computed every $K$ iterations. The histogram shows the
final distribution of insertion indices for all the nested samples, this may
not highlight specific problematic regions of the parameter space but if it is
not uniform, it is a clear indication that the sampler is consistently over-
or under-constraining. If the distribution of $p$-values in fig. 8 is non-
uniform then this is another clear indication of problems during sampling.

Figure 8: Example of the statistics that are tracked in our sampler as a
function of sampling iteration: (a) minimum (blue solid) and maximum (red
dashed) log-likelihood, (b) cumulative number of likelihood evaluations, (c)
log evidence $\log Z$ (blue solid) and fractional change in evidence
$\text{d}Z$ (red dashed), (d) proposal (blue solid) and population (red
dashed) acceptance and (e) $p$-value for cross-checks of every $K$ live
points. The iterations at which the normalising flow is trained are indicated
with vertical lines.

The acceptance is another important statistic to monitor during sampling since
we aim to develop a more efficient sampler. There are two acceptances we can
monitor in Nessai, the proposal acceptance and the population acceptance. The
first has a direct effect on the number of likelihood evaluations whilst the
second only affects the total run-time, both quantities are shown in fig. 8.
This figure also highlights how periodically retraining the normalising flow
leads to an increase in the proposal efficiency. It also shows how the
population process is typically inefficient which explains why on average 32
percent of the total run-time is spent on the population stage.

We also track the minimum and maximum log-likelihoods, number of log-
likelihood evaluations, log-evidence and fractional change in evidence. The
combination of these statistics allows the user to quickly understand the
current state of the sampler and identify potential issues such as plateaus in
the likelihood space and regions which are inefficient to sample. The complete
set of statistics is shown in fig. 8.

###  V.5 Tuning Nessai

Nessai includes various settings, a comprehensive list and description of each
can be found in the documentation [41]. In practice we find that a small
subset of the settings predominantly determine whether the algorithm converges
without any bias. We use the validation method described in section V.1 and
the diagnostics from section V.4 to understand how these settings affect
convergence.

As expected, the number of live points $K$ is an important setting but it is
even more crucial in Nessai since it limits the amount of training data
available. We find that a minimum of 1000 live points is required and for more
complex problems, such as gravitational-wave inference, at least 2000 live
points should be used.

There are large number of settings which relate to the complexity of the
normalising flow. Whilst tuning the sampler we found that the number of
coupling transformations greatly affected convergence. If too many transforms
were used the algorithm was prone to over-constraining the posterior
distribution. We attribute this to the complexity of the iso-likelihood
contour learnt by the flow, if the flow has too many trainable parameters it
can over-fit the distribution and exclude regions of the parameter space which
should be sampled. At the other extreme, if the model is too simple then
resulting contour can “smooth” fine details and more samples are drawn outside
of the initial likelihood constraint. These will not be accepted and the
sampling process is therefore less efficient. We use a similar logic for the
number of neurons and layers in the neural network that parameterises the flow
but we find that these parameters predominantly affect training time with a
lesser effect on overall convergence. Another parameter that is important to
consider is the batch size, during sampling the normalising flow can be
training upwards of 100 times. Hence, a larger batch size is recommended since
it can greatly reduce training time, we also recommend increasing the batch
size when using reparameterisations that increase the amount of training data,
such as the boundary inversion described in sections III.2 and B.

We note that the size of pool of new samples effects the efficiency of the
algorithm and the total run-time. If the pool-size is small then the
normalising flow is frequently retrained, in extreme case where the proposal
is inefficient due to, for example, the complexity of the parameter space,
then the normalising can be retrained multiple times during a single
iteration. Conversely, if the pool-size is large then if the flow is force-
ably retrained a number of points are discarded or, if the flow is only
retrained once the pool is empty, then the rejection sampling becomes in-
efficient since a large fraction of the potential new points will lie outside
the likelihood bound. We instead opt to inversely scale the pool-size given
the mean acceptance of the sampler since the last iteration the flow was
trained. We recommend setting the base pool-size to the number of live points,
only retraining the model when the pool is empty and setting the maximum pool-
size to be ten times the base pool-size.

As mentioned previously, approximately 32% of the run-time is spent on
populating the pool of new samples. This is directly attributable to the
efficiency of the rejection sampling required to ensure samples are
distributed according to the prior. In section III.1 we propose two methods
for drawing samples within the contour in the latent space, these produce
uniformly and normally distributed samples respectively. In practice we find
the two methods comparable in most cases with the exception of when the latent
radius lies in the tail of the chi-distribution that corresponds to the latent
prior $p_{\mathcal{Z}}$. In this case using the uniform distribution results
in lower population and proposal acceptances which leads to longer run-times.

##  VI Conclusions

We have proposed a novel method for sampling within a given iso-likelihood
contour according to the prior that can be incorporated into the standard
nested sampling algorithm. Our method employs normalising flows to learn the
density of current set of live points which, once trained, allows us to
produce points within the contour by sampling from a simple distribution and
using rejection sampling. The use of normalising flows allows us to avoid
using multiple bounding distributions and since new samples are independent of
the previous samples we eliminate the need to use a random walk. We implement
this proposal method in our sampler, Nessai, and conduct a series of tests to
verify that it recovers the correct Bayesian posteriors and then compare our
results to those obtained with another sampler to determine if our design does
in fact result in a more efficient sampler.

We apply our sampler to 128 four second duration simulated signals from the
coalescence of binary black hole systems sampled at 2048 Hz and we run two
separate analyses, one with distance marginalisation and another without. The
resulting P-P plots (fig. 3) show that our sampler more reliably recovers the
posterior distributions with distance marginalisation than without, however
both pass the P-P test. This indicates that our proposal method does not
introduce any inherent biases.

We use dynesty for the comparison, which has been shown to produce results
consistent with those used in previous LVK analyses [32]. We find that our
sampler returns evidences consistent with dynesty, which serves as further
verification of our results. Since we aim to produce a more efficient sampler
we also compare the likelihood evaluations required to reach convergence. When
not using distance marginalisation we find that Nessai requires $5.04\times
10^{6}$ likelihood evaluations, 2.32 times fewer than dynesty. When distance
marginalisation is enabled Nessai requires $7.22\times 10^{6}$, which, whilst
still 1.40 fewer than dynesty, is more than with the marginalisation disabled.
As such, we recommend using Nessai without distance marginalisation for
gravitational-wave inference.

However, this reduction in likelihood evaluations does not translate to a
direct reduction in total computation time because of the cost associated with
training the normalising flow and populating the pool of new samples. We find
that the fraction of the time spent of each process changes when using
distance marginalisation. Without the marginalisation, on average, 8.1% of the
total computation time is spent on training and a further 39.6% on population.
When using distance marginalisation this changes to 4.9% spent on training and
42.2% on population. We attribute the difference in population time to the
efficiency of the rejection sampling, which is improved when including the
reparameterisation for distance discussed in section III.3. We find that
without distance marginalisation the median run time for Nessai is 2.07 times
faster than dynesty. However when distance marginalisation is enabled we
observe that, on average, Nessai is only 1.34 times faster than dynesty. This
further reinforces our recommendation to use Nessai with distance
marginalisation disabled.

We also show how our sampler can make use of parallelised likelihood functions
by evaluating the likelihood of new live points during the population stage.
We repeat the previous analysis for a single injection without distance
marginalisation and parallelise the likelihood computation with increasing
number of threads up to 16. We observe that the reduction time evaluating the
likelihood does not quite match the theoretical values, indicating that there
is a small overhead associated with it. This also highlights how the limiting
factor is the time spent training the normalising flow and populating the pool
of new live point.

To aid in diagnosing potential biases during sampling, we include a series of
diagnostics in our sampler which allow us to easily identify under and over-
constraining. These diagnostics also help to tune the sampling settings and
highlight how periodically re-training the normalising flow during sampling
prevents the proposal from becoming inefficient during sampling.

We find that our algorithm is susceptible to under-sampling regions of the
parameter space which are close to the prior bounds. We consequently introduce
the previously described reparameterisations to mitigate this and a series of
diagnostics to aid in diagnosing biases and correctly tuning the settings. We
aim address this in further work with changes to the design of the normalising
flows we have used.

It is natural to compare this work to [24, 25, 26, 27] which use variational
autoeconders and normalising flows to produce posterior distributions. Our
approach differs from these in that it requires no prior computation since
training occurs during sampling and we do not introduce any assumptions about
the data other than those necessary to apply a nested sampling algorithm.
Nessai is therefore a drop-in replacement for existing sampling algorithms
that does not require changes to existing pipelines.

In future work we aim to evaluate our sampler using more expensive waveform
models including those for longer duration signals, such as those from binary
neutron star of neutron star-black hole system, and models which include
higher-order modes. We will also investigate the suitability of other types of
normalising flow transforms, such as the spline based transforms from [63] and
flows which allow for specifying a manifold [64].

In summary, we have proposed a novel variation of the standard nested sampling
algorithm that incorporates normalising flows specifically designed for
inference with computationally expensive likelihood functions. We have applied
our sampler to the problem of gravitational wave inference and shown that it
consistently recovers the Bayesian posteriors distributions and evidences with
2.32 times fewer total likelihood evaluations than dynesty, another commonly
used sampler, which translates to a 2.07 times reduction in computation time.
Our sampler therefore serves as a more efficient drop-in replacement for
existing samplers.

###### Acknowledgements.

The authors gratefully acknowledge the Science and Technology Facilities
Council of the United Kingdom. MJW is supported by the Science and Technology
Facilities Council [2285031]. JV and CM are supported by the Science and
Technology Research Council [ST/ L000946/1]. CM is also supported by the
European Cooperation in Science and Technology (COST) action [CA17137]. The
authors are grateful for computational resources provided by Cardiff
University, and funded by an STFC grant supporting UK Involvement in the
Operation of Advanced LIGO. Software: Nessai was initially develop using
cpnest [18] with permission from the authors and still shares a similar
interface and other core codes. Nessai is implemented in Python and uses NumPy
[65], SciPy [66], pandas [67, 68], nflows [43], PyTorch [42], matplotlib [69]
and seaborn [70]. Gravitational wave injections were generated using Bilby and
bilby_pipe [22]. Figures were prepared using matplotlib [69], seaborn [70],
Bilby [22] and corner [71].

##  Appendix A Loss function

A normalising flow applies a mapping $f:\mathcal{X}\to\mathcal{Z}$ conditioned
on its parameters $\bm{\theta}$ which are typically the trainable parameters
of a neural network. In this context the goal of training a normalising flow
is to approximate a target distribution $p^{*}_{\mathcal{X}}(x)$. The KL
divergence between the target distribution and the distribution of the flow
$p_{\mathcal{X}}(x|\theta)$ can be written as [29]:

$\begin{split}\mathcal{L}(\bm{\theta})=&D_{\textrm{KL}}[p^{*}_{\mathcal{X}}(x)||p_{\mathcal{X}}(x|\bm{\theta})]\\\
=&-\mathbb{E}_{p^{*}_{\mathcal{X}}(x)}[\ln
p_{\mathcal{X}}(x|\bm{\theta})]+\textrm{const.}\\\
=&-\mathbb{E}_{p^{*}_{\mathcal{X}}(x)}\left[\ln
p_{\mathcal{Z}}(f(x|\bm{\theta}))+\ln\left|\textrm{det}\frac{\partial
f(x|\bm{\theta})}{\partial x^{T}}\right|\right]+\textrm{const.}.\end{split}$
(14)

Then, assuming the set of K samples used for training is drawn from
$p^{*}_{\mathcal{X}}(x)$, the expectation value can be approximated as:

$\mathcal{L}(\bm{\theta})\approx-\frac{1}{K}\sum_{k=1}^{K}\ln
p_{\mathcal{Z}}(f(x|\bm{\theta}))+\ln\left|\textrm{det}\frac{\partial
f(x|\bm{\theta})}{\partial x^{T}}\right|.$ (15)

##  Appendix B Boundary inversion

In section III.2 we describe boundary inversion which we introduce to avoid
under-sampling regions which are close to the priors bounds. The user defines
which parameters the inversion can be applied to and before training the
sampler determines if it should be applied to each parameter using the
following steps:

  1. 1.

Compute the density of samples over the specified range and find the maximum
value.

  2. 2.

Compute the fraction of the density that lies within the initial and final
$m\%$ of the bounds, i.e. $[0,0.1]$ and $[0.9,1.0]$ if the parameter is
defined on $[0,1]$.

  3. 3.

Choose to apply inversion to bound with the highest density if it is at least
$n\%$ of the maximum density and the density at the bound is non-zero.

From our testing the percentages $m$ and $n$ default to 10% and 50%
respectively but can be changed by the user. We consider two methods for
applying the inversion:

  * •

duplication: duplicate the set of points and apply the inversion to the
duplicates,

  * •

splitting: randomly select half of the points to apply the inversion to.

We find that duplication generally provides more consistent results but at the
cost of the increasing the training time. As such we recommend using splitting
when inversion is applied to more than two parameters.

##  Appendix C Gravitational-wave priors

Table 1: Prior distributions used for each parameter for gravitational-wave parameter estimation. Their corresponding labels and the lower and upper bounds are included where applicable. Parameters | Label | Prior | Bounds
---|---|---|---
Chirp mass | $\mathcal{M}$ | Uniform | $[25,35]{\text{M}_{\odot}}$
Asymmetric mass ratio | $q$ | Uniform | $[0.125,1.0]$
Luminosity distance | $d_{\text{L}}$ | Uniform in co-moving volume | $[100,2000]{\text{Mpc}}$
Right ascension | $\alpha$ | Uniform | $[0,2\pi]$
Declination | $\delta$ | Cosine | -
Reference time at geocentre | $t_{\text{c}}$ | Uniform | $[-0.1,0.1]$
Inclination | $\theta_{\text{JN}}$ | Sine | -
Polarisation | $\psi$ | Uniform | $[0,\pi]$
Phase | $\phi_{c}$ | Uniform | $[0,2\pi]$
Dimensionless spin magnitudes | $\chi_{i}$ | Uniform | $[0,0.99]$
Spin tilt angles | $\theta_{i}$ | Sine | -
Difference between the azimuthal angles of each spin vector relative to the orbital angular momentum | $\phi_{12}$ | Uniform | $[0,2\pi]$
Difference between the azimuthal angles of the total and orbital angular momentum | $\phi_{\text{JL}}$ | Uniform | $[0,2\pi]$

##  Appendix D Nessai sampling settings

Table 2: Settings used for Nessai for gravitational-wave inference. These are
split into three categories: general settings which control aspects of the
sampler such as the choice of latent prior or pool-size, flow hyper-parameters
which determine the configuration of the normalising flow and flow training
settings which control the training process. For a complete description of
each see the documentation [41]. Nessai settings
---
General settings | Flow hyper-parameters | Flow training settings
Training frequency | None | Coupling transformations | 6 | Optimiser | Adam
Cooldown | 200 | Linear transformation | LU | Learning rate | 0.001
Base pool-size | 2000 | Network type | ResNet | Batch size | 2000
Update pool-size | True | Layers per network | 2 | Max. epochs | 500
Drawsize | 2000 | Neurons per layer | 32 | Patience | 50
Train on empty | True | Activation | ReLU |  |
Weights reset | False | Batch-Normalisation | Intra-transforms |  |
Latent prior | Truncated Gaussian |  |  |  |
Rescale | True |  |  |  |
Update bounds | True |  |  |  |

##  Appendix E P-P tests for dynesty

(a) With phase marginalisation (b) With phase and distance marginalisation

Figure 9: Probability-probability (P-P) plot showing the confidence interval
versus the fraction of the events within that confidence interval for the
posterior distributions obtained using dynesty for 128 simulated compact
binary coalescence signals produced with Bilby and bilby_pipe. The 1-, 2- and
3-$\sigma$ confidence intervals are indicated by the shaded regions and
$p$-values are shown for each of the parameters and the combined $p$-value is
also shown. We use the settings described in [32] with the exception of the
number of live points which we increase to 2000.

##  Appendix F Example corner plot

Figure 10: Corner plot comparing the posterior distributions produced with
dynesty (blue) and our sampler Nessai (red). The phase is marginalised and
remaining 14 parameters are shown, see appendix C for details on the
parameters. The respective 16% and 84% quantiles are also shown in the
$1$-dimensional marginalised posteriors.

## References

  * Abbott  _et al._ [2017a] B. Abbott _et al._ (LIGO Scientific, Virgo, Fermi GBM, INTEGRAL, IceCube, AstroSat Cadmium Zinc Telluride Imager Team, IPN, Insight-Hxmt, ANTARES, Swift, AGILE Team, 1M2H Team, Dark Energy Camera GW-EM, DES, DLT40, GRAWITA, Fermi-LAT, ATCA, ASKAP, Las Cumbres Observatory Group, OzGrav, DWF (Deeper Wider Faster Program), AST3, CAASTRO, VINROUGE, MASTER, J-GEM, GROWTH, JAGWAR, CaltechNRAO, TTU-NRAO, NuSTAR, Pan-STARRS, MAXI Team, TZAC Consortium, KU, Nordic Optical Telescope, ePESSTO, GROND, Texas Tech University, SALT Group, TOROS, BOOTES, MWA, CALET, IKI-GW Follow-up, H.E.S.S., LOFAR, LWA, HAWC, Pierre Auger, ALMA, Euro VLBI Team, Pi of Sky, Chandra Team at McGill University, DFN, ATLAS Telescopes, High Time Resolution Universe Survey, RIMAS, RATIR, SKA South Africa/MeerKAT), Multi-messenger Observations of a Binary Neutron Star Merger, Astrophys. J. Lett. 848, L12 (2017a), arXiv:1710.05833 [astro-ph.HE] .
  * Abbott  _et al._ [2018] B. Abbott _et al._ (LIGO Scientific, Virgo), GW170817: Measurements of neutron star radii and equation of state, Phys. Rev. Lett. 121, 161101 (2018), arXiv:1805.11581 [gr-qc] .
  * Abbott  _et al._ [2019a] B. Abbott _et al._ (LIGO Scientific, Virgo), A gravitational-wave measurement of the Hubble constant following the second observing run of Advanced LIGO and Virgo, arXiv e-prints (2019a), arXiv:1908.06060 [astro-ph.CO] .
  * Abbott  _et al._ [2020] B. Abbott _et al._ (KAGRA, LIGO Scientific, VIRGO), Prospects for Observing and Localizing Gravitational-Wave Transients with Advanced LIGO, Advanced Virgo and KAGRA, Living Reviews in Relativity 23, 10.1007/s41114-018-0012-9 (2020), arXiv:1304.0670 [gr-qc] .
  * Iyer  _et al._ [2011] B. Iyer, T. Souradeep, C. Unnikrishnan, S. Dhurandhar, S. Raja, and A. Sengupta,  _LIGO-India, Proposal of the Consortium for Indian Initiative in Gravitational-wave Observations (IndIGO)_, Tech. Rep. (LIGO Scientific Collaboration, 2011).
  * Abbott  _et al._ [2019b] B. Abbott _et al._ (LIGO Scientific, Virgo), GWTC-1: A Gravitational-Wave Transient Catalog of Compact Binary Mergers Observed by LIGO and Virgo during the First and Second Observing Runs, Phys. Rev. X 9, 031040 (2019b), arXiv:1811.12907 [astro-ph.HE] .
  * Abbott  _et al._ [2017b] B. Abbott _et al._ (LIGO Scientific, Virgo), GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral, Phys. Rev. Lett. 119, 161101 (2017b), arXiv:1710.05832 [gr-qc] .
  * Abbott  _et al._ [2020] R. Abbott _et al._ (LIGO Scientific, Virgo), GWTC-2: Compact Binary Coalescences Observed by LIGO and Virgo During the First Half of the Third Observing Run, arXiv e-prints , arXiv:2010.14527 (2020), arXiv:2010.14527 [gr-qc] .
  * Aasi  _et al._ [2013] J. Aasi _et al._ (LIGO Scientific, VIRGO), Parameter estimation for compact binary coalescence signals with the first generation gravitational-wave detector network, Phys. Rev. D 88, 062001 (2013), arXiv:1304.1775 [gr-qc] .
  * Abbott  _et al._ [2016a] B. Abbott _et al._ (LIGO Scientific, Virgo), Properties of the Binary Black Hole Merger GW150914, Phys. Rev. Lett. 116, 241102 (2016a), arXiv:1602.03840 [gr-qc] .
  * Veitch  _et al._ [2015] J. Veitch, V. Raymond, B. Farr, W. Farr, P. Graff, S. Vitale, B. Aylott, K. Blackburn, N. Christensen, M. Coughlin, _et al._ , Parameter estimation for compact binaries with ground-based gravitational-wave observations using the lalinference software library, Phys. Rev. D 91, 042003 (2015).
  * Brooks  _et al._ [2011] S. Brooks, A. Gelman, G. Jones, and X.-L. Meng,  _Handbook of markov chain monte carlo_ (CRC press, 2011).
  * Skilling  _et al._ [2006] J. Skilling _et al._ , Nested sampling for general bayesian computation, Bayesian analysis 1, 833 (2006).
  * Farr  _et al._ [2014] B. Farr, E. Ochsner, W. M. Farr, and R. O’Shaughnessy, A more effective coordinate system for parameter estimation of precessing compact binaries from gravitational waves, Phys. Rev. D 90, 024018 (2014).
  * Pürrer [2014] M. Pürrer, Frequency domain reduced order models for gravitational waves from aligned-spin compact binaries, Class. Quant. Grav. 31, 195010 (2014), arXiv:1402.4146 [gr-qc] .
  * Smith  _et al._ [2016] R. Smith, S. E. Field, K. Blackburn, C.-J. Haster, M. Pürrer, V. Raymond, and P. Schmidt, Fast and accurate inference on gravitational waves from precessing compact binaries, Phys. Rev. D 94, 044031 (2016).
  * Handley  _et al._ [2015] W. J. Handley, M. P. Hobson, and A. N. Lasenby, polychord: next-generation nested sampling, Monthly Notices of the Royal Astronomical Society 453, 4384 (2015), https://academic.oup.com/mnras/article-pdf/453/4/4384/8034904/stv1911.pdf .
  * Veitch  _et al._ [2021] J. Veitch, W. D. Pozzo, M. Williams, C. Talbot, M. Pitkin, G. Ashton, Cody, M. Hübner, A. Nitz, D. Macleod, G. Carullo, G. Davies, and Tony, johnveitch/cpnest: Version 0.11 (2021).
  * Smith  _et al._ [2020] R. J. E. Smith, G. Ashton, A. Vajpeyi, and C. Talbot, Massively parallel Bayesian inference for transient gravitational-wave astronomy, Monthly Notices of the Royal Astronomical Society 498, 4492 (2020), https://academic.oup.com/mnras/article-pdf/498/3/4492/33798799/staa2483.pdf .
  * Lange  _et al._ [2018] J. Lange, R. O’Shaughnessy, and M. Rizzo, Rapid and accurate parameter inference for coalescing, precessing compact binaries, arXiv e-prints , arXiv:1805.10457 (2018), arXiv:1805.10457 [gr-qc] .
  * Biwer  _et al._ [2019] C. M. Biwer, C. D. Capano, S. De, M. Cabero, D. A. Brown, A. H. Nitz, and V. Raymond, PyCBC Inference: A Python-based Parameter Estimation Toolkit for Compact Binary Coalescence Signal, PASP 131, 024503 (2019), arXiv:1807.10312 [astro-ph.IM] .
  * Ashton  _et al._ [2019] G. Ashton, M. Hübner, P. D. Lasky, C. Talbot, K. Ackley, S. Biscoveanu, Q. Chu, A. Divakarla, P. J. Easter, B. Goncharov, F. Hernandez Vivanco, J. Harms, M. E. Lower, G. D. Meadors, D. Melchor, E. Payne, M. D. Pitkin, J. Powell, N. Sarin, R. J. E. Smith, and E. Thrane, BILBY: A User-friendly Bayesian Inference Library for Gravitational-wave Astronomy, ApJS 241, 27 (2019), arXiv:1811.02042 [astro-ph.IM] .
  * Cuoco  _et al._ [2020] E. Cuoco, J. Powell, M. Cavaglià, K. Ackley, M. Bejger, C. Chatterjee, M. Coughlin, S. Coughlin, P. Easter, R. Essick, H. Gabbard, T. Gebhard, S. Ghosh, L. Haegel, A. Iess, D. Keitel, Z. Marka, S. Marka, F. Morawski, T. Nguyen, R. Ormiston, M. Puerrer, M. Razzano, K. Staats, G. Vajente, and D. Williams, Enhancing Gravitational-Wave Science with Machine Learning, arXiv e-prints , arXiv:2005.03745 (2020), arXiv:2005.03745 [astro-ph.HE] .
  * Gabbard  _et al._ [2019] H. Gabbard, C. Messenger, I. S. Heng, F. Tonolini, and R. Murray-Smith, Bayesian parameter estimation using conditional variational autoencoders for gravitational-wave astronomy, arXiv e-prints , arXiv:1909.06296 (2019), arXiv:1909.06296 [astro-ph.IM] .
  * Chua and Vallisneri [2020] A. J. K. Chua and M. Vallisneri, Learning Bayesian Posteriors with Neural Networks for Gravitational-Wave Inference, Phys. Rev. Lett. 124, 041102 (2020), arXiv:1909.05966 [gr-qc] .
  * Green  _et al._ [2020] S. R. Green, C. Simpson, and J. Gair, Gravitational-wave parameter estimation with autoregressive neural network flows, arXiv e-prints , arXiv:2002.07656 (2020), arXiv:2002.07656 [astro-ph.IM] .
  * Green and Gair [2020] S. R. Green and J. Gair, Complete parameter inference for GW150914 using deep learning, arXiv e-prints , arXiv:2008.03312 (2020), arXiv:2008.03312 [astro-ph.IM] .
  * Kobyzev  _et al._ [2019] I. Kobyzev, S. J. D. Prince, and M. A. Brubaker, Normalizing Flows: An Introduction and Review of Current Methods, arXiv e-prints , arXiv:1908.09257 (2019), arXiv:1908.09257 [stat.ML] .
  * Papamakarios  _et al._ [2019] G. Papamakarios, E. Nalisnick, D. Jimenez Rezende, S. Mohamed, and B. Lakshminarayanan, Normalizing Flows for Probabilistic Modeling and Inference, arXiv e-prints , arXiv:1912.02762 (2019), arXiv:1912.02762 [stat.ML] .
  * Feroz  _et al._ [2009] F. Feroz, M. P. Hobson, and M. Bridges, MULTINEST: an efficient and robust Bayesian inference tool for cosmology and particle physics, MNRAS 398, 1601 (2009), arXiv:0809.3437 [astro-ph] .
  * Speagle [2020] J. S. Speagle, DYNESTY: a dynamic nested sampling package for estimating Bayesian posteriors and evidences, MNRAS 493, 3132 (2020), arXiv:1904.02180 [astro-ph.IM] .
  * Romero-Shaw  _et al._ [2020] I. M. Romero-Shaw, C. Talbot, S. Biscoveanu, V. D’Emilio, G. Ashton, C. P. L. Berry, S. Coughlin, S. Galaudage, C. Hoy, M. Huebner, K. S. Phukon, M. Pitkin, M. Rizzo, N. Sarin, R. Smith, S. Stevenson, A. Vajpeyi, M. Arene, K. Athar, S. Banagiri, N. Bose, M. Carney, K. Chatziioannou, R. Cotesta, B. Edelman, C. Garcia-Quiros, A. Ghosh, R. Green, C. J. Haster, A. X. Kim, F. Hernandez-Vivanco, I. Magana Hernandez, C. Karathanasis, P. D. Lasky, N. De Lillo, M. E. Lower, D. Macleod, M. Mateu-Lucena, A. Miller, M. Millhouse, S. Morisaki, S. H. Oh, S. Ossokine, E. Payne, J. Powell, M. Puerrer, A. Ramos-Buades, V. Raymond, E. Thrane, J. Veitch, D. Williams, M. J. Williams, and L. Xiao, Bayesian inference for compact binary coalescences with BILBY: Validation and application to the first LIGO–Virgo gravitational-wave transient catalogue, arXiv e-prints , arXiv:2006.00714 (2020), arXiv:2006.00714 [astro-ph.IM] .
  * Jimenez Rezende and Mohamed [2015] D. Jimenez Rezende and S. Mohamed, Variational Inference with Normalizing Flows, arXiv e-prints , arXiv:1505.05770 (2015), arXiv:1505.05770 [stat.ML] .
  * Papamakarios  _et al._ [2017] G. Papamakarios, T. Pavlakou, and I. Murray, Masked Autoregressive Flow for Density Estimation, arXiv e-prints , arXiv:1705.07057 (2017), arXiv:1705.07057 [stat.ML] .
  * Huang  _et al._ [2018] C.-W. Huang, D. Krueger, A. Lacoste, and A. Courville, Neural Autoregressive Flows, arXiv e-prints , arXiv:1804.00779 (2018), arXiv:1804.00779 [cs.LG] .
  * Dinh  _et al._ [2016] L. Dinh, J. Sohl-Dickstein, and S. Bengio, Density estimation using Real NVP, arXiv e-prints , arXiv:1605.08803 (2016), arXiv:1605.08803 [cs.LG] .
  * Kingma and Dhariwal [2018] D. P. Kingma and P. Dhariwal, Glow: Generative Flow with Invertible 1x1 Convolutions, arXiv e-prints , arXiv:1807.03039 (2018), arXiv:1807.03039 [stat.ML] .
  * Durkan  _et al._ [2019a] C. Durkan, A. Bekasov, I. Murray, and G. Papamakarios, Neural Spline Flows, arXiv e-prints , arXiv:1906.04032 (2019a), arXiv:1906.04032 [stat.ML] .
  * Muller [1959] M. E. Muller, A note on a method for generating points uniformly on n-dimensional spheres, Commun. ACM 2, 19 (1959).
  * Marsaglia [1972] G. Marsaglia, Choosing a point from the surface of a sphere, Ann. Math. Statist. 43, 645 (1972).
  * Williams [2021a] M. J. Williams, Nessai documentation, https://nessai.readthedocs.io/ (2021a).
  * Paszke  _et al._ [2019] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala, Pytorch: An imperative style, high-performance deep learning library, in  _Advances in Neural Information Processing Systems 32_, edited by H. Wallach, H. Larochelle, A. Beygelzimer, F. d’Alché Buc, E. Fox, and R. Garnett (Curran Associates, Inc., 2019) pp. 8024–8035.
  * Durkan  _et al._ [2020] C. Durkan, A. Bekasov, I. Murray, and G. Papamakarios, nflows: normalizing flows in PyTorch (2020).
  * Dinh  _et al._ [2014] L. Dinh, D. Krueger, and Y. Bengio, NICE: Non-linear Independent Components Estimation, arXiv e-prints , arXiv:1410.8516 (2014), arXiv:1410.8516 [cs.LG] .
  * Ioffe and Szegedy [2015] S. Ioffe and C. Szegedy, Batch Normalization: Accelerating Deep Network Training by Reducing Internal Covariate Shift, arXiv e-prints , arXiv:1502.03167 (2015), arXiv:1502.03167 [cs.LG] .
  * He  _et al._ [2015] K. He, X. Zhang, S. Ren, and J. Sun, Deep Residual Learning for Image Recognition, arXiv e-prints , arXiv:1512.03385 (2015), arXiv:1512.03385 [cs.CV] .
  * He  _et al._ [2016] K. He, X. Zhang, S. Ren, and J. Sun, Identity Mappings in Deep Residual Networks, arXiv e-prints , arXiv:1603.05027 (2016), arXiv:1603.05027 [cs.CV] .
  * Kingma and Ba [2014] D. P. Kingma and J. Ba, Adam: A method for stochastic optimization, arXiv preprint arXiv:1412.6980 (2014).
  * Williams [2021b] M. J. Williams, Nessai: Nested sampling with artificial intelligence (2021b).
  * Abbott  _et al._ [2016b] B. Abbott _et al._ (LIGO Scientific, Virgo), Observation of Gravitational Waves from a Binary Black Hole Merger, Phys. Rev. Lett. 116, 061102 (2016b), arXiv:1602.03837 [gr-qc] .
  * Graff and Feroz [2013] P. Graff and F. Feroz, BAMBI: Blind Accelerated Multimodal Bayesian Inference (2013), ascl:1312.008 .
  * Levy  _et al._ [2017] D. Levy, M. D. Hoffman, and J. Sohl-Dickstein, Generalizing Hamiltonian Monte Carlo with Neural Networks, arXiv e-prints , arXiv:1711.09268 (2017), arXiv:1711.09268 [stat.ML] .
  * Hoffman  _et al._ [2019] M. Hoffman, P. Sountsov, J. V. Dillon, I. Langmore, D. Tran, and S. Vasudevan, NeuTra-lizing Bad Geometry in Hamiltonian Monte Carlo Using Neural Transport, arXiv e-prints , arXiv:1903.03704 (2019), arXiv:1903.03704 [stat.CO] .
  * Moss [2019] A. Moss, Accelerated Bayesian inference using deep learning, arXiv e-prints , arXiv:1903.10860 (2019), arXiv:1903.10860 [astro-ph.CO] .
  * Schmidt  _et al._ [2012] P. Schmidt, M. Hannam, and S. Husa, Towards models of gravitational waveforms from generic binaries: A simple approximate mapping between precessing and nonprecessing inspiral signals, Phys. Rev. D 86, 104063 (2012).
  * Khan  _et al._ [2019] S. Khan, K. Chatziioannou, M. Hannam, and F. Ohme, Phenomenological model for the gravitational-wave signal from precessing binary black holes with two-spin effects, Phys. Rev. D 100, 024059 (2019), arXiv:1809.10113 [gr-qc] .
  * Cook  _et al._ [2006] S. R. Cook, A. Gelman, and D. B. Rubin, Validation of software for bayesian models using posterior quantiles, Journal of Computational and Graphical Statistics 15, 675 (2006), https://doi.org/10.1198/106186006X136976 .
  * Talts  _et al._ [2018] S. Talts, M. Betancourt, D. Simpson, A. Vehtari, and A. Gelman, Validating Bayesian Inference Algorithms with Simulation-Based Calibration, arXiv e-prints , arXiv:1804.06788 (2018), arXiv:1804.06788 [stat.ME] .
  * Veitch and Vecchio [2010] J. Veitch and A. Vecchio, Bayesian coherent analysis of in-spiral gravitational wave signals with a detector network, Phys. Rev. D 81, 062003 (2010), arXiv:0911.3820 [astro-ph.CO] .
  * Fowlie  _et al._ [2020] A. Fowlie, W. Handley, and L. Su, Nested sampling cross-checks using order statistics, arXiv e-prints , arXiv:2006.03371 (2020), arXiv:2006.03371 [stat.CO] .
  * Smirnov [1948] N. Smirnov, Table for estimating the goodness of fit of empirical distributions, The annals of mathematical statistics 19, 279 (1948).
  * Arnold and Emerson [2011] T. B. Arnold and J. W. Emerson, Nonparametric goodness-of-fit tests for discrete null distributions., R Journal 3 (2011).
  * Durkan  _et al._ [2019b] C. Durkan, A. Bekasov, I. Murray, and G. Papamakarios, Neural Spline Flows, arXiv e-prints , arXiv:1906.04032 (2019b), arXiv:1906.04032 [stat.ML] .
  * Brehmer and Cranmer [2020] J. Brehmer and K. Cranmer, Flows for simultaneous manifold learning and density estimation, arXiv e-prints , arXiv:2003.13913 (2020), arXiv:2003.13913 [stat.ML] .
  * van der Walt  _et al._ [2011] S. van der Walt, S. C. Colbert, and G. Varoquaux, The numpy array: A structure for efficient numerical computation, Computing in Science Engineering 13, 22 (2011).
  * Virtanen  _et al._ [2020] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. Jarrod Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. Carey, İ. Polat, Y. Feng, E. W. Moore, J. Vand erPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and S. . . Contributors, SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python, Nature Methods 17, 261 (2020).
  * pandas development team [2020] T. pandas development team, pandas-dev/pandas: Pandas (2020).
  * Wes McKinney [2010] Wes McKinney, Data Structures for Statistical Computing in Python, in  _Proceedings of the 9th Python in Science Conference_, edited by Stéfan van der Walt and Jarrod Millman (2010) pp. 56 – 61.
  * Hunter [2007] J. D. Hunter, Matplotlib: A 2d graphics environment, Computing in Science & Engineering 9, 90 (2007).
  * Waskom and the seaborn development team [2020] M. Waskom and the seaborn development team, mwaskom/seaborn (2020).
  * Foreman-Mackey [2016] D. Foreman-Mackey, corner.py: Scatterplot matrices in python, The Journal of Open Source Software 1, 24 (2016).

  *[CBC]: compact binary coalescence
  *[GWs]: gravitational wave
  *[MCMC]: Markov Chain Monte Carlo
  *[i.i.d.]: independently and identically distributed
  *[GW]: gravitational wave