# Peregrine: Sequential simulation-based inference for gravitational wave
signals

Uddipta Bhardwaj<EMAIL_ADDRESS>GRAPPA Institute, Anton Pannekoek
Institute for Astronomy and Institute of High-Energy Physics,
University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands
James Alvey<EMAIL_ADDRESS>GRAPPA Institute, Institute for Theoretical
Physics Amsterdam,
University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands
Benjamin Kurt Miller  GRAPPA Institute, Institute for Theoretical Physics
Amsterdam,
University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands
AMLab, Institute for Informatics, LAB42, Science Park 900, 1098 XH Amsterdam
AI4Science, Institute for Informatics, LAB42, Science Park 900, 1098 XH
Amsterdam  Samaya Nissanke  GRAPPA Institute, Anton Pannekoek Institute for
Astronomy and Institute of High-Energy Physics,
University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands
Nikhef, Science Park 105, 1098 XG Amsterdam, The Netherlands  Christoph
Weniger  GRAPPA Institute, Institute for Theoretical Physics Amsterdam,
University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands

###### Abstract

The current and upcoming generations of gravitational wave experiments
represent an exciting step forward in terms of detector sensitivity and
performance. For example, key upgrades at the LIGO, Virgo and KAGRA facilities
will see the next observing run (O4) probe a spatial volume around four times
larger than the previous run (O3), and design implementations for e.g. the
Einstein Telescope, Cosmic Explorer and LISA experiments are taking shape to
explore a wider frequency range and probe cosmic distances. In this context,
however, a number of very real data analysis problems face the gravitational
wave community. For example, it will be crucial to develop tools and
strategies to analyse (amongst other scenarios) signals that arrive
coincidentally in detectors, longer signals that are in the presence of non-
stationary noise or other shorter transients, as well as noisy, potentially
correlated, coherent stochastic backgrounds. With these challenges in mind, we
develop peregrine, a new sequential simulation-based inference approach
designed to study broad classes of gravitational wave signal. In this work, we
describe the method and implementation, before demonstrating its accuracy and
robustness through direct comparison with established likelihood-based
methods. Specifically, we show that we are able to fully reconstruct the
posterior distributions for every parameter of a spinning, precessing compact
binary coalescence using one of the most physically detailed and
computationally expensive waveform approximants (SEOBNRv4PHM). Crucially, we
are able to do this using only 2% of the waveform evaluations that are
required in _e.g._ nested sampling approaches. Finally, we provide some
outlook as to how this level of simulation efficiency and flexibility in the
statistical analysis could allow peregrine to tackle these current and future
gravitational wave data analysis problems.

GitHub: The peregrine analysis and inference library will be made available
here.

##  I Introduction

Figure 1: A schematic illustration of the data analysis method developed in
this paper. We use Truncated Marginal Neural Ratio Estimation (TMNRE) to
perform full parameter estimation on gravitational wave signals. The algorithm
is implemented in our new code peregrine.

_Observational Status._ The gravitational wave (GW) sky is now much louder
than it was after the first detection of a gravitational wave signal from the
merger of two black holes in 2015 Abbott  _et al._ (2016). Indeed, source
catalogs are now complete and varied enough to be used for studying questions
in gravitational theory Abbott  _et al._ (2021a), cosmology Abbott  _et al._
(2021b); Nissanke  _et al._ (2010, 2011), and even the formation and internal
properties of black holes and neutron stars Abbott  _et al._ (2021c, 2018a).
To give some current context, the LIGO-Virgo collaboration (LVKC) has now
reported on 90 compact binary coalescence events111Additional catalogs such as
the deep extended catalog Abbott  _et al._ (2021d), the 3-OGC catalog Nitz
_et al._ (2021) and others Venumadhav  _et al._ (2019a, 2020) report several
other sources with different analysis pipelines or higher false alarm rate
thresholds. that are believed to be of astrophysical origin Abbott  _et al._
(2021e, f).

Along with more efficient search pipelines Allen  _et al._ (2012); Isi  _et
al._ (2018); Venumadhav  _et al._ (2019b), one of the major reasons for the
increase in size of the gravitational wave catalog is the large increase in
search volume probed by the detectors. For example, in the previous LVKC
observing run (O3), the projected comoving search volume for binary black hole
(BBH) mergers was $V_{\mathrm{BBH}}=3.4\,\times 10^{8}\,\mathrm{Mpc}^{3}$
Abbott  _et al._ (2018b). In comparison, however, the next observing run (O4)
that is planned for May 2023 has a projected search volume of around
$V_{\mathrm{BBH}}=1.5\,\mathrm{Gpc}^{3}$ Abbott  _et al._ (2018b); Petrov  _et
al._ (2022), representing over a 400% increase in volume and the corresponding
expected event rate222The equivalent volumes for binary neutron star (BNS)
coalescence events are $V_{\mathrm{BNS}}=3.3\,\times 10^{6}\,\mathrm{Mpc}^{3}$
and $V_{\mathrm{BNS}}=3.4\,\times 10^{8}\,\mathrm{Mpc}^{3}$ for O3 and O4
respectively Abbott  _et al._ (2018b); Petrov  _et al._ (2022).. This
observing run is predicted to have the advanced LIGO (aLIGO) detectors running
close to their design sensitivities (O5), and therefore, represents a key
benchmark for the current generation of gravitational wave observatories. Such
an improvement in detector sensitivity and corresponding event rate induces a
new era of computational challenges when it comes to gravitational wave data
analysis.

_Gravitational Wave Data Analysis._ In general, the data analysis pipeline for
gravitational waves splits into two parts: detection and parameter inference.
Whilst studying the optimality of search pipelines Allen  _et al._ (2012); Isi
_et al._ (2018), especially for new classes of signal, is an interesting and
important problem, in this work we will focus on the follow up step of high
precision parameter inference. The traditional (and current) approach to
gravitational wave parameter inference follows a Bayesian framework where the
goal is to infer either the full ($N$-dimensional) joint posterior over all
parameters or marginal posterior distributions over a (sub-)set of relevant
intrinsic (_e.g._ source properties) and extrinsic parameters (_e.g._ the
source-detector geometry) Cutler and Flanagan (1994); Thrane and Talbot
(2019). In addition, one must make some informed choice (ideally using some
astrophysically motivated range, or making a maximally uninformative choice in
the absence of this) of the prior distributions for the parameters of
interest. For current state-of-the-art waveform approximants/generators, which
are now capable of describing complex scenarios such as fully precessing and
spinning compact binary inspirals Ossokine  _et al._ (2020); Khan  _et al._
(2016); Pratten  _et al._ (2020); Edwards  _et al._ (2023), the number of
relevant parameters is around $N_{\mathrm{params}}\sim 12$–$17$. At a purely
practical level, this means that approaches such as profiling the
gravitational wave data likelihood are computationally infeasible, and it is
common therefore to use some form of stochastic sampling method to investigate
the parameter space.

Traditional stochastic algorithms such as Markov Chain Monte Carlo (MCMC)
Mackay (2003); Foreman-Mackey  _et al._ (2013) and nested sampling Skilling
(2006); Handley  _et al._ (2015); Ashton  _et al._ (2022); Speagle (2020) are
both implemented within the standard open source analysis software such as
LALInference Veitch  _et al._ (2015), PyCBC Biwer  _et al._ (2019), and bilby
Ashton  _et al._ (2019); Romero-Shaw  _et al._ (2020); Ashton and Talbot
(2021); Williams  _et al._ (2023). Both MCMC and nested sampling algorithms
function by generating samples from the full joint posterior (although they
achieve this in different ways), and then typically marginalise to obtain
posteriors over individual, or pairs of parameters. The computational cost in
these cases comes about since these algorithms carry out numerous likelihood
evaluations for each individual sample point. Indeed, a general rule of thumb
for current gravitational wave parameter inference is that generating around
$10^{4}$ independent posterior samples requires at least $10^{6}$ likelihood
evaluations333This is especially true for spinning and precessing compact
binary inspirals. (and corresponding forward evaluations of the waveform
generation model) Ashton and Talbot (2021); Williams  _et al._ (2023).
Depending on the signal duration and complexity of the waveform approximant,
this can result in inference times of hours to weeks for a single events, see
e.g. Ashton and Talbot (2021); Williams  _et al._ (2023). This is the context
in which concrete data analysis challenges will arise as, amongst other
things, _(i)_ the detection rate increases significantly Abbott  _et al._
(2018b), _(ii)_ more realistic waveform approximants that incorporate _e.g._
higher order modes become more expensive to evaluate Iacovelli  _et al._
(2022); Pratten  _et al._ (2020), and _(iii)_ future gravitational wave
detections involving multiple overlapping waveforms Pizzati  _et al._ (2022);
Antonelli  _et al._ (2021); Samajdar  _et al._ (2021); Janquart  _et al._
(2022); Langendorff  _et al._ (2022), non-stationary noise distributions, or
some stochastic background, see _e.g._ Caprini and Figueroa (2018);
Christensen (2019); Abbott  _et al._ (2019a), will render explicit likelihood
evaluations increasingly complicated and costly.

_Simulation-based Inference._ Recently, there have been significant advances
in the field of simulation-based inference (SBI) Cranmer  _et al._ (2020);
Brehmer and Cranmer (2020); Lueckmann  _et al._ (2021) (sometimes also known
as implicit likelihood inference), partly as a result of rapid developments in
machine learning, but also as a response to emerging data analysis challenges
such as those described above. The field of simulation-based inference is now
wide and varied, with numerous approaches available as fully implemented open-
source software Miller  _et al._ (2021, 2022a); Tejero-Cantero  _et al._
(2020); Alsing  _et al._ (2019). It has been shown to be successful in a
number of different physics settings, such as CMB analyses Cole  _et al._
(2022), strong lensing image analysis Montel  _et al._ (2022), point source
searches Anau Montel and Weniger (2022), and field-level cosmology Makinen
_et al._ (2021), along with other examples Dimitriou  _et al._ (2022); Gagnon-
Hartman  _et al._ (2023); Delaunoy  _et al._ (2020); Karchev  _et al._ (2022);
Lin  _et al._ (2022). Broadly, the key advantages of this set of methods are
_(i)_ they can be highly simulation-efficient compared to traditional methods,
and _(ii)_ they are “simulation-based” in the sense that they do not require
an explicit likelihood to be written down, only a realistic forward simulator
to be provided. Combined, this opens up the exciting possibility of using
higher fidelity forward models to take full advantage of even the most complex
data. In the context of gravitational waves, a particular version of
simulation-based inference known as neural posterior estimation (NPE) has
already been successfully applied to perform fully amortised
inference444Amortized inference refers to training a neural network to
approximate the posterior distribution over model parameters which allows for
fast inference on any input data. on compact binary mergers mergers Dax  _et
al._ (2021a); Wildberger  _et al._ (2022). This leads to impressive inference
performance across _e.g._ the LVKC catalog, where after the complex training
process (requiring millions of waveform evaluations) has been carried out,
almost real-time inference can be done Dax  _et al._ (2021b). It is worth
noting in this regard that the focus of our analysis pipeline is different. We
look to perform analysis on individual gravitational wave signals, where this
sort of global amortisation becomes infeasible. The prototypical example of
this is would be the overlapping waveform scenario Pizzati  _et al._ (2022);
Antonelli  _et al._ (2021); Samajdar  _et al._ (2021); Langendorff  _et al._
(2022); Janquart  _et al._ (2022), where the number of waveform calls required
to perform reliable amortised inference would likely be at least an order of
magnitude larger, if not more. In this current work, we however focus on
single BBH mergers to develop our method because this is where clear, robust
quantitative comparisons can be performed to validate our approach.

_Key Contributions._ In this work, we develop a new data analysis method that
is applicable to wide classes of gravitational wave signals. In particular, we
implement a particular, sequential simulation-based inference algorithm known
as Truncated Marginal Neural Ratio Estimation (TMNRE) Miller  _et al._ (2022a)
based on the swyft software Miller  _et al._ (2021). We demonstrate the
applicability, accuracy and robustness of this approach by studying two case
studies – a highly spinning, precessing binary black hole system at distances
of $200\,\mathrm{Mpc}$ (extremely high signal-to-noise ratio (SNR) of $\sim
100$ to act as a stringent test of the precision we can achieve) and
$900\,\mathrm{Mpc}$ (low SNR of $\sim 20$ to test our ability to reliably
distinguish between signal and noise). Using a state-of-the-art waveform
approximant SEOBNRv4PHM Ossokine  _et al._ (2020) (although in practice, any
waveform approximant for e.g. IMRPhenomXPHM Pratten  _et al._ (2021),
IMRPhenomPV2 or IMRPhenomD Hannam  _et al._ (2014); Khan  _et al._ (2016), can
easily be used), we demonstrate that we can fully reproduce the posterior
distributions generated by current likelihood-based methods. Importantly, we
show that we can achieve this agreement _with only 2% of the waveform
evaluations_ that are needed in standard likelihood-based approaches. We also
discuss the robustness tests that can be carried out to validate our results
independent of any other method. In addition, we will release our
implementation, which we call peregrine, as open-source, extendable analysis
software available to the community.

_Structure of the Work._ The remainder of the work is structured as follows:
we begin in Sec. II by briefly reviewing simulation-based inference approaches
before giving a detailed explanation of TMNRE and its application to
gravitational wave analysis. We apply this method in Sec. III to the two case
studies detailed above. Finally, in Sec. IV, we provide some discussion and
outline our conclusions regarding the application of our pipeline to current
and future gravitational wave data analysis challenges.

##  II Simulation-based Inference for Gravitational Waves

In this section, we will give a very brief overview of the general class of
simulation-based inference methods, before focusing our attention on the
particular implementation we use in this work. Finally, we will discuss some
specific design choices that are relevant to the analysis of gravitational
waves.

###  II.1 Overview of Simulation-based Inference

Over the last three to four years, there has been a significant increase in
the usage and development of so-called simulated-based inference approaches to
data analysis Cranmer  _et al._ (2020); Brehmer and Cranmer (2020)555As
discussed in the introduction, this is arguably for two main reasons, the
first is that the huge computational improvements in machine learning (ML)
methods allow for density estimators, classifiers and neural networks to be
efficiently trained on very high dimensional data. The second reason is then
that utilising these developments opens up the possibility to tackle new data
analysis challenges Pizzati  _et al._ (2022); Babak  _et al._ (2010); Schäfer
_et al._ (2023).. Whilst there are now multiple options in terms of particular
implementations and algorithms, every approach has a common setup. In
particular, the main question driving all SBI algorithms is: _can we still do
robust Bayesian inference if we are only given some generative model?_ Cranmer
_et al._ (2020). In other words, given some forward model $p(x,\theta)$ that
takes some underlying set of parameters $\theta$ and produces some simulated
data $x$, can we construct meaningful posteriors $p(\theta|x)$ for all, or
some subset, of the parameters $\theta$. In practice, how this is done is by
studying Bayes’ theorem applied to the posterior,

$p(\theta|x)=\frac{p(x|\theta)p(\theta)}{p(x)},$ (1)

where $p(\theta|x)$ is the posterior of $\theta$ given some observed or
simulated data $x$, $p(x|\theta)$ is the likelihood of a given set of data $x$
given some input parameters $\theta$, $p(\theta)$ is the (Bayesian) prior over
the parameters $\theta$, and $p(x)$ is the Bayesian evidence. The key
functionality of simulation-based inference is to realise that having a
forward generative model $p(x,\theta)=p(x|\theta)\,p(\theta)$ is equivalent to
being able to sample from the (simulated) likelihood. This is the origin of
the terms “likelihood-free” or “implicit likelihood” inference Cranmer  _et
al._ (2020).

At this point, the various SBI methods diverge somewhat in terms of how they
use this capability of sampling from the likelihood to construct posterior
densities or posterior samples. Broadly, they can be categorised in terms of
how they split up Bayes’ theorem in Eq. (1). There are three main ways to do
this given a prior $p(\theta)$ (which is explicitly known/chosen) and a
generative model $p(x,\theta)=p(x|\theta)\,p(\theta)$.

  * •

_Neural Posterior Estimation (NPE)._ The first class of methods attempts to
directly estimate the posterior density $p(\theta|x)$, typically through the
use of some flexible density estimator such as a normalising flow Papamakarios
and Murray (2016); Dax  _et al._ (2021a); Tejero-Cantero  _et al._ (2020);
Zeghal  _et al._ (2022). This has been used in a number of contexts, including
for the amortised analysis of gravitational waves Dax  _et al._ (2021b);
Wildberger  _et al._ (2022); Dax  _et al._ (2022).

  * •

_Neural Likelihood Estimation (NLE)._ In the second case, given some training
data generated from the forward model, one attempts to create a well-
controlled proxy for the simulated data likelihood $p(x|\theta)$ Papamakarios
and Murray (2016); Alsing  _et al._ (2019); Lin  _et al._ (2022). This can
then be used directly in traditional algorithms such as MCMC to generate
posterior samples.

  * •

_Neural Ratio Estimation (NRE)._ The final class of methods takes a different
approach and estimates the likelihood-to-evidence ratio $p(x|\theta)/p(x)$
Miller  _et al._ (2021); Rozet and Louppe (2021); Montel  _et al._ (2022);
Karchev  _et al._ (2022); Anau Montel and Weniger (2022); Gagnon-Hartman  _et
al._ (2023); Delaunoy  _et al._ (2022); Cole  _et al._ (2022); Dimitriou  _et
al._ (2022); Miller  _et al._ (2022b); Hermans  _et al._ (2019); Durkan  _et
al._ (2020). In this work, we will be using a particular version of this known
as Truncated Marginal Neural Ratio Estimation (TMNRE) Miller  _et al._
(2022a). Ratio estimation ultimately works by translating the Bayesian
inference problem into a binary classification task on joint and marginal
samples. We provide more details of this mapping below in the specific context
of gravitational waves.

###  II.2 The TMNRE Algorithm

In this work, we use a particular simulation-based inference algorithm known
as Truncated Marginal Neural Ratio Estimation (TMNRE) Miller  _et al._
(2022a), implemented within the swyft software package Miller  _et al._
(2021). There are two key features of the method that we look to take
advantage of for the highly simulation efficient analysis of gravitational
wave signals. For reference, the overall structure of the algorithm is
described in Fig. 1.

  * •

_The “T”._ TMNRE is a sequential method in the sense that the analysis is
performed in rounds (typically around 5 to 10). At the end of each round,
estimates of the posterior are generated for a _specific_ , targeted
observation $x_{0}$. The priors are then truncated before generating training
data for the next round, leading to a significant (empirically, of the order
of 50% to 90%) reduction in simulation budget Cole  _et al._ (2022). This, of
course, has the consequence that we are generally not developing an amortised
method, but rather one that can be applied to individual observations.

  * •

_The “M”._ The marginal aspect of TMNRE relates to the fact that the algorithm
estimates _marginal posteriors_. In other words, instead of estimating the
full ($N$-dimensional) joint posterior $p(\boldsymbol{\theta}|x)$ for all
parameters $\boldsymbol{\theta}\equiv(\theta_{1},\ldots,\theta_{N})$, and then
marginalising, we directly estimate the marginal posterior
$p(\theta_{i},\theta_{j},\cdots|x)$ for some parameter (or set of parameters)
of interest $\left\\{\theta_{i},\theta_{j},\ldots\right\\}$. This again has
the advantage that it typically reduces the simulation budget by a large
fraction and results in much simpler training procedures because the data-
parameter manifolds for the low dimensional marginals are much less complex.

More technically, TMNRE sets up the estimation of the ratio
$r(x;\vartheta=(\theta_{i},\theta_{j},\ldots))=p(x|\vartheta)/p(x)$, where
$\vartheta$ represents a single parameter or some subset of parameters, as a
binary classification problem through the following observation about Bayes’
theorem,

$r(x;\vartheta)=\frac{p(x|\vartheta)}{p(x)}=\frac{p(\vartheta|x)}{p(\vartheta)}=\frac{p(x,\vartheta)}{p(x)p(\vartheta)}.$
(2)

In other words, we can view this likelihood-to-evidence ratio in multiple
ways, _(i)_ as the posterior-to-prior ratio, which will ultimately allow us to
generate posterior samples by generating samples from the prior $p(\vartheta)$
and weighting them by the ratio $r(x;\vartheta)$, and _(ii)_ as the ratio
between the joint probability density $p(x,\vartheta)$ and marginal density
$p(x)p(\vartheta)$.

It is the last form that allows us to carry out a binary classification task.
In particular, suppose we have $N$ samples from our simulator (or joint
distribution $p(x,\vartheta)=p(x|\vartheta)p(\vartheta)$)
$\left\\{(\vartheta=(\theta_{i},\theta_{j},\ldots),x)_{i}\right\\}_{k=1...N}$,
then we can also construct samples from the marginal distribution
$p(x)p(\vartheta)$ by randomly shuffling the pair components. The method then
attempts to optimise a classifier $d_{\phi}(x,\vartheta)$ with some trainable
parameters $\phi$ to output a class $d_{\phi}(x,\vartheta)=0$ (say) when
$(x,\vartheta)$ is drawn marginally, and $d_{\phi}(x,\vartheta)=1$ when
$(x,\vartheta)$ is drawn jointly.

Figure 2: Two example observations that will be used in the case studies
below. The parameter values for the observations can be found in the inference
results below (see Figs. 4 and 5) _Left:_ A (typical) low signal-to-noise
ratio (SNR) signal as seen in the LIGO Hanford detector. The signal is shown
in pink, while the full detector response including noise is shown in grey.
_Right:_ A similar signal (shown in orange), but with a very high SNR, again
in the Hanford detector. _Upper and lower panels:_ The upper panels illustrate
the time-domain strain $h(t)$, whilst the lower panels show the corresponding
frequency domain response $\tilde{h}(f)$, both of which are used in our
analysis.

In terms of TMNRE, this is technically achieved by minimising the (binary
cross-entropy) loss function Miller  _et al._ (2021, 2022b); Delaunoy  _et
al._ (2022),

$\mathcal{L}[f_{\phi}]=-\int{\mathrm{d}x\mathrm{d}\vartheta}\,[p(x,\vartheta)\ln\left(\sigma(f_{\phi}(x,\vartheta))\right)\\\
+p(x)p(\vartheta)\ln\left(1-\sigma(f_{\phi}(x,\vartheta))\right)],$ (3)

where $\sigma(x)=[1+\exp(-x)]^{-1}$ is the sigmoid function and
$d_{\phi}(x,\vartheta)=\sigma(f_{\phi}(x,\vartheta))$ in the notation above.
Crucially, one of the key motivations for TMNRE is that one can actually
minimise this loss analytically to obtain the optimal classifier
$f_{\phi}^{\star}(x,\vartheta)\equiv\ln r(x;\vartheta)$ Miller  _et al._
(2021). This means that by successfully optimising this loss function, _we can
directly obtain full information about the likelihood-to-evidence ratio or
posterior-to-prior ratio._

In practice, this is where the “neural” aspect of TMNRE comes into play. Since
the data $x$ (and parameters $\vartheta$) can in principle be extremely high-
dimensional, we require both very flexible classifiers $d_{\phi}(x,\vartheta)$
that can process complex data structures such as images or long time series,
and also an efficient way to optimise the classifier parameters $\phi$ (here
neural network weights). Recent advances in machine learning architectures,
training, and hardware capabilities have now made this process extremely
scalable, and we make use of them here, as implemented in the swyft package
Miller  _et al._ (2021, 2022a).

It is important to understand how we can achieve marginal posterior estimates
with this algorithm, since it will be important for solving future sampling
challenges such as those associated to the analysis of multiple waveforms
Pizzati  _et al._ (2022); Antonelli  _et al._ (2021); Janquart  _et al._
(2022); Langendorff  _et al._ (2022); Samajdar  _et al._ (2021). In TMNRE, we
implicitly fully marginalise over every parameter
$\boldsymbol{\theta}=(\theta_{1},\theta_{2},\ldots)$ by varying them all in
each simulation. When we train the classifier(s), however, we only show the
network the small subset of parameters we are interested in, _i.e._
$(x,\vartheta=(\theta_{i},\theta_{j},\ldots))$. When one analytically
minimises the loss, one will find that you automatically obtain the properly
marginalised posterior $p(\vartheta|x)$. In practice, in this work, we will
only be interested in obtaining 1- and 2-dimensional marginal posteriors,
_i.e._ $\vartheta\equiv\theta_{k}$ for some $k$ or
$\vartheta\equiv(\theta_{i},\theta_{j})$ for some pair $(i,j)$.

We can finally turn to the inference and truncation part of TMNRE. With a
fully trained ratio estimator $\hat{r}(x;\vartheta)\simeq
p(\vartheta|x)/p(\vartheta)$, we can now target a specific observation $x_{0}$
(coming from _e.g._ an LVKC data stream, or some mock test observation). By
generating samples from the prior $p(\vartheta)$, we can obtain marginal
posterior inference results by weighting these samples with the ratio
estimator evaluated on the target observation $\hat{r}(x_{0};\vartheta)\simeq
p(\vartheta|x_{0})/p(\vartheta)$. Inevitably, this will lead to regions of
parameter space where the posterior density is high (covering the true values
if the algorithm is run successfully) and regions where the density is
vanishingly low. We use these regions of low density, specifically where the
posterior estimate drops below some threshold
$\hat{p}(\vartheta|x_{0})/\underset{\vartheta}{\mathrm{max}}\,\hat{p}(\vartheta|x_{0})<\epsilon$
to truncate our prior666This introduces a certain error into the marginal
posterior estimation proportional to the integrated quantity
$\int_{-\Gamma}\mathrm{d}\bar{\vartheta}\,p(\boldsymbol{\theta}|x_{0})$, where
$-\Gamma$ is the area outside the new proposal and $\bar{\vartheta}$ is the
set of parameters $\boldsymbol{\theta}$ not including the parameter(s) being
inferred. But, it is specifically in this region that the posterior density is
vanishingly low, so the error induced is small and controlled by $\epsilon$.
before generating samples for the next round of training and inference777We
typically take $\epsilon\sim 10^{-5}$ which approximately corresponds to a
$\sim 4.8\sigma$ exclusion if the posterior was a pure Gaussian. We checked
that the inference results are not sensitive to this choice provided
$\epsilon$ is small enough, only affecting the rate at which the truncation
occurs Miller  _et al._ (2022a). After a number of rounds, once the
sensitivity is at the level of statistical uncertainty, this truncation
procedure will converge, and the algorithm terminates. This process is
highlighted below for the specific case study presented in this work in Fig.
3.

To summarise, the TMNRE algorithm splits into four steps:

  1. 1.

Step 1: We generate a batch of simulations $(\vartheta,x)_{i=1...N}$ according
to our current proposal prior $p(\vartheta)$ using our forward
simulator.888Note that this step can be fully parallelised (and indeed is in
our implementation), and therefore made arbitrarily fast given hardware
access.

  2. 2.

Step 2: We train classifier(s)999In this work, we simultaneously train one
classifier for each gravitational wave parameter $\theta_{i}$, or pair
$(\theta_{i},\theta_{j})$, so we end up with a set of classifiers.
$f_{\phi}(x,\vartheta)$ to tell the difference between joint and marginal
samples for the parameter(s) of interest $\vartheta$. When optimised, this
gives a direct estimate of the ratio
$r(x;\vartheta)=p(\vartheta|x)/p(\vartheta)$ for the parameters of interest
$\vartheta=(\theta_{i},\theta_{j},\ldots)$.

  3. 3.

Step 3: We use the trained ratio estimator to perform inference on a specific
target observation $x_{0}$ and obtain marginal ratio estimates
$\hat{r}(x_{0};\vartheta)$.

  4. 4.

Step 4: By evaluating the ratios on a set of prior samples, we can truncate
the initial proposal distribution to exclude regions of vanishingly low
posterior density such that the variance in the training data for the next
round is significantly lower.

  5. 5.

We then go back to Step 1 and repeat until the truncation regions and
posterior estimates stabilise, after which we can perform the final inference
on the observation $x_{0}$.

This whole process is illustrated diagrammatically in Fig. 1.

###  II.3 Design Choices for Gravitational Waves

There are a number of concrete design choices and implementation steps that we
must make in order to apply the TMNRE algorithm to the analysis of
gravitational waves. Specifically we must _(i)_ choose/build the forward
simulator that generates the data $x$ given parameters $\theta$ and some noise
model, _(ii)_ implement a neural network that can parse and process the output
of the simulator for training the classifier, _(iii)_ choose the relevant
prior distributions for our parameters, and _(iv)_ make choices for the
relevant parameters in the TMNRE algorithm.

_Forward Simulator._ In terms of the generative model, in this work we choose
to study gravitational wave signals from compact binary black hole mergers
observed in LIGO and Virgo detectors. We analyse the signal as it appears in
all three detectors simultaneously101010Specifically, the Livingston (L1),
Hanford (H1) and Virgo (V1) detectors, although any combination or subset can
be analysed also., across both the time and frequency domains. To model the
waveforms from the merger, we use the functionality in the open-source code
bilby Ashton  _et al._ (2019); Romero-Shaw  _et al._ (2020); Ashton and Talbot
(2021). This allows us to inject a gravitational wave signal
$h(\boldsymbol{\theta}_{\mathrm{GW}})$ – typically split into the $+$ and
$\times$-polarisation strains and then projected – into a noisy detector with
coloured noise strain $n_{\mathrm{IFO}}$ generated from a power spectral
density (PSD) $S_{n}(f)$. Here,
$\boldsymbol{\theta}_{\mathrm{GW}}\equiv(q,\mathcal{M},\ldots)$ is some set of
parameters sampled from the prior $p(\boldsymbol{\theta}_{\mathrm{GW}})$ that
describe the intrinsic properties of the source under consideration such as
the mass ratio $q$, (redshifted) chirp mass of the system $\mathcal{M}$, as
well as extrinsic properties including the angular position on the sky
$(\alpha,\delta)$ and the luminosity distance $d_{L}$. All of these parameters
are provided as input into the waveform approximant of choice111111Here we use
the SEOBNRv4PHM waveform approximant Ossokine  _et al._ (2020) so as to
provide case studies on the most complex spinning and precessing BBH merger
systems. Importantly, the analysis pipeline remains identical regardless of
this choice, although the specific parameters in
$\boldsymbol{\theta}_{\mathrm{GW}}$ may change.. We provide a full list of
physical parameters in Tab. 1, including their prior choices. The full forward
model for a single detector is then defined by:

$p(x,\boldsymbol{\theta}_{\mathrm{GW}})=p(x=h+n_{\mathrm{IFO}}|h_{+,\times},\boldsymbol{\theta}_{\mathrm{GW}})\\\
\times
p(h_{+,\times}|\boldsymbol{\theta}_{\mathrm{GW}})p(\boldsymbol{\theta}_{\mathrm{GW}}),$
(4)

where $h$ here is the projection of the $h_{+,\times}$ strains onto the
detector frame according to the extrinsic parameters in
$\boldsymbol{\theta}_{\mathrm{GW}}$. The full implementation of the forward
model can be found in the Simulator class of the peregrine package. Two
example observations
$x_{0}=h(\boldsymbol{\theta}_{\mathrm{GW}})+n_{\mathrm{IFO}}$ which are the
direct output of the simulator and will form the basis of our case studies
below are shown in Fig. 2.

Parameter | Prior Choice | Injection
---|---|---
Mass ratio, $q$ | $\mathrm{U}(0.125,1)$ | 0.8858
Chirp mass $\mathcal{M}\,\mathrm{[M}_{\odot}\mathrm{]}$ | $\mathrm{U}(25,100)$ | 32.14
Inclination angle $\theta_{jn}\,\mathrm{[rad]}$ | $\mathrm{sine}(0,\pi)^{\dagger}$ | 0.4432
Polarisation angle $\psi\,\mathrm{[rad]}$ | $\mathrm{U}(0,\pi)$ | 1.100
Phase $\phi_{c}\,\mathrm{[rad]}$ | $\mathrm{U}(0,2\pi)$ | 5.089
Tilt angles $\theta_{1},\theta_{2}\,\mathrm{[rad]}$ | $\mathrm{sine}(0,\pi)^{\dagger}$ | 1.497, 1.102
Dimensionless spins $a_{1},a_{2}$ | $\mathrm{U}(0.05,1)$ | 0.9702, 0.8118
Spin angles $\phi_{12},\phi_{jl}\,\mathrm{[rad]}$ | $\mathrm{U}(0,2\pi)$ | 6.220, 1.885
Right ascension $\alpha\,\mathrm{[rad]}$ | $\mathrm{U}(0,2\pi)$ | 5.556
Declination $\delta\,\mathrm{[rad]}$ | $\mathrm{cosine}(-\pi/2,\pi/2)^{\dagger}$ | 0.071
Merger time $t_{c}\,\mathrm{[GPS\ s]}$ | $\mathrm{U}(-0.1,0.1)$ | 0.000
Luminosity Distance $d_{\mathrm{L}}\,\mathrm{[Mpc]}$ | $\mathrm{U}_{\mathrm{vol.}}(100,2000)^{\star}$ | 900.0 (C1), 200.0 (C2)

Table 1: Definitions and prior choices for all relevant gravitational wave
parameters $\boldsymbol{\theta}_{\mathrm{GW}}$ in this work.
† Note that these are the default priors used in $\mathrm{bilby}$ analyses for
BBH systems (subject to calibration ranges of the waveform approximants).
⋆ Specifically, the luminosity distance prior is taken to be uniform in
comoving volume in the _source frame_.

_Inference Network._ To carry out the ratio estimation step in the TMNRE
algorithm (Step 2 above), we need to construct a flexible enough
parameterisation of the classifier $f_{\phi}(x,\vartheta)$. In principle (at
least in some infinite training data limit), any sufficiently general
parameterisation should be able to minimise the loss in Eq. (3). In practice,
making sensible design choices about network architectures that are well
matched to the data format (particularly the signal structure) leads to huge
increases in performance, robustness and general applicability of the method.

In the context of gravitational waves, the signal is a complex time series or
frequency domain strain – _e.g._ at a sampling frequency of
$2048\,\mathrm{Hz}$ a $4\,\mathrm{s}$ signal will consist of over $8000$ data
points. Most of the information about the parameters
$\boldsymbol{\theta}_{\mathrm{GW}}$ is encoded in various non-linear ways. As
such, there is no good analytic method for fully processing noisy data
streams, or deriving an optimal architecture. On the other hand, in the
context of TMNRE, it is empirically the case that the precise network
architecture is not so important provided it is _(a)_ sensibly adjusted to the
data format – _e.g._ it would be non-optimal to first flatten an image,
totally removing all spatial structure – and _(b)_ deep enough (i.e., with
enough network parameters or structure) to be able to flexibly fit the
parameter and data manifold. This relative simplicity comes about as a result
of the specific targeting both of single observations and marginal posteriors,
significantly reducing the complexity of the fit we are implicitly trying to
perform. An interesting comparison is with the NPE implementation of
gravitational wave inference where a highly optimised and bespoke network
architecture is required to perform fully amortised inference on the high
dimensional joint posterior Dax  _et al._ (2021b).

Practically, we choose a network architecture that acknowledges a number of
things about the physics behind gravitational wave signals:

  * •

_Signal segments._ Broadly, compact binary coalescence gravitational wave
signals can be split into 3 components (both in the time and frequency
domains) corresponding to different phases of their evolution: the inspiral,
merger and ringdown phases. The various parameters in
$\boldsymbol{\theta}_{\mathrm{GW}}$ affect these components differently, for
example, the coalescence time $t_{c}$ defines the time of merger, typically
associated to the peak of the amplitude in the time domain signal $h(t)$. We
want to choose a network that does not treat each part of the waveform
identically – such as a convolutional structure that simply applies the same
kernel to each part. Instead we construct a 1d-analogue of the well-known unet
architecture Ronneberger  _et al._ (2015), which is excellent for image (or
signal here) segmentation and subsequent analysis on the various segments.

  * •

_Domain choice._ Information about gravitational wave parameters is encoded
differently in the time and frequency domains, although they are related
directly by a Fourier transform. At the end of the day, we are interested in
extracting some optimal set of features from the signal that are correlated
with various physical parameters. Naturally, the effect of some parameters is
much more direct in one domain vs. the other. For example, it is very simple
to see the effect of varying the time of coalescence on the time domain strain
$h(t)$ – it is simply a horizontal shift and is easy to pick up with any
simple network. On the other hand, in the frequency domain, the waveform gets
modified as $\tilde{h}(f;t_{c})=\tilde{h}(f;0)\exp(2\pi ift_{c})$, which
affects the full structure of the frequency domain strain and is much more
complicated to extract. In practice we take advantage of the full flexibility
of simulation-based methods and present _both_ the time and frequency domain
strains to the network, analysing both of them independently to extract
features before combining these summary statistics to perform ratio
estimation.

  * •

_Multiple detectors._ The current generation of detectors is set up to provide
complementary information about the same gravitational wave event at each of
the relevant facilities (_e.g._ LIGO Hanford and Virgo). We model the detector
response in all (or a subset) of these locations and stack the time (and
corresponding frequency) domain strains such that the same geocentric times
(or corresponding frequency bins) are aligned. Again, this is taking advantage
of the freedom to use any data representation that fully encodes our knowledge
about the relevant timing delays in the detectors.

The above describes the rationale for our choice of inference network used to
process our data. Its general structure consists of two parallel unet networks
which act on the time and frequency domain strains independently, before
compressing down to a smaller set of summary features in each case. These
features, which importantly are _automatically_ optimised by the classifier
during training are then concatenated and passed to a simple multi-layer
perceptron that is the default implemented in swyft which performs the ratio
estimation. More generally, one can understand the network as first performing
some optimised data compression into a set of summary statistics121212There
are some clear analogues with approximate Bayesian computation here Sunnaker
_et al._ (2013), although the summary statistics are trained to be optimal
rather than hand crafted. on which the ratio estimation is then performed with
(a subset of) the physical parameters
$\vartheta\subset\boldsymbol{\theta}_{\mathrm{GW}}$. The full details of the
network architecture can be found both in Appendix A as well as in the
peregrine implementation. In terms of generalisation, we expect this network
to perform well on any compact merger signal, including _e.g._ , overlapping
waveforms, although in principle it can also be modified arbitrarily to search
for other classes of gravitational waves without any other change in the data
analysis pipeline.

_Prior Choices._ The aim of this work is to establish and develop the TMNRE
algorithm in the context of gravitational wave analysis. As such, we choose
priors $p(\boldsymbol{\theta}_{\mathrm{GW}})$ on the physical parameters
$\boldsymbol{\theta}_{\mathrm{GW}}=(q,\mathcal{M},\cdots)$ that represent the
standard choices in the literature for analysing this type of BBH merger
signal Thrane and Talbot (2019); Veitch  _et al._ (2015); Ashton  _et al._
(2019); Biwer  _et al._ (2019). This allows for a direct comparison to
inference results from traditional algorithms that are the focus of the next
section. For reference, these can be found in _e.g._ the bilby documentation
Ashton  _et al._ (2019); Romero-Shaw  _et al._ (2020); Ashton and Talbot
(2021); Handley  _et al._ (2015); Speagle (2020). All the parameter
definitions and prior distributions relevant to this work can be found in Tab.
1.

TMNRE Setting | Value
---|---
Number of rounds | 7⋆
Simulation schedule | 30k, 60k, 90k, 120k, 120k,
| 150k, 150k
Bounds threshold $\epsilon$ | $10^{-5}$
Noise shuffling | True
Min./Max. training epochs | 30/200
Early stopping | 7
Initial learning rate | $5\times 10^{-4}$
Training/Validation batch size | 256/256
Train : Validation ratio | 0.9 : 0.1

Table 2: Choices for the hyperparameters and settings for the TMNRE algorithm
in this work.
⋆ This is the minimum number of rounds, if the algorithm has not converged, we
continue rounds of inference until the truncation procedure terminates.

_TMNRE Setup._ Finally, there are a number of choices that need to be made
when running the TMNRE algorithm. These choices are fully detailed in Tab. 2,
or in the configuration files supplied with peregrine. The settings broadly
split into two categories: those that control the training of the inference
network, and choices for the TMNRE algorithm itself. In the former case, we
can define how many epochs to train the network for (min./max. training
epochs), how many epochs to wait for the validation loss to
decrease131313Specifically, the training keeps track of the loss on the
training and validation sets, then tests whether the validation loss has
increased relative to the previous round. This can be a sign of overfitting,
which should of course be avoided, but we are also training via stochastic
gradient descent so we instead look for a number of epochs (early stopping)
where no decrease in the validation loss has been seen before stopping. We
then re-initialise the network to the state with the lowest validation loss.
(early stopping), the initial learning rate (initial learning rate), the batch
sizes shown to the network (training/validation batch size), and finally the
split in the data between training and validation (train : validation ratio).
In addition, we have the TMNRE settings such as the (minimum) number of rounds
(number of rounds), the number of simulations to use in each round141414In
later rounds, it is typically the case that a larger simulation budget is
required to obtain the necessary precision in the parameter inference.
(simulation schedule), and the bounds threshold for performing the truncation
($\epsilon$). The last setting to discuss is the noise shuffling switch. This
is closely related to the question of overfitting to training data, where for
small simulation batches it can be the case that the network essentially
“remembers” the training data including the noise realisations. This leads to
a very low value of the loss on the training data set, but poor performance on
validation data. We found that one extremely effective way of overcoming this
is to shuffle the noise realisations in each batch. This effectively shows the
network brand new examples every epoch, but with the same signal component
relevant for inference. This strategy should be applicable for _any_ additive
noise model (related to gravitational waves, or not) and can reduce the
simulation budget by potentially an order of magnitude.

To summarise this section, we have discussed broadly the field and
applicability of simulation-based inference techniques, before describing in
detail the development of our new gravitational wave data analysis pipeline
peregrine. In doing so, we explained the technicalities of our particular
implementation of simulation-based inference, TMNRE, that is applicable to
targeted parameter inference on a specific observation of interest. In the
next section, we will demonstrate the application of our analysis method to
two case studies of highly spinning, precessing BBH mergers, which represent a
current state-of-the-art analysis and modelling challenge.

##  III Results: Two Case Studies

Figure 3: Illustration of the truncation procedure as applied to the high SNR
case study (C2). Each row corresponds to one round of inference for five
separate parameters (chirp mass $\mathcal{M}$, right ascension $\alpha$,
inclination angle $\theta_{\mathrm{jn}}$, polarisation angle $\psi$, and
primary spin $a_{1}$, from left to right). The insets are zoomed in versions
of the constrained region to illustrate the coverage of the true injected
values (shown by black dashed lines). The solid blue lines show the prior
truncation bounds derived during the inference step.

We now turn our attention to the application of the method to two concrete
case studies that represent the state-of-the-art for parameter estimation in
current generation gravitational wave detectors. Specifically, in this section
we consider two noisy mock observations with waveforms generated by the
SEOBNRv4PHM approximant Ossokine  _et al._ (2020). This time-domain waveform
approximant is tuned to describe BBH merger systems that are both spinning and
precessing. It employs more accurate aligned-spin as well as precessing-spin
two-body dynamics than its predecessors.151515In particular, it includes the
$(2,\pm 2),(2,\pm 1),(3,\pm 3),(4,\pm 4)$ and $(5,\pm 5)$ modes in the co-
precessing frame of the binary instead of the commonly used $(2,\pm 2)$ and
$(2,\pm 1)$ modes, and displays improved accuracy in calibration fits with
numerical relativity simulations across a larger parameter space. When
generating aligned-spin waveforms, SEOBNRv4PHM reduces to SEOBNRv4HM. It takes
as input 15 relevant physical parameters, which are detailed in Tab. 1. We do
note, however, that our pipeline makes no assumption about the choice of
waveform approximant, so any other model could be easily integrated.

We generate the noise according to the LIGO-O4 power spectral density for each
detector (H1, L1, V1) Ashton and Talbot (2021). We assume the detector noise
to be gaussian-distributed, stationary and coloured by the estimated PSD of
the upcoming LIGO-Virgo observing run.161616In full technicality, the detector
noise originates from a number of different continuous and transient sources
which lead it to be non-stationary and coloured. However, it has been shown
that away from these transient noise artefacts a gaussian distribution
coloured by an online estimate of the PSD is an excellent approximation Abbott
_et al._ (2020). In practice, particularly because we are looking to analyse
individual observations, the noise model could be modified to include _e.g._ ,
any online estimate of the PSD at the time of an event Abbott  _et al._ (2016,
2019b). Whilst of course, the precise posteriors obtained would be different
as a result of the different noise realisation, the data analysis pipeline
would remain identical and we expect our approach to generalise to any noise
model.

In the rest of the section we will describe the case studies in more detail
before presenting our results when analysing the signals using the TMNRE
algorithm implemented in peregrine. Importantly, we show that we can achieve
excellent agreement compared to traditional established methods across all
parameters in the 1-dimensional and 2-dimensional marginals. To obtain this
comparison, we run the nested sampling code dynesty Speagle (2020) that comes
as default with the bilby package171717The results for the posteriors obtained
from dynesty Handley  _et al._ (2015); Speagle (2020) can be highly sensitive
to the sampling settings. In this work, we follow other analyses in the
literature Romero-Shaw  _et al._ (2020) taking nlive $=2000$ points as well as
n_act $=10$, and perform analytic phase marginalisation which is crucial for
stable results. We do not analytically marginalise over luminosity distance or
merger time, although this should not affect the results. Ashton  _et al._
(2019); Romero-Shaw  _et al._ (2020); Ashton and Talbot (2021). We also
discuss additional consistency checks that we can perform on simulation-based
inference methods that are independent of – and typically also infeasible to –
classical methods.

###  III.1 Case Study Description

We test our data analysis methodology on two separate case studies, both
within the context of the SEOBNRv4PHM waveform approximant Ossokine  _et al._
(2020), which is currently one of the most detailed and computationally
intensive models available (although, again, we could use any model, see e.g.
Khan  _et al._ (2016); Pratten  _et al._ (2020)). In particular we choose an
injected signal with the true values given in Tab. 1, which are chosen because
they are one set of best fit parameters for the interesting source recently
discussed in Refs. Hannam  _et al._ (2022); Payne  _et al._ (2022); Husa  _et
al._ (2016). In particular, the source comprises two objects that appear to be
both precessing and highly spinning (with dimensionless spins
$a_{i}>0.8$)181818In general, this interesting source highlights the need for
development in waveform modelling, as well as an investigation surrounding the
implications of noise artifacts on inference Romero-Shaw  _et al._ (2023); Xu
and Hamilton (2022). This is additional motivation to develop flexible,
targeted inference methods so that the impact of individual components of the
setup (waveform approximant, noise model etc.) on the final inference can be
reliably tested. This is very challenging to achieve with fully amortised
methods.. We place the systems at two different luminosity distances to set up
our two case studies:

  * •

C1. $d_{L}=900\,\mathrm{Mpc}$ – we study a quieter source (with an SNR of
$\sim$20) to ensure that our method can deal reliably with inference on a data
stream where there is significant amounts of noise. In other words, we test
the ability of our method to distinguish between signal and noise in a
gravitational wave event.

  * •

C2. $d_{L}=200\,\mathrm{Mpc}$ – conversely, we analyse a very loud source
(with an SNR of $\sim$100) to determine our ability to achieve high precision
inference. Whilst this source is arguably slightly unphysical according to
current merger rates (in terms of how loud it is), the aim of this test is to
confirm that we can achieve extremely high precision inference through the
application of the TMNRE algorithm.

The detector response for both of these target observations are shown in Fig.
2 above for the LIGO Hanford facility.

###  III.2 Parameter Estimation with TMNRE

We perform parameter inference on both of the case studies using the
implementation of the TMNRE algorithm in peregrine. In particular, we fix the
algorithm settings to the ones given in Tab. 2. The key results of this
section are shown in Figs. 4 and 5, while the truncation procedure is
illustrated for a number of parameters in Fig. 3. A full set of 2d-posterior
distributions can be found in the Appendix in Figs. 6 and 9.

_Discussion._ There are a number of levels to discuss these results. Firstly
we consider the TMNRE algorithm and its applicability to gravitational wave
data analysis. We can see clearly from the results in Figs. 4 and 5 that
across all parameters we have excellent reconstruction of the injected
parameters, including _e.g._ , the high primary spin which can be typically
challenging to analyse Romero-Shaw  _et al._ (2023); Xu and Hamilton (2022).
More importantly, we see that we achieve very close agreement with the
traditional nested sampling method dynesty Speagle (2020), both in terms of
accuracy of parameter reconstruction, as well as the precision of the
posteriors. To be quantitative, we compute the Jensen-Shannon divergence (JSD)
between the two sets of posterior distributions. In Appendix B, we provide a
full set of these values, averaged over all physical parameters across a
number of different samplers including dynesty Speagle (2020), but also cpnest
Veitch  _et al._ (2021) and ptemcee Foreman-Mackey  _et al._ (2013); Vousden
_et al._ (2016), as well as a complete description of our methodology for
computing the JSD191919We use nlive $=2000$, n_act $=10$ for dynesty, nlive
$=2000$ for cpnest and nsamples $=2000$ for ptemcee.. Given these results,
there are two key points to highlight here. Firstly, we see that the variation
between samplers as quantified by the JSD is an observation specific quantity
and that in this case the variance is significantly larger for the high SNR
case study (on average around $(10-30)\times 10^{-3}\,\mathrm{nat}$) compared
to the low SNR one (on average around $(3-9)\times 10^{-3}\,\mathrm{nat}$).
This trend follows at the level of variation between posterior distributions
obtained from the same sampler (e.g. dynesty) run twice with different random
seeds. Secondly, and more relevantly in terms of the peregrine posteriors, we
see that for both case studies, our results are fully consistent with the
variations between individual samplers (although of course we will happen to
show better agreement with one compared to another). In addition, it is worth
noting that we also obtain significantly better coverage statistics than
_e.g._ , bilby-MCMC and dynesty (see Refs. Ashton  _et al._ (2019); Ashton and
Talbot (2021) and the discussion below). As a final point, we can see that the
main discrepancies are on parameters that are nonetheless poorly constrained,
and certainly not measured to any precision – for example the secondary
dimensionless spin $a_{2}$, or the phase $\phi_{c}$ which is pre-marginalised
in the sampling approach anyway.

A crucial element to achieving this level of agreement was to ensure that the
same underlying noise realisation, and corresponding observation was analysed
in both the TMNRE and nested sampling analysis. Performing the parameter
inference on the same signal with a different noise realisation will
(correctly) lead to shifts in the posterior estimates of the parameter that
only _on average_ will be centred on the true value. This is something we
actually take advantage of in the next section to perform additional quality
checks of our inference. These are known broadly as coverage tests and are
becoming common practice in the simulation-based inference literature to
validate inference results Hermans  _et al._ (2021); Lemos  _et al._ (2023).

Figure 4: Full set of 1d marginal posteriors for all parameters in the low SNR
case study (C1). The inference results from the TMNRE algorithm are shown in
pink, with the $1\sigma$, $2\sigma$, and $3\sigma$ contours shown underlaid
behind. The purple histograms are the dynesty results for the analysis
performed on the same observation (see Fig. 2). The black dashed lines
indicate the true injected values. Figure 5: Full set of 1d marginal
posteriors for all parameters in the high SNR case study (C2). The inference
results from the TMNRE algorithm are shown in orange, with the $1\sigma$,
$2\sigma$, and $3\sigma$ contours shown underlaid behind. The purple
histograms are the dynesty results for the analysis performed on the same
observation (see Fig. 2). The black dashed lines indicate the true injected
values.

Two of the important features of the TMNRE algorithm that we highlighted in
Sec. II were the truncation and marginal elements. In Fig. 3, we highlight the
former of these and illustrate the impact of truncation on various classes of
parameters. For example, the first three columns of panels illustrate the
convergence of the posterior estimates over the rounds of truncation and re-
simulation for parameters that are (extremely) well measured, such as the
chirp mass or inclination angle. In the final two panels, we see another class
of parameters (here $\psi$ and $a_{1}$ as examples) which are not well
constrained but are nonetheless reconstructed. We see for this second class of
parameter that in the initial rounds, the marginal posterior estimate of
_e.g._ $p(\psi|x_{0})$ is rather poor – and indeed, one should not interpret
this as the true posterior, because the TMNRE algorithm has not converged.
Once the other well measured parameters are tightly constrained, however, the
quality of the inference of $\psi$ and $a_{1}$ increases significantly,
eventually converging to agree with the sampling result. This is a general
trend when running TMNRE algorithms that convergence in other parameters over
the course of the rounds ultimately leads to later improvements in posteriors
for different parameters, see _e.g._ Karchev  _et al._ (2022). It is also
worth briefly commenting on multi-modality here: whilst the case studies we
have chosen have well-constrained parameters such as the chirp mass that are
single modes, multi-modal inference in gravitational waves is a common
occurrence. In the context of peregrine, the same level of efficiency could be
realised in these scenarios by modifying how we perform the truncation of the
prior proposal. Specifically, instead of bounding in some rectangular region,
we could instead impose a threshold on the ratio and only sample in
(potentially disconnected) regions of the prior that exceed this. This is
currently under active development through _e.g._ , nested sampling applied to
the prior distributions, and we emphasise it is just technical hurdle, rather
than any limitation of the method. The second feature we highlighted was the
benefit of directly estimating marginal posteriors $p(\vartheta|x_{0})$ as
opposed to first generating the full ($N$-dimensional) joint distribution.
This plays a crucial role in simulation efficiency at no cost to statistical
accuracy – indeed the comparison with nested sampling highlights that we
obtain the correct marginal distributions even though we do not marginalise
from a joint distribution.202020Of course, as discussed in Sec. II, we may not
explicitly marginalise, but we do so implicitly by varying every parameter
across all simulations

We can explicitly highlight this simulation efficiency by comparing the number
of waveform generation steps that were required in the nested sampling case
versus in the TMNRE algorithm. In the low SNR case study, we required around
720k waveform evaluations, whilst the corresponding dynesty run required over
44 million. This represents an over 98% reduction in simulation budget for
identical marginal posterior reconstruction. The difference is actually larger
at the level of simulation time, since in each round of TMNRE, the waveform
calls can be fully parallelised, whilst evaluations in nested sampling or MCMC
broadly need to be performed sequentially. This parallelisation is implemented
as default in our peregrine code. This is not an entirely fair comparison,
since we are also required to train the classifiers at each stage. This should
be compared to _e.g._ , the total likelihood evaluation time for a meaningful
runtime comparison. Again, taking the low SNR case study as an example, the
dynesty run took around 43 hours versus around 12 hours for the TMNRE
algorithm on the same architecture (including network training).212121In
future code versions, there is room for this time to be substantially improved
through _e.g._ accelerated waveform evaluations, reduction in network size, or
a more efficient scheme for sampling from the truncated prior proposal
distributions. For scenarios where the cost per simulation is much higher,
however, this trade-off and benefit from simulation efficiency will become
even more pronounced.

As a final note, whilst at the level of parameter estimation the 1d marginal
posteriors are an informative set of results, it is often interesting to
explore the parameter degeneracies. To do so, we take advantage of the
flexibility of TMNRE to estimate only the marginal posteriors we are
interested in testing and construct the full set of 2d posteriors. We make a
comparison to the nested sampling contours in both case studies and find
similarly good agreement. To achieve this, we required an additional 500k
simulations in both the low and high SNR cases on which to train the set of
marginal posterior classifiers. These are shown in the Appendix in Figs. 6 and
9.

###  III.3 Additional Consistency Checks

The results and comparisons to traditional methods discussed above are a
helpful and meaningful way to test the data analysis pipeline that we have
developed. On the other hand, the key motivation for developing peregrine is
to apply it to scenarios where the application of traditional methods is
either impossible or extremely costly. This is slightly different motivation
to current research efforts to perform _e.g._ , online parameter inference for
LVKC data streams Dax  _et al._ (2021b); Wong  _et al._ (2023). Given this, it
is important to have additional ways to cross check and validate our inference
results that do not rely on comparisons to sampling. This is an active area of
ongoing research in simulation-based inference Hermans  _et al._ (2021);
Lueckmann  _et al._ (2021); Cannon  _et al._ (2022), however, a number of
established techniques now exist.

The most common are known as coverage tests, which implement the idea
discussed above that in the context of Bayesian parameter estimation, repeated
inference over different noise/statistical realisations of the same signal
should result in posteriors that shift relative to the true value. The
rationale behind this sort of expected coverage test is that – simply as a
result of statistical fluctuations – the $x$% credible interval for the
posterior should contain the simulation-truth value $x$% of the time. For a
more general discussion of this in the context of simulation-based inference,
see _e.g._ Hermans  _et al._ (2021). In practice, we can take a large
sample222222In the tests shown below, we generated 2000 extra simulations from
the final round truncated prior region. of observations generated from the
truncated prior (where we have converged posterior estimates), and perform
inference on each observation. In each case we can note how often the injected
value was contained inside the $x$% confidence interval and construct a
cumulative distribution. A well calibrated posterior distribution will be a
totally diagonal line when the expected coverage is plotted against the
empirical findings. Due to the fact that full inference must be done for each
observation, this type of test is typically infeasible to do for sequential
sampling methods such as MCMC. It should be noted, however, that this test is
diagnostic in the sense that a failure implies a poorly calibrated posterior,
but success is not a definitive guarantee.

The coverage test results for all 15 parameters in both the high and low SNR
case studies are provided in the Appendix (see Figs. 7, 8, 10 and 11). We see
that across every parameter, we achieve extremely good posterior coverage,
something that adds additional validation to our agreement with dynesty.

##  IV Conclusions and Outlook

In this final section, we present the key conclusions of this work, before
giving some outlook as to the possible use cases for our data analysis
strategy. Ultimately the motivation for our method development comes from the
analysis challenges facing the field of gravitational waves as we move towards
the next generation of detectors Maggiore  _et al._ (2020); Reitze  _et al._
(2019); Amaro-Seoane  _et al._ (2017). Specifically, we will need methods that
are highly scalable, simulation efficient, and flexible to take full advantage
of future data. In this context, the key research contributions from this work
are:

  * •

_Simulation-based inference pipeline._ The key development put forward in this
work is a sequential simulation-based inference pipeline for analysing
gravitational waves (see Sec. II for a discussion). Specifically, we
implemented an algorithm known as Truncated Marginal Neural Ratio Estimation
(TMNRE) Miller  _et al._ (2022a). The motivation and advantages associated to
this choice for gravitational waves are broadly centred around simulation
efficiency and its applicability to individual observations. We showed for
example that in a standard analysis for a binary black hole merger, we reduced
the simulation budget (number of waveform evaluations) by over 98% compared to
likelihood-based methods. This scaling is achieved primarily through the
ability of the TMNRE algorithm to directly estimate marginal posteriors, and
is the key argument for the use of our method in currently intractable
sampling problems such as the analysis of overlapping sources Pizzati  _et
al._ (2022); Antonelli  _et al._ (2021); Samajdar  _et al._ (2021); Janquart
_et al._ (2022); Langendorff  _et al._ (2022).

  * •

_Case studies._ We studied two cases where we analysed highly spinning
precessing binary black hole mergers, one with a high SNR of around 100, and
another with a relatively low SNR of approximately 20. In each case, we
performed full parameter inference using our method on observations generated
using the state-of-the-art SEOBNRv4PHM waveform approximant Ossokine  _et al._
(2020) (see Fig. 2), although the method can be applied to any waveform model.
We compared our results for the 1d- and 2d-marginal posteriors to posterior
samples obtained from nested sampling and demonstrated the excellent agreement
across all parameters (see Figs. 4 and 5). In addition, we carried out
expected coverage tests to validate the behaviour of our posterior estimates
(see _e.g._ Figs. 10 and 11 for the high SNR case).

  * •

_Robust and flexible method._ Our implementation of TMNRE is not the first
method to apply modern machine learning techniques to the analysis of
gravitational waves Gabbard  _et al._ (2022); Green and Gair (2021); Chua and
Vallisneri (2020); Álvares  _et al._ (2020); Dax  _et al._ (2021b); Delaunoy
_et al._ (2020); Dax  _et al._ (2022). On the other hand, these other
impressive parameter estimation pipelines typically have a different goal in
mind – rapid or instantaneous parameter inference. This has great relevance
when attempting to _e.g._ optimise signal triggering or follow up quickly on
electromagnetic counterparts. However, this comes at the cost of complicated
architectures, or bespoke likelihood designs Wong  _et al._ (2023); Zackay
_et al._ (2018); Leslie  _et al._ (2021) which render these approaches less
flexible for studying arbitrary classes of signals. In contrast, our
motivation is to develop a method that is robust and flexible enough to be
used to analyse a wide range of signal classes in specific observations. In
order to study a new class of signal – only a new forward simulation model is
required, whilst the rest of the pipeline will remain identical. Importantly,
in cases where traditional methods are not available, the coverage tests
described in Sec. III will still be applicable, so we can still validate our
results.

  * •

_Code release._ Alongside this work, we have developed an implementation of
our analysis method, known as peregrine (built on top of swyft Miller  _et
al._ (2021, 2022a)). This is a highly modular and scalable library that allows
the user to implement their own Simulator class. Given the simulation
efficiency and relatively small network size, it is usable on an individual
basis with reasonable access to consumer hardware. Both the algorithm as well
as the coverage tests are fully developed and tested, with examples available.
Access to peregrine can be requested at this link .

###  IV.1 Outlook

As discussed in the introduction to this work, we are currently in an exciting
era of gravitational wave experimental development. With future detectors such
as the Einstein Telescope Maggiore  _et al._ (2020) or LISA Amaro-Seoane  _et
al._ (2017) taking shape, the future gravitational wave sky could be extremely
loud, varied, and informative. On the other hand, there are clear data
analysis challenges associated to fully utilising this data Babak  _et al._
(2010); Schäfer  _et al._ (2023). Two important examples of these challenges
include _(i)_ the analysis of two (or more) overlapping gravitational wave
signals, which will arise as a result of large increases in detector
sensitivity Pizzati  _et al._ (2022); Antonelli  _et al._ (2021); Samajdar
_et al._ (2021); Janquart  _et al._ (2022); Langendorff  _et al._ (2022), and
_(ii)_ the identification and analysis of stochastic signals that will be
relevant to space-based observatories such as LISA Amaro-Seoane  _et al._
(2017). In the former case, the challenge is a sampling one, with current
methods potentially taking months to analyse a single event. In the latter
case, there is also potentially a large sampling challenge in an attempt to
fit multiple signal templates or noise models, but there are also statistical
challenges that could arise from correlated signal or detector noise, as well
as significant simulation costs.

It is this class of data analysis challenges that we had in mind when
developing peregrine.232323Although in addition, motivated by the observed
simulation efficiency, it would be interesting to explore the extent to which
the inference can be fully amortised by significantly expanding the initial
simulation budget, after which approximately real time inference could be
performed. We will leave this aspect to future work, however, along with
investigations as to how to include the effect of non-stationary noise in such
an approach. There are a number of reasons to believe that the framework put
forward here can approach these problems and extract the maximum
science/physics results. Firstly, the method is designed to be highly
simulation efficient, something we have demonstrated explicitly in the case of
binary black hole mergers where we saw a $\sim$98% reduction in simulation
budget. This is one of the keys to solving sampling based problems, where we
can directly perform inference at the level of marginal distributions rather
than solve the full joint problem. Assuming the same simulation efficiency
holds, this could lead to a reduction in analysis time for a problem such as
parameter inference for multiple overlapping waveforms from
$\mathcal{O}$(months) Janquart  _et al._ (2022) to less than a day. Similarly,
TMNRE and other simulation-based inference approaches are by construction
implicit-likelihood methods. This is equivalent to the statement that one only
needs a forward simulator to perform inference, rather than an explicit
likelihood. In the case of _e.g._ , possible correlated noise in the
stochastic gravitational wave background Caprini and Figueroa (2018);
Christensen (2019), this feature could unlock a number of physics questions
that are currently intractable. It is also worth noting that non-stationary
noise, which is typically challenging to deal with in a likelihood-based
framework is simple in the context of simulation-based inference and
peregrine. In particular, the analysis would remain identical if _e.g._ , the
PSD was time-dependent, provided this was included in the detector aspect of
the simulator. We will leave these to future works to analyse these situations
in detail.

To summarise, it is clear that in the era of high-sensitivity gravitational
wave experiments, flexible and scalable data analysis strategies will be
crucial to uncovering exciting fundamental physics discoveries. In this work,
we have presented peregrine as a community-based tool that can hopefully take
some steps towards that goal.

## Acknowledgements

This work is part of the project CORTEX (NWA.1160.18.316) of the research
programme NWA-ORC which is (partly) financed by the Dutch Research Council
(NWO). Additionally, CW acknowledges funding from the European Research
Council (ERC) under the European Union’s Horizon 2020 research and innovation
programme (Grant agreement No. 864035). UB is supported through the CORTEX
project of the NWA-ORC with project number NWA.1160.18.316 which is partly
financed by the Dutch Research Council (NWO). JA is supported through the
research program “The Hidden Universe of Weakly Interacting Particles” with
project number 680.92.18.03 (NWO Vrije Programma), which is partly financed by
the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (Dutch Research
Council). SN acknowledges support from the NWO Projectruimte grant (Samaya
Nissanke). BKM is funded by the University of Amsterdam Faculty of Science
(FNWI), Informatics Institute (IvI), and the Institute of Physics (IoP). BKM
is also affiliated with the European Laboratory for Learning and Intelligent
Systems (ELLIS Society). The main analysis for this work was carried out on
the Lisa and Snellius Computing Clusters at SURFsara. We are very grateful to
Thomas Edwards and Sam Witte for their careful reading of the draft and useful
comments. We also thank Andrew Williamson and Pablo Bosch for helpful
conversations about this project.

## References

  * Abbott  _et al._ (2016) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 116, 061102 (2016), arXiv:1602.03837 [gr-qc] .
  * Abbott  _et al._ (2021a) R. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. D 103, 122002 (2021a), arXiv:2010.14529 [gr-qc] .
  * Abbott  _et al._ (2021b) B. P. Abbott _et al._ (LIGO Scientific, Virgo, VIRGO), Astrophys. J. 909, 218 (2021b), arXiv:1908.06060 [astro-ph.CO] .
  * Nissanke  _et al._ (2010) S. Nissanke, D. E. Holz, S. A. Hughes, N. Dalal, and J. L. Sievers, Astrophys. J. 725, 496 (2010), arXiv:0904.1017 [astro-ph.CO] .
  * Nissanke  _et al._ (2011) S. Nissanke, J. Sievers, N. Dalal, and D. Holz, Astrophys. J. 739, 99 (2011), arXiv:1105.3184 [astro-ph.CO] .
  * Abbott  _et al._ (2021c) R. Abbott _et al._ (LIGO Scientific, Virgo), Astrophys. J. Lett. 913, L7 (2021c), arXiv:2010.14533 [astro-ph.HE] .
  * Abbott  _et al._ (2018a) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 121, 161101 (2018a), arXiv:1805.11581 [gr-qc] .
  * Abbott  _et al._ (2021d) R. Abbott _et al._ (LIGO Scientific, VIRGO), (2021d), arXiv:2108.01045 [gr-qc] .
  * Nitz  _et al._ (2021) A. H. Nitz, C. D. Capano, S. Kumar, Y.-F. Wang, S. Kastha, M. Schäfer, R. Dhurkunde, and M. Cabero, Astrophys. J. 922, 76 (2021), arXiv:2105.09151 [astro-ph.HE] .
  * Venumadhav  _et al._ (2019a) T. Venumadhav, B. Zackay, J. Roulet, L. Dai, and M. Zaldarriaga, Phys. Rev. D 100, 023011 (2019a).
  * Venumadhav  _et al._ (2020) T. Venumadhav, B. Zackay, J. Roulet, L. Dai, and M. Zaldarriaga, Phys. Rev. D 101, 083030 (2020).
  * Abbott  _et al._ (2021e) R. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. X 11, 021053 (2021e), arXiv:2010.14527 [gr-qc] .
  * Abbott  _et al._ (2021f) R. Abbott _et al._ (LIGO Scientific, VIRGO, KAGRA), (2021f), arXiv:2111.03606 [gr-qc] .
  * Allen  _et al._ (2012) B. Allen, W. G. Anderson, P. R. Brady, D. A. Brown, and J. D. E. Creighton, Phys. Rev. D 85, 122006 (2012), arXiv:gr-qc/0509116 .
  * Isi  _et al._ (2018) M. Isi, R. Smith, S. Vitale, T. J. Massinger, J. Kanner, and A. Vajpeyi, Phys. Rev. D 98, 042007 (2018), arXiv:1803.09783 [gr-qc] .
  * Venumadhav  _et al._ (2019b) T. Venumadhav, B. Zackay, J. Roulet, L. Dai, and M. Zaldarriaga, Phys. Rev. D 100, 023011 (2019b), arXiv:1902.10341 [astro-ph.IM] .
  * Abbott  _et al._ (2018b) B. P. Abbott _et al._ (KAGRA, LIGO Scientific, Virgo, VIRGO), Living Rev. Rel. 21, 3 (2018b), arXiv:1304.0670 [gr-qc] .
  * Petrov  _et al._ (2022) P. Petrov, L. P. Singer, M. W. Coughlin, V. Kumar, M. Almualla, S. Anand, M. Bulla, T. Dietrich, F. Foucart, and N. Guessoum, Astrophys. J. 924, 54 (2022), arXiv:2108.07277 [astro-ph.HE] .
  * Cutler and Flanagan (1994) C. Cutler and E. E. Flanagan, Phys. Rev. D 49, 2658 (1994), arXiv:gr-qc/9402014 .
  * Thrane and Talbot (2019) E. Thrane and C. Talbot, PASA 36, e010 (2019), arXiv:1809.02293 [astro-ph.IM] .
  * Ossokine  _et al._ (2020) S. Ossokine _et al._ , Phys. Rev. D 102, 044055 (2020), arXiv:2004.09442 [gr-qc] .
  * Khan  _et al._ (2016) S. Khan, S. Husa, M. Hannam, F. Ohme, M. Pürrer, X. Jiménez Forteza, and A. Bohé, Phys. Rev. D 93, 044007 (2016), arXiv:1508.07253 [gr-qc] .
  * Pratten  _et al._ (2020) G. Pratten, S. Husa, C. Garcia-Quiros, M. Colleoni, A. Ramos-Buades, H. Estelles, and R. Jaume, Phys. Rev. D 102, 064001 (2020), arXiv:2001.11412 [gr-qc] .
  * Edwards  _et al._ (2023) T. D. P. Edwards, K. W. K. Wong, K. K. H. Lam, A. Coogan, D. Foreman-Mackey, M. Isi, and A. Zimmerman, (2023), arXiv:2302.05329 [astro-ph.IM] .
  * Mackay (2003) D. Mackay,  _Informaiton Theory, Inference, and Learning Algorithms_ (Cambridge University Press, Cambridge, 2003).
  * Foreman-Mackey  _et al._ (2013) D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Goodman, Publ. Astron. Soc. Pac. 125, 306 (2013), arXiv:1202.3665 [astro-ph.IM] .
  * Skilling (2006) J. Skilling, Bayesian Analysis 1, 833 (2006).
  * Handley  _et al._ (2015) W. J. Handley, M. P. Hobson, and A. N. Lasenby, Mon. Not. Roy. Astron. Soc. 453, 4384 (2015), arXiv:1506.00171 [astro-ph.IM] .
  * Ashton  _et al._ (2022) G. Ashton _et al._ , Nature 2 (2022), 10.1038/s43586-022-00121-x, arXiv:2205.15570 [stat.CO] .
  * Speagle (2020) J. S. Speagle, Mon. Not. Roy. Astron. Soc. 493, 3132 (2020), arXiv:1904.02180 [astro-ph.IM] .
  * Veitch  _et al._ (2015) J. Veitch _et al._ , Phys. Rev. D 91, 042003 (2015), arXiv:1409.7215 [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, Publ. Astron. Soc. Pac. 131, 024503 (2019), arXiv:1807.10312 [astro-ph.IM] .
  * Ashton  _et al._ (2019) G. Ashton _et al._ , Astrophys. J. Suppl. 241, 27 (2019), arXiv:1811.02042 [astro-ph.IM] .
  * Romero-Shaw  _et al._ (2020) I. M. Romero-Shaw _et al._ , Mon. Not. Roy. Astron. Soc. 499, 3295 (2020), arXiv:2006.00714 [astro-ph.IM] .
  * Ashton and Talbot (2021) G. Ashton and C. Talbot, Mon. Not. Roy. Astron. Soc. 507, 2037 (2021), arXiv:2106.08730 [gr-qc] .
  * Williams  _et al._ (2023) M. J. Williams, J. Veitch, and C. Messenger, (2023), arXiv:2302.08526 [astro-ph.IM] .
  * Iacovelli  _et al._ (2022) F. Iacovelli, M. Mancarella, S. Foffa, and M. Maggiore, Astrophys. J. 941, 208 (2022), arXiv:2207.02771 [gr-qc] .
  * Pizzati  _et al._ (2022) E. Pizzati, S. Sachdev, A. Gupta, and B. Sathyaprakash, Phys. Rev. D 105, 104016 (2022), arXiv:2102.07692 [gr-qc] .
  * Antonelli  _et al._ (2021) A. Antonelli, O. Burke, and J. R. Gair, Mon. Not. Roy. Astron. Soc. 507, 5069 (2021), arXiv:2104.01897 [gr-qc] .
  * Samajdar  _et al._ (2021) A. Samajdar, J. Janquart, C. Van Den Broeck, and T. Dietrich, Phys. Rev. D 104, 044003 (2021), arXiv:2102.07544 [gr-qc] .
  * Janquart  _et al._ (2022) J. Janquart, T. Baka, A. Samajdar, T. Dietrich, and C. Van Den Broeck, (2022), arXiv:2211.01304 [gr-qc] .
  * Langendorff  _et al._ (2022) J. Langendorff, A. Kolmus, J. Janquart, and C. Van Den Broeck, (2022), arXiv:2211.15097 [gr-qc] .
  * Caprini and Figueroa (2018) C. Caprini and D. G. Figueroa, Class. Quant. Grav. 35, 163001 (2018), arXiv:1801.04268 [astro-ph.CO] .
  * Christensen (2019) N. Christensen, Rept. Prog. Phys. 82, 016903 (2019), arXiv:1811.08797 [gr-qc] .
  * Abbott  _et al._ (2019a) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. D 100, 061101 (2019a), arXiv:1903.02886 [gr-qc] .
  * Cranmer  _et al._ (2020) K. Cranmer, J. Brehmer, and G. Louppe, Proc. Nat. Acad. Sci. 117, 30055 (2020), arXiv:1911.01429 [stat.ML] .
  * Brehmer and Cranmer (2020) J. Brehmer and K. Cranmer, (2020), arXiv:2010.06439 [hep-ph] .
  * Lueckmann  _et al._ (2021) J.-M. Lueckmann, J. Boelts, D. S. Greenberg, P. J. Gonçalves, and J. H. Macke, (2021), arXiv:2101.04653 [stat.ML] .
  * Miller  _et al._ (2021) B. K. Miller, A. Cole, P. Forré, G. Louppe, and C. Weniger, in  _35th Conference on Neural Information Processing Systems_ (2021) arXiv:2107.01214 [stat.ML] .
  * Miller  _et al._ (2022a) B. K. Miller, A. Cole, C. Weniger, F. Nattino, O. Ku, and M. W. Grootes, J. Open Source Softw. 7, 4205 (2022a).
  * Tejero-Cantero  _et al._ (2020) A. Tejero-Cantero, J. Boelts, M. Deistler, J.-M. Lueckmann, C. Durkan, P. J. Gonçalves, D. S. Greenberg, and J. H. Macke, Journal of Open Source Software 5, 2505 (2020).
  * Alsing  _et al._ (2019) J. Alsing, T. Charnock, S. Feeney, and B. Wandelt, Mon. Not. Roy. Astron. Soc. 488, 4440 (2019), arXiv:1903.00007 [astro-ph.CO] .
  * Cole  _et al._ (2022) A. Cole, B. K. Miller, S. J. Witte, M. X. Cai, M. W. Grootes, F. Nattino, and C. Weniger, JCAP 09, 004 (2022), arXiv:2111.08030 [astro-ph.CO] .
  * Montel  _et al._ (2022) N. A. Montel, A. Coogan, C. Correa, K. Karchev, and C. Weniger, Mon. Not. Roy. Astron. Soc. 518, 2746 (2022), arXiv:2205.09126 [astro-ph.CO] .
  * Anau Montel and Weniger (2022) N. Anau Montel and C. Weniger, in  _36th Conference on Neural Information Processing Systems_ (2022) arXiv:2211.04291 [astro-ph.IM] .
  * Makinen  _et al._ (2021) T. L. Makinen, T. Charnock, J. Alsing, and B. D. Wandelt, JCAP 11, 049 (2021), arXiv:2107.07405 [astro-ph.CO] .
  * Dimitriou  _et al._ (2022) A. Dimitriou, C. Weniger, and C. A. Correa, (2022), arXiv:2206.11312 [astro-ph.CO] .
  * Gagnon-Hartman  _et al._ (2023) S. Gagnon-Hartman, J. Ruan, and D. Haggard, Mon. Not. Roy. Astron. Soc. 520, 1 (2023), arXiv:2301.05241 [astro-ph.CO] .
  * Delaunoy  _et al._ (2020) A. Delaunoy, A. Wehenkel, T. Hinderer, S. Nissanke, C. Weniger, A. R. Williamson, and G. Louppe, (2020), arXiv:2010.12931 [astro-ph.IM] .
  * Karchev  _et al._ (2022) K. Karchev, R. Trotta, and C. Weniger,  (2022), 10.1093/mnras/stac3785, arXiv:2209.06733 [astro-ph.CO] .
  * Lin  _et al._ (2022) K. Lin, M. von Wietersheim-Kramsta, B. Joachimi, and S. Feeney, (2022), arXiv:2212.04521 [astro-ph.CO] .
  * Dax  _et al._ (2021a) M. Dax, S. R. Green, J. Gair, M. Deistler, B. Schölkopf, and J. H. Macke, (2021a), arXiv:2111.13139 [cs.LG] .
  * Wildberger  _et al._ (2022) J. Wildberger, M. Dax, S. R. Green, J. Gair, M. Pürrer, J. H. Macke, A. Buonanno, and B. Schölkopf, (2022), arXiv:2211.08801 [gr-qc] .
  * Dax  _et al._ (2021b) M. Dax, S. R. Green, J. Gair, J. H. Macke, A. Buonanno, and B. Schölkopf, Phys. Rev. Lett. 127, 241103 (2021b), arXiv:2106.12594 [gr-qc] .
  * Pratten  _et al._ (2021) G. Pratten _et al._ , Phys. Rev. D 103, 104056 (2021), arXiv:2004.06503 [gr-qc] .
  * Hannam  _et al._ (2014) M. Hannam, P. Schmidt, A. Bohé, L. Haegel, S. Husa, F. Ohme, G. Pratten, and M. Pürrer, Phys. Rev. Lett. 113, 151101 (2014), arXiv:1308.3271 [gr-qc] .
  * Babak  _et al._ (2010) S. Babak _et al._ (Mock LISA Data Challenge Task Force), Class. Quant. Grav. 27, 084009 (2010), arXiv:0912.0548 [gr-qc] .
  * Schäfer  _et al._ (2023) M. B. Schäfer _et al._ , Phys. Rev. D 107, 023021 (2023), arXiv:2209.11146 [astro-ph.IM] .
  * Papamakarios and Murray (2016) G. Papamakarios and I. Murray, (2016), arXiv:1605.06376 [stat.ML] .
  * Zeghal  _et al._ (2022) J. Zeghal, F. Lanusse, A. Boucaud, B. Remy, and E. Aubourg, (2022), arXiv:2207.05636 [astro-ph.IM] .
  * Dax  _et al._ (2022) M. Dax, S. R. Green, J. Gair, M. Pürrer, J. Wildberger, J. H. Macke, A. Buonanno, and B. Schölkopf, (2022), arXiv:2210.05686 [gr-qc] .
  * Rozet and Louppe (2021) F. Rozet and G. Louppe, (2021), arXiv:2110.00449 [cs.LG] .
  * Delaunoy  _et al._ (2022) A. Delaunoy, J. Hermans, F. Rozet, A. Wehenkel, and G. Louppe, (2022), arXiv:2208.13624 [stat.ML] .
  * Miller  _et al._ (2022b) B. K. Miller, C. Weniger, and P. Forré, (2022b), arXiv:2210.06170 [stat.ML] .
  * Hermans  _et al._ (2019) J. Hermans, V. Begy, and G. Louppe, (2019), arXiv:1903.04057 [stat.ML] .
  * Durkan  _et al._ (2020) C. Durkan, I. Murray, and G. Papamakarios, (2020), arXiv:2002.03712 [stat.ML] .
  * Ronneberger  _et al._ (2015) O. Ronneberger, P. Fischer, and T. Brox, (2015), arXiv:1505.04597 [cs.CV] .
  * Sunnaker  _et al._ (2013) M. Sunnaker, A. Busetto, E. Numminen, J. Corander, M. Foll, and C. Dessimoz, PLoS Comput Biol 9, e1002803 (2013).
  * Abbott  _et al._ (2020) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Class. Quant. Grav. 37, 055002 (2020), arXiv:1908.11170 [gr-qc] .
  * Abbott  _et al._ (2019b) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. X 9, 031040 (2019b), arXiv:1811.12907 [astro-ph.HE] .
  * Hannam  _et al._ (2022) M. Hannam _et al._ , Nature 610, 652 (2022), arXiv:2112.11300 [gr-qc] .
  * Payne  _et al._ (2022) E. Payne, S. Hourihane, J. Golomb, R. Udall, R. Udall, D. Davis, and K. Chatziioannou, Phys. Rev. D 106, 104017 (2022), arXiv:2206.11932 [gr-qc] .
  * Husa  _et al._ (2016) S. Husa, S. Khan, M. Hannam, M. Pürrer, F. Ohme, X. J. Forteza, and A. Bohé, Phys. Rev. D 93, 044006 (2016).
  * Romero-Shaw  _et al._ (2023) I. M. Romero-Shaw, D. Gerosa, and N. Loutrel, Mon. Not. Roy. Astron. Soc. 519, 5352 (2023), arXiv:2211.07528 [astro-ph.HE] .
  * Xu and Hamilton (2022) Y. Xu and E. Hamilton, (2022), arXiv:2211.09561 [gr-qc] .
  * Veitch  _et al._ (2021) J. Veitch, W. D. Pozzo, A. Lyttle, M. Williams, C. Talbot, M. Pitkin, G. Ashton, Cody, M. Hübner, A. Nitz, D. Macleod, G. Carullo, G. Davies, and Tony, “cpnest: v0.11.3,”  (2021).
  * Vousden  _et al._ (2016) W. D. Vousden, W. M. Farr, and I. Mandel, Mon. Not. Roy. Astron. Soc. 455, 1919 (2016), arXiv:1501.05823 .
  * Hermans  _et al._ (2021) J. Hermans, A. Delaunoy, F. Rozet, A. Wehenkel, V. Begy, and G. Louppe, (2021), arXiv:2110.06581 [stat.ML] .
  * Lemos  _et al._ (2023) P. Lemos, A. Coogan, Y. Hezaveh, and L. Perreault-Levasseur, (2023), arXiv:2302.03026 [stat.ML] .
  * Wong  _et al._ (2023) K. W. K. Wong, M. Isi, and T. D. P. Edwards, (2023), arXiv:2302.05333 [astro-ph.IM] .
  * Cannon  _et al._ (2022) P. Cannon, D. Ward, and S. M. Schmon, (2022), arXiv:2209.01845 [stat.ML] .
  * Maggiore  _et al._ (2020) M. Maggiore _et al._ , JCAP 03, 050 (2020), arXiv:1912.02622 [astro-ph.CO] .
  * Reitze  _et al._ (2019) D. Reitze _et al._ , Bull. Am. Astron. Soc. 51, 035 (2019), arXiv:1907.04833 [astro-ph.IM] .
  * Amaro-Seoane  _et al._ (2017) P. Amaro-Seoane, H. Audley, S. Babak, J. Baker, E. Barausse, P. Bender, E. Berti, _et al._ , (2017), arXiv:1702.00786 [astro-ph.IM] .
  * Gabbard  _et al._ (2022) H. Gabbard, C. Messenger, I. S. Heng, F. Tonolini, and R. Murray-Smith, Nature Phys. 18, 112 (2022), arXiv:1909.06296 [astro-ph.IM] .
  * Green and Gair (2021) S. R. Green and J. Gair, Mach. Learn. Sci. Tech. 2, 03LT01 (2021), arXiv:2008.03312 [astro-ph.IM] .
  * Chua and Vallisneri (2020) A. J. K. Chua and M. Vallisneri, Phys. Rev. Lett. 124, 041102 (2020), arXiv:1909.05966 [gr-qc] .
  * Álvares  _et al._ (2020) J. a. D. Álvares, J. A. Font, F. F. Freitas, O. G. Freitas, A. P. Morais, S. Nunes, A. Onofre, and A. Torres-Forné (2020) arXiv:2011.10425 [gr-qc] .
  * Zackay  _et al._ (2018) B. Zackay, L. Dai, and T. Venumadhav, (2018), arXiv:1806.08792 [astro-ph.IM] .
  * Leslie  _et al._ (2021) N. Leslie, L. Dai, and G. Pratten, Phys. Rev. D 104, 123030 (2021), arXiv:2109.09872 [astro-ph.IM] .
  * Paszke  _et al._ (2019) A. Paszke, S. Gross, F. Massa, _et al._ , (2019), arXiv:1912.01703 [cs.LG] .

##  Appendix A Network Architecture

The network architecture used for the analysis in peregrine is shown in the
figure below. As described in the text in Sec. II, we use the flexibility of
the TMNRE implementation to automatically extract a small set of summary
statistics from both the time and frequency domain strains after processing
and compressing the input data. This summary can then be passed to the ratio
estimator along with the physical parameters
$\boldsymbol{\theta}_{\mathrm{GW}}$ to perform the inference. The full pytorch
Paszke  _et al._ (2019) implementation of the network components shown below
(specifically the unet and linear compression networks) can be found in the
InferenceNetwork class of peregrine.

##  Appendix B Comparison with Likelihood Based Samplers

In order to validate our results in Sec. III, we quantitatively compare the
peregrine posteriors against those obtained using traditional sampling methods
such as nested sampling and MCMC. Specifically, we compute the Jensen-Shannon
divergence amongst the 1d posteriors obtained using peregrine, dynesty,
ptemcee, and cpnest. Our implementation including all the data relevant to
this article are available publicly on our Zenodo page (link.). section of the
peregrine repository, along with files that can regenerate the JSD values
(both on average and on a parameter-by-parameter basis) given below in Tab. 3.

$\mathbf{JSD}\,\mathrm{[}10^{-3}\,\mathrm{nat]}$ | peregrine | dynesty |  dynesty (re-run) | ptemcee | cpnest
---|---|---|---|---|---
peregrine |  | $18.0$ | $13.1$ | $21.0$ | $20.6$
dynesty | $9.75$ |  | $4.02$ | $10.0$ | $35.2$
dynesty (re-run) | $10.9$ | $0.55$ |  | $10.6$ | $29.9$
ptemcee | $16.5$ | $7.01$ | $7.56$ |  | $35.9$
cpnest | $6.82$ | $3.52$ | $4.46$ | $8.60$ |

Table 3: Average Jensen-Shannon divergences over all physical parameters
$\boldsymbol{\theta}_{\mathrm{GW}}=(q,\mathcal{M},\ldots)$ for different
sampler combinations. Values below the diagonal (in pink boxes) refer to the
low SNR case study (C1), whilst those above the diagonal (in orange boxes) are
for the high SNR case (C2).

Figure 6: Full set of 2d marginal posteriors in the low SNR case study (C1).
The results from the TMNRE analysis are shown in pink, whilst the
corresponding dynesty analysis is shown in blue. The darker and subsequently
lighter contours in the 2d marginals indicate the $1\sigma$, $2\sigma$ and
$3\sigma$ confidence intervals respectively. The sky map shows the
($\alpha$,$\delta$) contours centered at the injection value. Figure 7:
Coverage results for the low SNR case study (C1) for all parameters. This
compares the expected coverage of the true value as a percentage on the
$x$-axis against the actual coverage of our ratio estimator on the $y$-axis.
The pink line indicates the average coverage, while the blue contour indicates
the $68$% confidence interval on the coverage. Figure 8: Coverage results for
the low SNR case study (C1) for all parameters. This is the same information
as Fig. 10, but with $z_{p}$ defined by
$p=\int_{-z_{p}}^{z_{p}}{\mathrm{d}z\,1/\sqrt{2\pi}\exp(-z^{2}/2)}$. This
places more emphasis on the behaviour of the posteriors in the tail regions.
The pink line indicates the average coverage, while the blue contour indicates
the $68$% confidence interval on the coverage. Figure 9: Full set of 2d
marginal posteriors in the high SNR case study (C2). The results from the
TMNRE analysis are shown in orange, whilst the corresponding dynesty analysis
is shown in blue. The darker and subsequently lighter contours in the 2d
marginals indicate the $1\sigma$, $2\sigma$ and $3\sigma$ confidence intervals
respectively. The sky map shows the ($\alpha$,$\delta$) contours centered at
the injection value. Figure 10: Coverage results for the high SNR case study
(C2) for all parameters. This compares the expected coverage of the true value
as a percentage on the $x$-axis against the actual coverage of our ratio
estimator on the $y$-axis. The pink line indicates the average coverage, while
the blue contour indicates the $68$% confidence interval on the coverage.
Figure 11: Coverage results for the high SNR case study (C2) for all
parameters. This is the same information as Fig. 10, but with $z_{p}$ defined
by $p=\int_{-z_{p}}^{z_{p}}{\mathrm{d}z\,1/\sqrt{2\pi}\exp(-z^{2}/2)}$. This
places more emphasis on the behaviour of the posteriors in the tail regions.
The pink line indicates the average coverage, while the blue contour indicates
the $68$% confidence interval on the coverage.