#  Eryn††thanks: https://github.com/mikekatz04/Eryn : A multi-purpose sampler
for Bayesian inference

Nikolaos Karnesis,1 Michael L. Katz,2 Natalia Korsakova,3 Jonathan R. Gair,2
Nikolaos Stergioulas1
1 Department of Physics, Aristotle University of Thessaloniki, Thessaloniki
54124, Greece
2Max-Planck-Institut für Gravitationsphysik, Albert-Einstein-Institut, Am
Mühlenberg 1, 14476 Potsdam-Golm, Germany
3Astroparticule et Cosmologie, CNRS, Université Paris Cité, F-75013 Paris,
France E-mail<EMAIL_ADDRESS>
###### Abstract

In recent years, methods for Bayesian inference have been widely used in many
different problems in physics where detection and characterization are
necessary. Data analysis in gravitational-wave astronomy is a prime example of
such a case. Bayesian inference has been very successful because this
technique provides a representation of the parameters as a posterior
probability distribution, with uncertainties informed by the precision of the
experimental measurements. During the last couple of decades, many specific
advances have been proposed and employed in order to solve a large variety of
different problems. In this work, we present a Markov Chain Monte Carlo (MCMC)
algorithm that integrates many of those concepts into a single MCMC package.
For this purpose, we have built Eryn, a user-friendly and multipurpose toolbox
for Bayesian inference, which can be utilized for solving parameter estimation
and model selection problems, ranging from simple inference questions, to
those with large-scale model variation requiring trans-dimensional MCMC
methods, like the LISA global fit problem. In this paper, we describe this
sampler package and illustrate its capabilities on a variety of use cases.

###### keywords:

methods: data analysis – software: data analysis – software: development –
gravitational waves.

††pubyear: 2023

##  1 Introduction

In physics, and in science in general, one of the most encountered problems is
the one of model calibration and comparison. We test our models of the
physical world against the measured data, to estimate their parameters and to
robustly determine the most suitable model that describes our observations. A
crucial first step in this direction, is to efficiently explore the posterior
distribution of the parameters given the measured data. Markov Chain Monte
Carlo (MCMC) algorithms have proven to be very successful in this regard
(Metropolis et al., 1953; Hastings, 1970; Hitchcock, 2003; Gilks et al.,
1995), being one of the few methods which can efficiently perform Bayesian
inference when the posterior is not analytically tractable and without solving
exactly for the marginal likelihood. This is compared to, for example, grid
methods, which are often computationally unfeasible. This is especially true
in the field of Gravitational-Wave (GW) astronomy, where MCMC methods have
been extensively used in order to find physical parameters for signals buried
in the data (see e.g. Ashton & Talbot (2021); Biwer et al. (2019); Abbott et
al. (2021a); Collaboration & the Virgo Collaboration (2020); Shawhan (2003);
Abbott et al. (2021b)), as well as to hierarchically infer the properties of
the underlying astrophysical populations (e.g. Ashton & Prix (2018); Thrane &
Talbot (2019); Isi et al. (2022)). MCMC approaches can also compute the
marginal likelihood or evidence (see section 2), by using techniques such as
Thermodynamic Integration (see section 2.4). In a Bayesian framework, the
evidence difference between two models can be used to compute the Bayes
Factor, which is used to select between different models that could describe
the observations.

Thermodynamic integration (or other approximations, see section 2.5 or Cornish
& Littenberg (2007); Gelman et al. (1996)), is ideal for cases where the
number of competing models is small. However, in situations where the number
of potential models becomes too large, the task of iteratively and
hierarchically computing the marginal likelihood can become computationally
inefficient, or even practically unachievable. Such is the case for future
signal-dominated GW observatories, such as the Laser Interferometer Space
Antenna (LISA) (Amaro-Seoane et al., 2017) or other proposed space-borne GW
observatories (Luo et al., 2016; Kawamura et al., 2008; Ren et al., 2023).
LISA will observe different types of GW sources, the most numerous of them
being the Ultra Compact Binaries (UCBs) within the Milky Way (Amaro-Seoane et
al., 2017, 2012; Littenberg et al., 2020; Littenberg & Cornish, 2023; Strub et
al., 2022; Zhang et al., 2021; Karnesis et al., 2021; Crowder & Cornish,
2007). Those are mostly comprised of a population of Double White Dwarfs
(DWD), although a small fraction of neutron star - white dwarf (NS-WD) or
double neutron stars (NS-NS) binaries are expected (Breivik et al., 2020;
Nelemans et al., 2001). In fact, LISA is going to detect GWs from the complete
population of $\mathcal{O}(10^{7})$ sources simultaneously, but only a small
fraction of them are going to be individually resolvable
($\mathcal{O}(10^{4})$). The large majority of signals will generate an
anisotropic and non-stationary “confusion” type of signal, which will dominate
the LISA band between $0.05$ and $\sim 0.2$ $\mathrm{mHz}$.

In the above context, computing the marginal likelihood for such a large
parameter space and for all possible numbers of events that could be in the
population becomes computationally prohibitive. Instead, we can employ
dynamical trans-dimensional MCMC methods (Green, 1995). This family of methods
can be quite challenging to tune, but it has proven to yield satisfactory
results, even for such demanding problems as the LISA data (Littenberg et al.,
2020; Littenberg & Cornish, 2023). There are also implementation challenges,
which arise from technical aspects of the algorithm; one example being the
dimension matching requirement when proposing moves between models with
different dimensionality. In terms of algorithm efficiency, it is also crucial
to choose proposal distributions that allow smooth transitions on a dynamical
parameter space, a task which in many cases requires substantial effort. For
these technical reasons, all the available software tools of this kind have
been specifically developed for the particular problem they intend to solve.

In this work, we present Eryn, a reversible jump MCMC algorithm, capable of
efficiently sampling dynamical parameter spaces, while remaining generic and
usable by a large community. We build upon various ideas from statistics,
astronomy, etc., in order to develop an efficient statistical toolbox that can
be applied to the majority of problems involving detection and
characterization of signals. Our primary goal, however, is to utilize Eryn as
a basic ingredient for a data analysis pipeline to perform the LISA Global Fit
(Littenberg et al., 2020; Littenberg & Cornish, 2023). The Global Fit is a
data analysis strategy required to tackle the problem of multiple source
detection, separation, and characterization in LISA data. For demonstration
purposes in this work, we use Eryn to analyze a “reduced” scenario of the LISA
data in section 4.

This paper is organized as follows: in section 2, we begin with explaining the
foundations of the MCMC algorithm, as well as some of the relevant methods
that we have adopted for our implementation. In section 3, we describe how the
methods introduced in section 2 are combined into the actual toolbox
implemented in Eryn. In section 3.3, we demonstrate the capabilities of Eryn
through some toy examples, while in section 4 we apply our machinery to more
demanding applications in Gravitational-Wave astronomy. In particular, we
demonstrate Eryn by performing model selection on a simulated population of
Ultra Compact Galactic Binaries (UCBs) as measured by the future LISA
observatory. Finally, in section 5, we summarize our work and discuss future
applications. We should state again here, that Eryn is available as open
source software in https://github.com/mikekatz04/Eryn.

##  2 Markov Chain Monte Carlo algorithms

Nowadays, MCMC methods are considered to be a cornerstone of Bayesian
Inference, being very effective in finding solutions to problems encountered
across wide-ranging disciplines (e.g. Hobson et al., 2009; Ashton & Talbot,
2021; Biwer et al., 2019; Hogg & Foreman-Mackey, 2018; Sharma, 2017; Johannes
& Polson, 2009). These include the sampling of the posterior densities of
parameters of interest, the numerical marginalisation over nuisance
parameters, and providing a framework to compute the marginal posterior
distributions (or evidences) that can be used for model selection. The
Bayesian framework is based around Bayes’ Theorem:

$p(\vec{\theta}|y,\mathcal{M})=\frac{p(y|\vec{\theta},\mathcal{M})p(\vec{\theta}|\mathcal{M})}{p(y|\mathcal{M})},$
(1)

where $y$ is the measured data and $\mathcal{M}$ our chosen model of analysis.
The $p(\vec{\theta}|y)$ term is the posterior distribution of the parameter
set $\vec{\theta}$, which is related to the likelihood function of the data
$p(y|\vec{\theta},\mathcal{M})$ and the prior densities of the parameters
$p(\vec{\theta}|\mathcal{M})$. The evidence $p(y|\mathcal{M})$ is the marginal
posterior over the parameter space $\vec{\theta}\in\Theta$:

$\displaystyle\mathcal{Z}\equiv p(\vec{y}|\mathcal{M})$
$\displaystyle=\int_{\Theta}p(\vec{\theta},\vec{y}|\mathcal{M})\mathrm{d}\vec{\theta}$
$\displaystyle=\int_{\Theta}p(\vec{y}|\vec{\theta},\mathcal{M})p(\vec{\theta}|\mathcal{M})\mathrm{d}\vec{\theta}.$
(2)

For parameter estimation purposes, the evidence acts as a normalization
constant and can be ignored. However, it is really important if one wants to
perform model selection over the measured data. We shall describe in detail
how one can numerically approximate the integral of eq. (2) in section 2.5.

MCMC algorithms work by constructing a Markov Chain sequence, whose elements,
$\vec{\theta}(t_{i})$, for $i=0,1,\ldots$, are (asymptotically) independent
samples from the target distribution, $f(\vec{\theta})$. Under fairly general
assumptions, the distribution of samples in the chain will converge to the
target distribution provided the algorithm satisfies detailed balance:

$f(\vec{\theta})p(\vec{\theta}\rightarrow\vec{\theta}^{\prime})=f(\vec{\theta}^{\prime})p(\vec{\theta}^{\prime}\rightarrow\vec{\theta}).$
(3)

Here $p(\vec{\theta}\rightarrow\vec{\theta}^{\prime})$ is the probability that
the Markov chain moves from point $\vec{\theta}$ to point
$\vec{\theta}^{\prime}$. The most widely-used MCMC algorithm is Metropolis-
Hastings (Metropolis et al., 1953; Hastings, 1970), which is explained in
algorithm box 1. The first step of the algorithm is to define an initial
state, $\vec{\theta}(t_{0})$. Then, at each subsequent step $i$, a new state
is proposed, by randomly drawing from a given proposal distribution
$q(\vec{\theta}^{\prime}|\vec{\theta}(t_{i}))$. The newly proposed state is
then accepted with a certain probability, given by eq. (4). If the move is
accepted we set $\vec{\theta}(t_{i+1})=\vec{\theta}^{\prime}$, otherwise we
set $\vec{\theta}(t_{i+1})=\vec{\theta}(t_{i})$. Any reasonable choice of the
proposal density will generate a Markov chain with the correct stationary
distribution. However, a good choice of $q$ is critical for its efficiency,
i.e. achieving the convergence of the MCMC chains within a reasonable
computational time. For the special case of a symmetric proposal
distribution111That is, a proposal with the property
$p(\vec{\theta}\rightarrow\vec{\theta}^{\prime})=p(\vec{\theta}^{\prime}\rightarrow\vec{\theta})$.,
such as the widely used multivariate Gaussian distribution centered around
$\vec{\theta}(t_{i})$, the ratio of eq. (4) in algorithm box 1 becomes simply
the ratio of the target densities at the current $\vec{\theta}(t_{n})$ and
proposed $\vec{\theta}^{\prime}$ points. For high-dimensional problems, the
multivariate Gaussian proposal can be tuned during the burn-in period of
sampling to improve efficiency (Haario et al., 2001; Frenkel, 2006; Andrieu &
Thoms, 2008; Gelman et al., 1996; Roberts & Rosenthal, 2009; Christensen &
Meyer, 2022), or even scaled according to the Cramer-Rao bound, estimated from
the Information matrix (Vallisneri, 2008).

1:Set the initial state of chain: $\vec{\theta}(t_{0})\equiv\vec{\theta}$.

2:At each subsequent step $n$, draw a new state $\vec{\theta}^{\prime}$ from a
proposal distribution $q(\vec{\theta}^{\prime}|\vec{\theta}(t_{n-1}))$.

3:Take
$\vec{\theta}(t_{n})=\left\\{\begin{array}[]{lr}\vec{\theta}^{\prime}\mbox{
with probability }\alpha(\vec{\theta}(t_{n-1}),\vec{\theta}^{\prime})\\\
\vec{\theta}\mbox{ with probability
}1-\alpha(\vec{\theta}(t_{n-1}),\vec{\theta}^{\prime})\end{array}\right.,$

4:where  $\alpha(x,y)=1\wedge\left\\{\frac{f(y)q(x|y)}{f(x)q(y|x)}\right\\}.$
(4)

5:Go to 2 and repeat until equilibrium is reached, and enough independent
samples have been drawn from the target distribution.

Algorithm 1 The Metropolis-Hastings algorithm (Metropolis et al., 1953;
Hastings, 1970). $f$ is the target density to be sampled.

Although the MH algorithm has been quite successful in tackling inference
problems, there are practical implementation issues to overcome. Improving
acceptance rate is crucial for convergence, and sometimes improvements in the
proposal distribution are not sufficient to efficiently sample the parameter
space. To tackle these issues, various MCMC enhancements have been proposed. A
prime example is the Hamiltonian Monte Carlo (HMC) algorithms that utilize
local gradients in order to generate proposal points (Neal, 2011; Betancourt &
Girolami, 2013). One variant of HMC is the No-U-Turn sampler which automates
part of the required tuning of the HMC Hoffman & Gelman (2011) sampler.
Another alternative to the “standard” MH is the Gibbs sampling algorithm,
which is particularly useful if the conditional distributions of the
parameters of the model are known (Ritter & Tanner, 1992; Muller, 1991; Gilks
& Wild, 1992). All of the above developments have been shown to be useful in
various disciplines (Gelman et al., 2004; Hogg & Foreman-Mackey, 2018; Brooks,
1998; Robert & Casella, 2011; Joseph, 2019; Sharma, 2017; Kendall et al.,
2005; Sivia & Skilling, 2006; Sorensen & Gianola, 2007; Johannes & Polson,
2010; Baio, 2012). Finally, there have recently been numerous proposals that
aim to enhance sampling with machine learning techniques. At their core, many
of these methods optimize the exploration of the likelihood surface, either by
learning it directly (see for example Hermans et al. (2019)) or by sampling it
in a simpler latent space (for example Wu et al. (2020)).

In this work, we introduce Eryn, which is built around the emcee package
(Foreman-Mackey et al., 2013), enhanced with a variety of sampling mechanisms
that allow us to perform inference on dynamical parameter spaces with minimal
tuning. We expand on the most important features in the sections below.

###  2.1 Affine-invariant samplers

An affine transformation is one of the form
$\vec{\theta}\rightarrow\vec{\zeta}=A\vec{\theta}+b$, where $A$ and $b$ are a
constant matrix and vector respectively. Under an affine transformation a
probability density $p(\vec{\theta}|y)$ transforms to

$p_{A,\,b}(\vec{\zeta}|y)=p(A^{-1}(\vec{\zeta}-b)|y)/{\rm det}(A).$ (5)

Such transformations can help to transform difficult-to-sample distributions
into easier-to-sample ones. A simple example is a multivariate Normal
distribution. If the dynamical range of the eigenvalues of the covariance
matrix is very large, then sampling can be difficult, but any multivariate
Normal distribution can be transformed into a spherical distribution via an
affine transformation.

Affine-invariant MCMC is a class of samplers that are designed to have equal
sampling efficiency for all distributions that are related by an affine
transformation (Foreman-Mackey et al., 2013; Goodman & Weare, 2010). The
sequence of samples in a Markov chain, $\\{X(t)\\}$, can be written
deterministically as a function of a sequence of random variables, $\xi(t)$,
which represent the random draws used to propose new points and evaluate the
accept/reject decision. Specifically we can always write

$X(t+1)=R(X(t),\xi(t),p)$ (6)

where $p$ denotes the target density. An affine-invariant sampler has the
property

$R(AX(t)+b,\xi(t),p_{A,b})=AR(X(t),\xi(t),p)+b,$ (7)

i.e., the sequence of points visited when sampling an affine-transformed
density are the affine transformations of the states visited when sampling the
original density. If an affine transformation exists that maps the given
target density to one which is more straightforward to sample from, an affine-
invariant sampler should sample it as efficiently as it could the simpler
distribution, so the convergence of affine-invariant samplers is less affected
by correlations between the parameters (Foreman-Mackey et al., 2013).

In practice, this goal is achieved by following an ensemble of points, called
walkers, and basing proposed moves on the distribution of other points in the
ensemble. In Foreman-Mackey et al. (2013), the primary update move is the so-
called stretch-move proposal. Each walker at state $X_{i}(t)$ is updated by
randomly selecting another walker $j$ and proposing a new value
$Y=X_{j}(t)+Z[X_{i}(t)-X_{j}(t)]$, where $Z$ is a random variable drawn from
the distribution (Goodman & Weare, 2010)

$g(z)\propto\left\\{\begin{array}[]{ll}\frac{1}{\sqrt{z}}&z\in\left[\frac{1}{a},a\right]\\\
0&\mbox{otherwise}\end{array}\right..$ (8)

The parameter $a$ can be tuned to improve convergence, but $a=2$ works well in
the majority of applications (Foreman-Mackey et al., 2013). The proposed point
is accepted with probability

$\alpha(X_{i}(t),Y)=1\wedge\left\\{z^{d-1}\frac{p(Y)}{p(X_{i}(t))}\right\\}$
(9)

where $p$ is the target density and $d$ is the dimension of the parameter
space. This acceptance probability is specific to the stretch proposal
distribution given by Eq. (8). For other stretch proposals, the term $z^{d-1}$
must be replaced by $z^{d-2}g(1/z)/g(z)$. Following this scheme, detailed
balance is maintained, and it can be proven that affine-invariant samplers
converge faster to their target distribution (Foreman-Mackey et al., 2013).
Below in section 3 we discuss the extension of the stretch-move proposal to
Reversible-Jump MCMC methods. The benefits of running MCMC chains in parallel,
combined with a proposal distribution that requires almost no tuning, have
contributed to an increasing popularity of affine-invariant samplers. In
particular, the emcee package (Foreman-Mackey et al., 2013), has been used
widely in Astrophysics and Cosmology (e.g Virtanen et al., 2020; Wong et al.,
2020; McMillan, 2017; De et al., 2018).

###  2.2 Delayed Rejection

The Delayed Rejection (DR) scheme of sampling was devised in order to improve
two aspects of MCMC algorithms. First, it allows for improvements in the
acceptance rate of the proposals, yielding “healthier” parameter chains, with
better mixing. Secondly, it is more robust against becoming trapped in local
maxima of the posterior surface (Trias et al., 2009; Mira, 2001; Green & Mira,
2001; Haario et al., 2006). The strategy, as the name suggests, is, at each
iteration, instead of immediately rejecting the newly proposed point based on
algorithm 1, we keep proposing new points while maintaining detailed balance
by computing both the forward and backward transition probabilities. Suppose
we are at a point $\vec{\theta}_{0}$ and use a proposal
$q(\vec{\theta}_{1}|\vec{\theta}_{0})$ to propose a new point
$\vec{\theta}_{1}$. The usual acceptance probability, following the notation
of eq. (1), is

$\alpha_{1}(\vec{\theta}_{0},\vec{\theta}_{1})=1\wedge\left\\{\frac{p(\vec{\theta}_{1}|y)q(\vec{\theta}_{0}|\vec{\theta}_{1})}{p(\vec{\theta}_{0}|y)q(\vec{\theta}_{1}|\vec{\theta}_{0})}\right\\},$
(10)

as per eq. (4). If $\vec{\theta}_{1}$ is rejected, then instead of going back
to step 1 of algorithm 1, we propose a new point, $\vec{\theta}_{2}$, drawn
from a proposal distribution
$q(\vec{\theta}_{2}|\vec{\theta}_{1},\vec{\theta}_{0})$. This proposal
distribution may depend only on $\vec{\theta}_{1}$, but we write it more
generally here to allow for the case that the proposal is adapted based on the
sequence of steps that have been rejected. The acceptance probability for
$\vec{\theta}_{2}$,
$\alpha_{2}(\vec{\theta}_{0},\vec{\theta}_{1},\vec{\theta}_{2})$, is computed
using

$\displaystyle\alpha_{2}(\vec{\theta}_{0},\vec{\theta}_{1},\vec{\theta}_{2})=$
(11) $\displaystyle
1\wedge\left\\{\frac{p(\vec{\theta}_{2}|y)q(\vec{\theta}_{1}|\vec{\theta}_{2})q(\vec{\theta}_{0}|\vec{\theta}_{1},\vec{\theta}_{2})\left[1-\alpha_{1}\left(\vec{\theta}_{2},\vec{\theta}_{1}\right)\right]}{p(\vec{\theta}_{0}|y)q(\vec{\theta}_{1}|\vec{\theta}_{0})q(\vec{\theta}_{2}|\vec{\theta}_{1},\vec{\theta}_{0})\left[1-\alpha_{1}\left(\vec{\theta}_{0},\vec{\theta}_{1}\right)\right]}\right\\}.$

If $\vec{\theta}_{2}$ is rejected, further steps can be included and each step
adds additional proposal and rejection-probability terms to the numerator and
denominator of the acceptance probability. For example, the three step
acceptance probability,
$\alpha_{3}(\vec{\theta}_{0},\vec{\theta}_{1},\vec{\theta}_{2},\vec{\theta}_{3})$
is the minimum of one and

$\displaystyle\frac{p(\vec{\theta}_{3}|y)q(\vec{\theta}_{2}|\vec{\theta}_{3})q(\vec{\theta}_{1}|\vec{\theta}_{2},\vec{\theta}_{3})q(\vec{\theta}_{0}|\vec{\theta}_{1},\vec{\theta}_{2},\vec{\theta}_{3})}{p(\vec{\theta}_{0}|y)q(\vec{\theta}_{1}|\vec{\theta}_{0})q(\vec{\theta}_{2}|\vec{\theta}_{1},\vec{\theta}_{0})q(\vec{\theta}_{3}|\vec{\theta}_{2},\vec{\theta}_{1},\vec{\theta}_{0})}$
$\displaystyle\hskip
14.22636pt\times\frac{\left[1-\alpha_{1}\left(\vec{\theta}_{3},\vec{\theta}_{2}\right)\right]\left[1-\alpha_{2}\left(\vec{\theta}_{3},\vec{\theta}_{2},\vec{\theta}_{1}\right)\right]}{\left[1-\alpha_{1}\left(\vec{\theta}_{0},\vec{\theta}_{1}\right)\right]\left[1-\alpha_{2}\left(\vec{\theta}_{0},\vec{\theta}_{1},\vec{\theta}_{2}\right)\right]}$
(12)

The proposal $q$ can be different at each step, as long as the relevant
proposal density is used in eq. (11). For example, in Trias et al. (2009) the
proposal is built upon a Gaussian mixture model that tries further points in
the parameter space with the aim of efficiently exploring multiple modes of
the posterior distribution. As the number of steps in the DR scheme becomes
arbitrarily large the acceptance probability slowly approaches zero. This
algorithm is also limited in practice by high computational requirements,
since at every delayed rejection step we need to evaluate a new likelihood and
compute the backwards probability (the
$\alpha_{1}(\vec{\theta}_{2},\vec{\theta}_{1})$ from eq. (11)). Nevertheless,
the DR scheme offers many advantages, and despite the computational cost, it
is very useful when the posterior surface exhibits high dimensionality, and
when acceleration techniques are available. These, for example, might include
the use of Graphical Processing Units (GPUs), and/or heterodyned likelihoods
(Cornish, 2021). In our implementation here, we follow closely the one in
Trias et al. (2009), for improving the acceptance rate of the between-model
step of the Reversible Jump algorithm (see section 2.6). As already mentioned,
the Reversible Jump MCMC allows for sampling dynamical parameter spaces. In
the special case of nested models, such as the case of searching multiple
signals in the LISA data, proposing the ‘birth’ of a signal out of a very wide
prior can be very inefficient. A delayed rejection scheme alleviates this
problem, by effectively performing a small search around the first set of
rejections, increasing the chances of finding a good signal candidate, and
thus improving the mixing of the chains.

###  2.3 Multiple Try Metropolis

The Multiple Try Metropolis (MTM) (Martino, 2018; Liu et al., 2000; Martino et
al., 2011; Bédard et al., 2012) is a sub-class of the implementation of the MH
algorithm, which is based on the idea of generating a number of proposals for
each individual current state, and then selecting one of them based on their
importance weight. In proposing a move from $\vec{\theta}_{t-1}$, a set of $N$
possible new points, $\\{y_{k}\\}$, are drawn from a proposal distribution
$q(y|\vec{\theta}_{t-1})$ and are assigned weights
$w_{i}=w(y_{i}|\vec{\theta}_{t-1})$ using a weight function
$w(\cdot|\vec{\theta}_{t-1})$. One of these proposed points, $y_{J}$, is
selected with probability given by the normalised weight

$\bar{w}_{i}=\frac{w_{i}}{\sum_{k=1}^{N}w_{k}}.$ (13)

To compute the acceptance probability we need to draw $N-1$ points,
$\\{x_{i},i=1,\ldots,N-1\\}$, for the reverse move from the proposal
$q(x|y_{J})$, and assign weights $w(x|y_{J})$. We then set
$\vec{\theta}_{t}=y_{J}$ with probability

$\alpha(\vec{\theta}_{t-1},y_{J})=1\wedge\left\\{\frac{w(y_{J}|\vec{\theta}_{t-1})+\sum_{k=1,k\neq
J}^{N}w(y_{k}|\vec{\theta}_{t-1})}{w(\vec{\theta}_{t-1}|y_{J})+\sum_{k=1}^{N-1}w(x_{k}|y_{J})}\right\\}.$
(14)

and set $\vec{\theta}_{t}=\vec{\theta}_{t-1}$ otherwise (Martino, 2018). This
procedure will satisfy detailed balance if the weight function is chosen such
that

$p(\vec{\theta}_{0}|y)q(\vec{\theta}_{1}|\vec{\theta}_{0})w(\vec{\theta}_{1}|\vec{\theta}_{0})=p(\vec{\theta}_{1}|y)q(\vec{\theta}_{0}|\vec{\theta}_{1})w(\vec{\theta}_{0}|\vec{\theta}_{1}).$
(15)

This will be satisfied by a weight function of the form

$w(\vec{\theta}_{t}|\vec{\theta}_{t-1})=p(\vec{\theta}_{t}|y)q(\vec{\theta}_{t-1}|\vec{\theta}_{t})\xi(\vec{\theta}_{t-1},\vec{\theta}_{t}),$
(16)

where $\xi(\vec{\theta}_{t-1},\vec{\theta}_{t})$ is any symmetric function,
i.e.,
$\xi(\vec{\theta}_{t-1},\vec{\theta}_{t})=\xi(\vec{\theta}_{t},\vec{\theta}_{t-1})$,
$\forall\vec{\theta}_{t},\vec{\theta}_{t-1}\in{\cal
D}\subseteq\mathbb{R}^{d}$, with $d$ being the dimensionality of the problem
at hand. The detailed balance condition can also be satisfied by a weight
function of the form

$w(\vec{\theta}_{t}|\vec{\theta}_{t-1})=\frac{p(\vec{\theta}_{t}|y)}{q(\vec{\theta}_{t}|\vec{\theta}_{t-1})}.$
(17)

Making this choice and additionally using a proposal function that is
independent of the current point,
$q(\vec{\theta}_{t}|\vec{\theta}_{t-1})=q(\vec{\theta}_{t})$ only, we obtain
the Independent MTM algorithm (Martino, 2018). When using the independent MTM
algorithm detailed balance is maintained when the same set of points is used
for the reverse proposal as for the forward proposal, which saves the
evaluation of $N-1$ posterior densities.

The base MTM is currently implemented in Eryn with options for the Independent
MTM algorithm and symmetric proposals. For a symmetric proposal distribution,
$q(\vec{\theta}_{t-1}|\vec{\theta}_{t})=q(\vec{\theta}_{t}|\vec{\theta}_{t-1})$,
eq. (15) can be satisfied using the weight function
$w(\vec{\theta}_{1}|\vec{\theta}_{0})=p(\vec{\theta}_{1}|y)$. In this case, we
still need to draw separate samples for the reverse step (unlike in the
Independent MTM case).

Generating a large number of candidate points yields certain advantages. As
expected, the first advantage is the fact that there is usually very good
coverage of the parameter space. The second is that the implementation of the
MTM usually results in chain states with very low correlation between them.
Nevertheless, as for Delayed Rejection, this algorithm requires increased
computational resources, since multiple likelihoods/posterior densities have
to be evaluated at each iteration of the chain. This cost can be offset in
cases where the computations can be parallelized, for example using either CPU
or GPU acceleration.

###  2.4 Adaptive Parallel Tempering

The concept of Parallel Tempering was introduced in order to efficiently
sample surfaces with high multi-modality (Swendsen & Wang, 1986; Hukushima &
Nemoto, 1996; Vousden et al., 2016). The idea is based on a transformation of
the posterior density to a density with a different temperature, $T$, defined
by

$p_{T}(\vec{\theta}|y)\propto p(y|\vec{\theta})^{1/T}p(\vec{\theta}).$ (18)

For $T=1$ this is the target posterior density. In the limit
$T\rightarrow\infty$ it is the prior density. Intermediate temperatures
“smooth out” the posterior by reducing the contrast between areas of high and
low likelihood.

In parallel tempering, a set of Markov chains are constructed in parallel,
each one sampling the transformed posterior for a different temperature $T$.
These chains periodically exchange information. The idea is that the hottest
chains explore the parameter space more widely, and information about areas of
high likelihood that they encounter propagate to the colder chains.
Information is exchanged by proposing swaps of the states between the
different chains. If two chains are sampling from target densities
$p_{1}(\vec{\theta})$ and $p_{2}(\vec{\theta})$ respectively, then the
transition probability for chain $1$ in the swap is
$p_{1}(\vec{\theta}_{0}\rightarrow\vec{\theta}_{1})=p_{2}(\vec{\theta}_{1})\alpha(\vec{\theta}_{0},\vec{\theta}_{1})$.
Detailed balance is thus maintained by accepting the swap with probability

$\alpha(\vec{\theta}_{0},\vec{\theta}_{1})=1\wedge\left\\{\frac{p_{1}(\vec{\theta}_{1})p_{2}(\vec{\theta}_{0})}{p_{1}(\vec{\theta}_{0})p_{2}(\vec{\theta}_{1})}\right\\},$
(19)

which for the specific case of swapping between two tempered chains $i$ and
$j$ when doing parallel tempering is

$\alpha_{i,j}=1\wedge\left\\{\left(\frac{p(y|\vec{\theta}_{i})}{p(y|\vec{\theta}_{j})}\right)^{\beta_{j}-\beta_{i}}\right\\},$
(20)

with $\beta_{i}=1/T_{i}$ being the inverse temperature, and $\vec{\theta}_{i}$
the given parameter state for the $i$-th chain.

The temperature ladder $T_{i}$ should be chosen in order to maximize the
information flow between chains of different temperatures, so as to encourage
the efficient exploration of the complete parameter space. Typically, this
ladder can be static or dynamically adjusted during the sampling procedure. In
Eryn we have adopted the procedure of Vousden et al. (2016), which adapts the
temperature ladder based on the swap acceptance rate calculated directly from
the chains. Ideally, one should aim for equal acceptance ratio between every
pair of neighboring tempered chains, thus tuning their log-temperature-
difference $S_{i}\equiv\log(T_{i}-T_{i-1})$, according to the swap acceptance
rate from eq. (20):

$\frac{{\mathrm{d}}S_{i}}{{\mathrm{d}}t}=\kappa(t)\left[\alpha_{i,i-1}(t)-\alpha_{i+1,i}(t)\right],$
(21)

where $\kappa(t)$ tunes the timescale of the evolution of the temperatures.
The function $\kappa(t)$ can be chosen depending on the desired behavior of
the procedure. In Vousden et al. (2016) a hyperbolic dependence on the $t$
state is chosen, in order to suppress large dynamic adjustments on long
timescales. This setup is the default option in Eryn, but it can be
customized. This process is more straightforward for ensemble samplers, where
multiple walkers are used, simply because one can get an estimate of the
acceptance rate directly from the particular state of the walkers at any given
time step $t$. Otherwise, the acceptance rate is computed after iterating for
a predefined number of steps, which can be chosen by the user for the given
problem at hand. It can be proven (Vousden et al., 2016), that the temperature
ladder will converge to a particular stable configuration. One should only use
this scheme for the initial burn-in stage of sampling, and then continue with
a stationary ladder for the rest of the analysis.

###  2.5 Marginal posterior calculation for model selection

One of the most frequently encountered problems in physics, and in science in
general, is that of model or variable selection, i.e., identifying the model
best supported by the observed data. Working in a Bayesian framework,
comparison between different hypotheses may be done by computing their
evidences or marginal posteriors (Gelman et al., 2004) and evaluating the
Bayes Factor:

$B_{12}=\frac{p(\vec{y}|\mathcal{M}_{1})p(\mathcal{M}_{1})}{p(\vec{y}|\mathcal{M}_{2})p(\mathcal{M}_{2})},$
(22)

where the term $p(\mathcal{M}_{i})$, is the prior probability assigned to the
model ${\cal M}_{i}$.

The marginal posterior density, or evidence, is given by the integral of eq.
(2) and is in general quite challenging to compute. For some high signal-to-
noise ratio (SNR) cases it can be approximated if the covariance matrix
$\Sigma$ of the parameters for all candidate models $\mathcal{M}$ are known.
This approach is called the Laplace approximation (Kass & Raftery, 1995;
Gelman et al., 2004). However, this is only an approximation, and it sometimes
fails for models with weak support from the data (Cornish & Littenberg, 2007)
(in particular when the posterior cannot be approximated by a multivariate
Gaussian at $\vec{\theta}_{\mathrm{MAP}}$).

When using parallel tempering 2.4, it is possible to compute the evidence by a
procedure known as thermodynamic integration (Lartillot & Philippe, 2006). We
define a continuous distribution of evidences based on the target distribution
for a chain with inverse temperature $\beta=1/T$ via

$\mathcal{Z}_{i,\beta}=\int
p(y|\vec{\theta},\mathcal{M}_{i})^{\beta}p(\vec{\theta})\,{\rm
d}\vec{\theta}.$ (23)

For $\beta=0$ the chain is sampling the prior and therefore
$\log\mathcal{Z}_{i,0}=0$. For $\beta=1$ we are sampling the target
distribution and $\log\mathcal{Z}_{i,1}=\log\mathcal{Z}_{i}$. Additionally we
have

$\displaystyle\frac{{\rm d}\log\mathcal{Z}_{\beta}}{{\rm d}\beta}$
$\displaystyle=\int\log[p(y|\vec{\theta},\mathcal{M}_{i})]\,p(y|\vec{\theta},\mathcal{M}_{i})^{\beta}p(\vec{\theta}){\rm
d}\vec{\theta}$ $\displaystyle\equiv\mathbb{E}_{\beta}[\log
p(y|\vec{\theta},\mathcal{M}_{i})].$ (24)

From this we deduce

$\log\mathcal{Z}_{i}=\int_{0}^{1}\mathbb{E}_{\beta}\left[\log
p\left(y|\vec{\theta},\mathcal{M}_{i}\right)\right]d\beta.$ (25)

The expectation value is over the distribution being sampled by the chain at
temperature $\beta$ and so can be computed by averaging over the posterior
samples (Lartillot & Philippe, 2006; Goggans & Chi, 2004; Vousden et al.,
2016). The integral can then be evaluated using standard methods, for example
the trapezium rule, using the full ladder of temperatures. This approach
generates reliable evidences, with accuracy limited only by the number of
temperatures being sampled, and the efficiency of the sampling of the
parameter space $\Theta$ by the chains. Since its first introduction, there
have been many applications of this approach, and in particular, there is
extensive usage in GW astronomy (Vousden et al., 2016; Katz et al., 2022;
Maturana-Russel et al., 2019; Littenberg & Cornish, 2009, 2010).

The thermodynamic integral in Eq. (25) can be thought of as computing the
evidence as a sum of differences between the evidences at different
temperatures. An alternative approach, called the stepping-stone algorithm
(Xie et al., 2010), writes the evidence as a product of the ratios of
evidences at different temperatures

$\mathcal{Z}_{i}\equiv\frac{\mathcal{Z}_{i,1}}{\mathcal{Z}_{i,0}}=\prod_{k=1}^{N_{T}-1}\frac{\mathcal{Z}_{i,\beta_{k+1}}}{\mathcal{Z}_{i,\beta_{k}}}$
(26)

where $\beta_{k}$ denotes the inverse temperature of chain $k$, $N_{T}$ is the
number of different temperatures being sampled, and we assume $\beta_{1}=0$
and $\beta_{N_{T}}=1$. Each evidence ratio can be written as a posterior
integral

$\displaystyle\frac{\mathcal{Z}_{i,\beta_{k+1}}}{\mathcal{Z}_{i,\beta_{k}}}$
$\displaystyle=\int
p(y|\vec{\theta},\mathcal{M}_{i})^{\beta_{k+1}-\beta_{k}}\frac{p(y|\vec{\theta},\mathcal{M}_{i})^{\beta_{k}}p(\vec{\theta})}{\mathcal{Z}_{\beta_{k}}}\,{\rm
d}\vec{\theta}$
$\displaystyle=\mathbb{E}_{\beta}\left[p\left(y|\vec{\theta},\mathcal{M}_{i}\right)^{\beta_{k+1}-\beta_{k}}\right]$
$\displaystyle\approx\frac{1}{n}\sum_{i=1}^{n}p\left(y|\vec{\theta}^{i}_{k},\mathcal{M}_{i}\right)^{\beta_{k+1}-\beta_{k}},$
(27)

where $n$ is the number of posterior samples in each chain, and
$\vec{\theta}^{i}_{k}$ denotes the $i$’th sample at temperature $\beta_{k}$.
This leads to the final expression

$\log\mathcal{Z}=\sum_{k=1}^{N_{T}-1}\log\sum_{i=1}^{n}p(y|\vec{\theta}^{i}_{k},\mathcal{M})^{\beta_{k+1}-\beta_{k}}-(N_{T}-1)\log
n.$ (28)

In challenging situations, e.g., where the number of tempered chains is
relatively small, the stepping-stone algorithm has been shown to produce more
accurate estimates of the marginal likelihood than methods that use
thermodynamic integration (Xie et al., 2010). This was also demonstrated with
practical examples drawn from Gravitational Wave Astronomy in Maturana-Russel
et al. (2019).

###  2.6 Reversible Jump

Another approach to the model selection problem is to follow a Reversible Jump
(RJ) MCMC strategy, which can dynamically estimate the most probable
hypotheses given the data (Green, 1995). The RJ-MCMC is a generalization of
the MH algorithm that allows trans-dimensional proposals. Thus, the model
order is considered a free parameter which is fitted together with the
parameters of the individual models. The most widespread variation of the
algorithm uses a two-stage procedure. The first stage or in-model step, uses
the standard MH algorithm to update all the parameters $\vec{\theta}_{k}$ for
the given model $k$. The second stage or between-model step proposes to update
the model state $k$ to a new model state $l$. Parameters $\vec{\theta}_{l}$
for the new model are also proposed. The newly proposed state $l$, is accepted
with a probability defined by (Godsill, 2001)

$\alpha^{\prime}=1\wedge\left\\{\frac{p(\vec{\theta}_{l}|l,y)g^{\prime}(u_{l})}{p(\vec{\theta}_{k}|k,y)g(u_{k})}|{\bf
J}|\,\right\\},$ (29)

where

$p(\vec{\theta}_{k}|k,y)=p(y|\vec{\theta}_{k},k)p(\vec{\theta}_{k},k)p(k)$
(30)

with $p(y|\vec{\theta}_{k},k)$ the likelihood for model $k$,
$p(\vec{\theta}_{k}|k)$ the prior on the parameters $\vec{\theta}_{k}$ in
model $k$, and $p(k)$ the prior for the model state $k$. The term
$g^{\prime}(u_{l})|{\bf J}|/g(u_{k})$ arises because of the need for dimension
matching between the different model states. In general, we can define a move
between model states in terms of a deterministic invertible mapping,
$\vec{\theta}_{k}=q(\vec{\theta}_{l},u_{l})$ with inverse
$\vec{\theta}_{l}=q^{\prime}(\vec{\theta}_{k},u_{k})$, that is a function of
the parameters and two sets of random variables, $u_{k}$ and $u_{l}$, drawn
from distributions $g(u_{k})$ and $g^{\prime}(u_{l})$ (Godsill, 2001; Green,
1995). The term $|{\bf J}|$ is the Jacobian defined by this invertible mapping

$|{\bf
J}|=\left|\frac{\partial(\vec{\theta}_{l},u_{l})}{\partial(\vec{\theta}_{k},u_{k})}\right|$
(31)

and the term $g^{\prime}(u_{l})/g(u_{k})$ plays the role of the proposal ratio
in the standard MH acceptance probability. Dimension mapping means that ${\rm
dim}(\vec{\theta}_{l})+{\rm dim}(u_{l})={\rm dim}(\vec{\theta}_{k})+{\rm
dim}(u_{k})$.

Using RJ-MCMC introduces additional computational cost at each MCMC iteration,
as well as technical challenges in implementation. Luckily, implementation can
be easier when sampling nested models. This refers to problems where more
complicated models contain their simpler counterparts. Examples of such cases
are fitting polynomial models, which differ only in the highest order to be
determined, or detection problems where multiple similar signals are
potentially present the data. In such cases, the between-model step can always
be formulated such that the Jacobian of eq. (31) becomes unity, and eq. (29)
simplifies to the ratio of posteriors accounting for any differences in prior
and proposal volumes (Dellaportas et al., 2002; Lopes & West, 2004; Littenberg
et al., 2020).

After running RJ-MCMC, the Bayes Factor can be approximated by the ratio of
the number of iterations spent within each model:

$B_{\rm
12}=\frac{\rm\\#\>of\>iterations\>in\>model\>{\mathcal{M}}_{1}}{\rm\\#\>of\>iterations\>in\>model\>{\mathcal{M}}_{2}}.$
(32)

This algorithm has proven to be robust for evaluating high-dimensional
competing models, and has been quite successful in tackling data analysis
problems in GW astronomy (Littenberg et al., 2020; Cornish & Littenberg, 2007;
Karnesis et al., 2014) as well as areas spanning physics and signal processing
(e.g., Marrs, 1997; Khan et al., 2005; Yu et al., 2021). However, designing an
efficient RJ-MCMC algorithm can be quite challenging. The first challenge is
to choose suitable proposal distributions, which can greatly affect the
convergence of the algorithm. In situations where the models are nested, it is
both tempting and convenient to take the proposal to be the same as the prior
distribution of the parameters. As an illustrative example, we refer to
section 3.3.1, which describes a toy problem of searching for Gaussian pulses
in noisy data. There, the parameters of the individual pulses are the
amplitude and location of the pulse described by their $(x,\,y)$ coordinates.
In order to search for those signals, the prior on their location must be wide
enough to include the complete data set (see figure 1(a)). A birth proposal
based on the prior would inevitably be quite inefficient, simply because the
chance of proposing a good source candidate is small, especially if the
proposal distribution is flat across $(x,\,y)$. We treat the above problem as
a motivation to adopt efficient proposals with minimal tuning in Eryn, which
we further discuss in section 3. The second major challenge, which of course
depends on the given problem at hand, arises from the samplers’ capability to
explore a multi-modal dynamical parameter space. We discuss our strategy to
overcome that challenge in section 3.

##  3 Eryn: Gathering all the pieces together

All the different algorithms described in previous sections can be extremely
useful in tackling different kinds of problems that require sampling. In
Gravitational Wave Astronomy, we encounter such problems far too often, where
dynamical parameter spaces require vast computational resources in order to be
explored efficiently. Motivated by those problems, we have implemented a new
toolbox that combines all these techniques to enhance the capabilities of an
MCMC sampler. We have named this package Eryn 222Eryn Vorn (Sindarin for
Blackwood) was a wooded cape in Eriador, and a region of dark pine trees.
Located in western Minhiriath, Eryn Vorn (likely named so by the Númenóreans )
was originally part of the vast ancient treescape that covered most of north-
western Middle-earth, borrowing the name from the Tolkien mythos
(TolkienGateway). The analogy has its basis in the idea of a forest: within a
forest you have trees which correspond to different walkers, a.k.a. Ents, in
an ensemble MCMC sampler. On each tree there are branches that represent the
various types of models used to fit the data. For example, in the case of GW
global fitting for LISA, we can imagine using the Galactic binaries as one
branch and massive black hole binaries as another branch. Each branch has
leaves, which represent the individual instances of each model. In the LISA
example, leaves would represent the individual Galactic binaries or massive
black hole binaries. And finally, to zoom out, when adding tempering
capabilities, we can think of groups of walkers in each temperature taking the
form of many forests (of walkers) located within different temperature
climates.

We adopt the architecture of ensemble samplers, and in particular the one of
emcee (Foreman-Mackey et al., 2013). Having multiple walkers running in
parallel is ideal for efficiently sampling the parameter spaces using
techniques such as parallel tempering, as described in 2.4 (also see section
2.1). In this setting, we evolve $n_{w}$ walkers per temperature $T_{i}$,
where each walker follows a Reversible Jump MCMC (see section 2.6), mapping a
parameter space of altering dimensionality. In practice, walkers in higher
temperatures sample the dynamic parameter space with fewer model components as
the penalty from higher prior volume is not compensated by the smoother
annealed likelihood. In other words higher temperatures have a sharper Occam’s
razor: the data can be explained with models that are simpler, or lower-
dimensional. The highest temperature chain samples the prior on the model
space (provided that $T_{\mathrm{max}}=\infty$). More details will be given in
section 3.3.

As already mentioned in section 2.6, Reversible Jump algorithms are extremely
challenging to tune, even for simpler classes of problems. One of the major
challenges is the low acceptance rate for the between-model proposal, i.e.,
when we propose a new state where the parameter dimensionality differs. In
cases of signal search and detection (which is a nested model situation), it
is convenient to set the proposal corresponding to a “birth” move to be the
same as the prior distribution. In order to accommodate all possibilities for
the signals present, the prior densities are usually quite wide, and thus
accepting a new higher dimensional state becomes quite improbable. For that
reason, within Eryn, we have implemented a Delayed Rejection scheme with the
aim of improving this acceptance ratio. When proposing $\vec{\theta}_{l}$ for
a higher-dimensional model $l$, we do not reject immediately, but rather make
new delayed rejection proposals around the first rejected point
$\vec{\theta}_{l}$, using the given in-model step update proposal. This, in
principle, allows the sampler to explore around $\vec{\theta}_{l}$ before
rejecting the new state (Trias et al., 2009), which in turn improves the _b_
etween-model step acceptance rate and produces healthier MCMC chains.

The Delayed Rejection scheme, as described in section 2.2, requires a serial
computation of the delayed rejection acceptance ratio for walkers where the
newly proposed state has been rejected. This scheme of calculating costly
likelihoods sequentially in a loop during the between-model step, can lead to
a computational bottleneck of the MCMC process. This is especially true for
the LISA Global Fit problem, where multiple binary waveform signals are
present in the data-stream. Then, the computational time for each RJ-MCMC
iteration is significantly increased, since the progress will be halted until
all walkers have gone through their respective Delayed Rejection process,
which requires evaluation of new waveforms at each step. For the reasons
summarized here, we have not used the Delayed-Rejection scheme for our
analysis in section 4, and have resorted to the Multiple Try scheme. However,
the Delayed-Rejection scheme, as explained in section 2.2, has been
implemented in the Eryn package.

The Multiply Try scheme was essentially implemented in order to facilitate use
of a parallelized likelihood framework. Parallelization is naturally
compatible with Multiple Try MCMC as multiple proposals are made for each
individual walker, which allows for the parallelized evaluation of proposal
distributions, likelihood functions, and acceptance probabilities. Under these
parallelized settings, one proposal can act as many when compared to the usual
serial evaluation of proposals, allowing for better chain mixing in situations
where proposals are infrequently accepted. That being said, it is still
important to choose a good proposal distribution, for both the in-model and
between-model RJ-MCMC steps, which we discuss further below.

###  3.1 Choosing efficient Proposal distributions

In the previous sections, we briefly discussed some of the challenges in
choosing efficient proposal distributions for both the in-model and between-
model steps of the RJ algorithm. For the in-model case, the challenge arises
from the fact that it is sometimes impractical, or even unfeasible, to define
a well-tuned proposal for each of the possible models that could represent the
data. Using again the example of LISA data, one would need to tediously design
an effective proposal distribution for the thousands of overlapping binary
signals in the data. On the other hand, for the case of the between-model
step, choosing proposals from the prior distribution, especially if it is
highly uninformative, can be very inefficient for Reversible Jump sampling.
For Eryn, in order to tackle those issues, we have implemented the Group
Proposals explained below for addressing the within-model proposals in RJ, as
well as a scheme to design an efficient proposal for birth moves during the
between-model step, which is based on normalizing flows.

####  3.1.1 Group Proposals

In section 2.1, the stretch-move proposal was introduced and discussed. One of
the obvious advantages of such a scheme of proposing new MCMC samples is that
it requires minimal tuning (Foreman-Mackey et al., 2013). However, it does not
extend well in its simplest form to the generalized Reversible Jump MCMC. The
stretch proposal is based on the idea that the ensemble of points ($X_{j}$) is
sitting on the same posterior mode as the current point ($X_{i}$). In a nested
model situation where both the model count and the individual model parameters
change, each point may lie on a different posterior mode representing a
different set of leaves (sources) in the data. This can be alleviated by
applying the stretch move to individual leaves within each branch of each
walker, but there is still an issue of identifying leaves in different walkers
that lie on the same posterior mode. The stretch proposal will technically
still work when mixing leaves in different posterior modes, but the acceptance
fraction will be negatively affected. However, within the stretch proposal
formalism, the choice of $X_{j}$ is customizable. The key to maintaining
detailed balance is that $X_{j}$ cannot depend on $X_{i}$, and $X_{j}$ cannot
be updated in the same iterative step as $X_{i}$ (Goodman & Weare, 2010).

We leverage this property to design a new type of stretch move that can handle
Reversible Jump setups while maintaining a small number of tuning parameters.
We call this proposal the “group” proposal. The mathematics that govern the
group proposal are equivalent to that of the original stretch proposal. The
key difference is that the group proposal chooses $X_{j}$ (see Section 2.1)
from a stationary group that is fixed for many proposed updates. This is in
contrast to the original stretch proposal that always uses the current set of
points in the ensemble sampler to draw $X_{j}$. The stationary group is
updated after a large number of sampler iterations and we make sure that
detailed balance is maintained during the update. We update every leaf within
every branch of every walker at each iteration and repeat many iterations
between updates of the stationary group.

The appropriate stationary group varies from problem to problem. The goal is
to set a group that resembles as best as possible the posterior modes of the
current leaves and then draw from it strategically so that the drawn point is
likely (but not guaranteed) to lie on the same mode as the leaf that is
currently being updated, $X_{i}$. In the example of the LISA Galactic binaries
analysis, we set our stationary group to the full set of leaves (binaries)
across all walkers at a specific temperature of the sampler at the end of a
given iteration. Then, at proposal time, we efficiently locate the $\sim
n_{w}$ points in the stationary group that are closest to $X_{i}$ from based
on their initial frequency parameter. We then draw $X_{j}$ from this group.
The hope is that some percentage of the $n_{w}$ drawn points will lie on the
posterior mode on which $X_{i}$ sits. The exact percentage will vary depending
on how close the posterior modes are to each other and how many model
instances exist in the sampler that include this specific mode. For low SNR
binaries, for example, a source may exist in some walkers and not others,
making it harder to map its posterior mode with the current group of
stationary points.

The performance of group proposals is highly situation- and/or model-
dependent. With individual source posterior modes that are well separated and
easy to define in terms of separation, the performance will approach the
performance of the base stretch proposal in non-Reversible Jump MCMC because
the stationary group will well represent the specific posterior mode on which
$X_{i}$ is located. As the parameter space becomes more crowded and/or
separation (distance) metrics become harder to define, the performance of
group proposals will worsen.

####  3.1.2 Learning from the data

The second improvement concerns the between-model step of the RJ-MCMC. As
mentioned earlier, for the case of nested models, it is often convenient to
draw “birth” candidates directly from the prior distribution of parameters of
the given model. This practice can be quite ineffective in terms of acceptance
rate. As an example we can again use the LISA data-set case. The UCBs are
distributed within the Galactic disk, congregated mostly around its Center
(Amaro-Seoane et al., 2022), therefore, adopting a proposal based on an
uninformative uniform prior across the sky, would waste computational
resources exploring a part of the parameter space with low probability mass. A
proposed solution is to use an informative prior derived from the spatial
distribution of binaries in the Galaxy (Littenberg et al., 2020). In our work
here, we have chosen an alternative data-driven route, based on the actual
residual data after a burn-in period of the RJ-MCMC, which we describe below.

(a)

(b)

(c)

Figure 1: Searching for 2D Gaussian pulses in the presence of Gaussian noise.
Panel (a): the simulated data, which consists of injections of $25$ pulses in
Gaussian noise with $\sigma_{n}=0.2$. Panel (b): the distribution of the model
order, obtained by exploring the dynamical parameter space with Eryn. The true
value is marked with a dashed red line. For this toy investigation, the
correct number of simulated components is recovered. Panel (c): the 2D
posterior densities for the parameters of the $k$ Gaussian peaks (see text for
details).

After a sufficient number of RJ-MCMC iterations, we can extract a subset of
sources from nested models which are constantly present in almost all walkers
of our cold chain. In other words, we can find and subtract the brightest
sources from the data, and then allow for another burn-in period on the
resulting residuals. This allows the sampler to explore the remaining
parameter space more easily, thus providing a good initial estimate for the
weak signals possibly buried in the noisy residual data. We can then use those
samples to construct a proposal density which will help us search for good
candidates for those weak signals, without excluding the rest of the parameter
space. This can be accomplished by fitting the distribution to the residual
data described above. The most efficient way to fit to the generic
distribution, is to use an invertible transformation, such as a normalizing
flow (for example Dinh et al., 2016; Durkan et al., 2019). The method works in
the following way: we sample from the base distribution (which is usually
chosen to be Normal $\mathcal{N}(\mathbf{z};0,1)$) and transform samples to
the desired distribution $p(\bm{\theta})$ by applying the change of variable
equation

$p(\bm{\theta})=N(f^{-1}(\bm{\theta}))|J_{f^{-1}}(\bm{\theta})|.$ (33)

Here function $f(\mathbf{z})$ is a bijection which we fit by optimizing the
Kullback-Leibler divergence,
$D_{\text{KL}}[p(f(\bm{\theta}))||\mathcal{N}(\mathbf{z};0,1)]$, between a
normal distribution and the inverse transform of the distribution that we want
to estimate. After the fit has converged we can draw samples from the Normal
distribution and transform these to samples from the distribution fitted to
all residuals and use it as a proposal. We will cover this method in more
detail in a separate paper which is being prepared.

####  3.1.3 Convergence Diagnostics

Many standard approaches for assessing convergence of MCMC chains can be
applied to the RJMCMC chains generated by Eryn. Trace plots, both for
individual parameters within models, and the indices labelling different
models, will show that the chains are mixing well and that the initial burn-in
phase has finished. Repeating runs with different numbers of walkers,
different starting positions, different numbers of parallel chains or
different choices of temperatures and comparing results can be used to assess
convergence. Similarly, posteriors produced from randomly selected subsets of
chain points can be compared (using for example statistical tools such as the
Jensen-Shannon divergence (Menéndez et al., 1997)).

Additionally, Eryn computes and outputs the Potential Scale Reduction Factor
(PSRF) for each tempered chain, which provides a quantitative diagnostic of
the convergence of the result (Gelman & Rubin, 1992; Brooks & Gelman, 1998).
For a single parameter, the PSRF index is given by

$\widehat{R}=\sqrt{\frac{(d+3)\widehat{V}}{(d+3)W}},$ (34)

where $W$ the mean of the empirical variance within each chain, and
$\widehat{V}$ is the estimated variance of the all chains, assuming that the
target distribution is Gaussian. The degrees of freedom are estimated by the
method of moments as $d\approx
2\widehat{V}/\widehat{\mathrm{var}}(\widehat{V})$ (Brooks & Gelman, 1998).

The PSRF index takes values $\widehat{R}\geq 1$, with values close to unity
indicating converged MCMC chains. In the end, the quantity that really matters
is the PSRF value for the cold chain, but converged chains across for all
temperatures is a good indicators for healthy mixing and efficient temperature
swaps.

###  3.2  Implementation

In this section, we discuss the main implementation details of Eryn. We refer
the interested reader, or user, to the Eryn documentation for more exhaustive
information and examples (Documentation).

The goal of Eryn is to produce a sampler that can handle all (or most) cases
of MCMC sampling ranging from basic, non-tempered, single-model type, single-
model instance posterior estimation to the full reversible jump MCMC with
tempering, multiple model types, and adjustable model counts, as well as
everywhere in between. In the basic case, Eryn aims to be a close replica of
emcee trying to maintain as much simplicity as possible. At the complicated
end of the spectrum, Eryn attempts to provide a common interface and
underlying infrastructure for the variety of problems that may arise, allowing
the user to maintain usage of the majority of the code from project to
project, focusing on changing only the specific parts of the code that are
difficult to implement or require special treatment for each specific problem.
Since Eryn is effectively an enhanced version of the emcee package, the
overall structure of emcee is strongly maintained. Like in emcee, “State”
objects move coordinate and likelihood information around the ensemble
sampler, storing information in a similar back end object either in memory or
HDF5 files. Additionally, the interface used for adding proposals has
remained.

The various enhancements discussed in this work, including tempering,
reversible jump moves, multiply try MCMC, etc., are all implemented within the
emcee-like structure. This involved two main changes. First, the State objects
have been scaled to hold information necessary for reversible jump MCMC:
temperature information, prior information, and efficient and concise
containers for multiple types of models with an adjustable number of
individual model instances. Second, the reversible-jump proposal has been
added as a proposal base, similar to the use of the ‘MH’ or ‘RedBlue’ moves
within emcee. Beyond these main enhancements, there are also a variety of
smaller, but useful, additions to Eryn that help the user build a variety of
analysis pipelines. These include stopping or convergence functions, functions
to periodically update the sampler setup while running, objects to carry
special information through the sampler, and aids for coordinate
transformation.

###  3.3 Toy Examples

In this section, we present a series of working examples for Eryn. We begin
with simple problems, such as searching for simple signals in noisy data, with
the aim of demonstrating the performance of this toolbox in a dynamical
parameter space. The impact of the different enhancements discussed in section
3 will be assessed and discussed. Finally, in section 4, we will apply this
machinery to more realistic problems in Gravitational-Wave astronomy.

####  3.3.1  Searching for pulse signals in Gaussian noise

In this first example, we explore the capabilities of Eryn in a simplified
application, commonly encountered in physical sciences. We perform a search
for Gaussian pulses in a simulated 2D data-set, in the presence of Gaussian
noise with variance $\sigma_{n}=0.2$. We generate $25$ pulses randomly
distributed on the $x-y$ plane with all pulses contained within
$x,y\in[-10,10]$ (see figure 1(a)), and amplitude uniformly drawn from
$\mathcal{U}[0.7,\,1.5]$. The amplitude $A_{k}$ of each pulse, labelled by
$k$, is considered a free parameter to be estimated, in addition to the
Cartesian coordinates of their centres. The pulses’ width was kept fixed to
$\sigma_{p}\times\delta_{ij}$, with $\sigma_{p}=0.2$, for the sake of
simplicity. Thus, we are required to estimate $N_{p}$, the total number of
pulses in the data, and also estimate the parameters for each individual
signal $k$: $\vec{\theta}_{k}=\\{A_{k},\,x_{k},\,y_{k}\\}$. The noise variance
$\sigma_{n}$ is estimated as part of the fit.

Figure 2:  Top panel: The evolution of the temperature chains running in
parallel for the toy problem of searching for 2D Gaussian pulses in Gaussian
noise. The different colors indicate the the initial temperature chain index.
Following the parallel tempering recipe of (Vousden et al., 2016), the
temperature ladder is tuned according to eq. (21), and the chains start to
converge after $\sim 10^{4}$ iterations. Middle panel: The evolution of the
swap acceptance rate $\alpha_{i,j}$ described in eq. (20), per number of
walkers $n_{w}$. For this run we have used $n_{w}=10$ walkers. After $10^{5}$
iterations, the swap acceptance rate converges to a single (different) value
for every temperature chain. Bottom panel: The “in-model” acceptance rate per
temperature chain, given by eq. (4).

The analysis of this problem is performed using the adaptive parallel
tempering scheme of section 2.4 and the Reversible Jump MCMC proposals
(section 2.6). The in-model proposals for each model component are Gaussians,
with a diagonal covariance matrix $\Sigma=10^{-4}\delta_{ij}$. This proposal
is not tuned during sampling. The priors for the parameters are quite wide,
covering the entire range of the data, while the prior on the number of pulses
$k$ is set to $k\sim\mathcal{U}[0,\,50]$. With the above settings, we obtain
the results summarized in figures 1 and 2. In figure 1(b) we plot the most
probable number of Gaussian pulses present in the data, or in other words, the
most probable model given this particular data-set. It is clear that for the
given level of noise, it is straightforward to recover the true number of
signals. The noise variance is also estimated accurately as $\sigma_{n}=0.2\pm
2\times 10^{-3}$. In figure 1(c) we plot the posterior densities for the
parameters of all pulses recovered, while we also mark the true injected
values. Figure 1(c) shows the trans-dimensional MCMC chains “stacked” over all
samples of both model order and model-parameters. As already mentioned, in
this simplified scenario all signals have similar value for the amplitude,
thus the almost uni-modal marginal on $A_{k}$. This illustrative example is
useful as an introductory application to the more complicated case of
detection in Gravitational Wave astronomy presented below, in section 4.

In figure 2, three diagnostic quantities for this run are shown. In the top
panel, the evolution of temperatures is presented. Following the recipe of
Vousden et al. (2016), we control the distances between each temperature chain
based on their in-between swap acceptance rate, computed from eq. (21). The
tuning term $\kappa(t)$ is set to $\kappa(t)=t_{0}/\left(\nu(t+t_{0})\right)$,
with the adaptation lag $t_{0}=10^{4}$ and the adaptation time $\nu=10^{2}$.
The middle panel shows the evolution of the swap acceptance rate per number of
walkers between the chains running at different temperatures. After $\sim
10^{5}$ sampler iterations, the system converges to an equilibrium, where the
rate of swapping states reaches a single value across the temperature range.
In the bottom panel we show the acceptance rate for the in-model step of the
algorithm, for all temperatures. As expected, after temperature equilibrium at
$\sim 10^{5}$ samples, the acceptance rate converges to a (different) value
for each temperature, which is higher for higher temperatures (smoother
posterior surfaces are easier to explore).

Finally, it is interesting to investigate how the sampled dimensionality of
the problem varies at at different temperatures. In figure 3, we plot the
posterior on the number of pulses at each temperature. As expected, higher
temperature chains tend to favour lower model dimensionality and the
$T_{\infty}$ chain samples the prior on $k$. This can be attributed to the
choice of priors and “birth proposal” distributions for both the signal
parameters and $k$. The likelihood is down-weighted at higher temperatures,
making it harder to overcome the Occam penalty from including extra parameters
in the model. This means quieter sources are less likely to be added and the
preferred models have fewer sources.

Figure 3: Posterior on the number of Gaussians, $k$, at each temperature
$T_{i}$, for the toy problem of section 3.3.1. The different colors indicate
the the initial temperature chain index. Darker colors corresponds to colder
chains and vice versa.

####  3.3.2  Modelling power spectra: Searching for the optimal number of
B-spline knots.

One of the most common problems in signal processing is the characterization
of the spectra of the data. This is often done by adopting spectral models and
fitting the spectra directly in the frequency domain. This methodology is used
when the signal of interest has stochastic properties. Examples from GW
astronomy, include the measurement of stochastic signals with astrophysical,
or cosmological origin (Abbott et al., 2021c; Auclair et al., 2023). There are
many examples of possible stochastic signals for LISA (Amaro-Seoane et al.,
2017, 2012; Auclair et al., 2023; Amaro-Seoane et al., 2022). Searching for
signals with stochastic properties requires flexible spectral models, both for
the observatory instrumental noise, and the measured stochastic signal. For
these reasons, it is sometimes convenient to adopt a versatile model, such as
one that is based on B-spline interpolation schemes.

B-splines are a geometrical modeling tool, and have proven to be very useful
for modelling or generating smoother representations of data. They are piece
wise polynomial curves with a certain number of continuous derivatives, and
can be parameterised in various ways (Piegl & Tiller, 1996). For this
application, we follow Baghi et al. (2023), and choose to work with cubic-
spline interpolation, using the corresponding SciPy library (Virtanen et al.,
2020). The procedure starts by selecting a number of control points, or knots,
with a given position and amplitude, which the smooth polynomial curve crosses
and at which there is a change in the first non-continuous derivative. One of
the challenging problems using such methods, is to choose the optimal number
of spline knots for fitting the data, without overfitting. This is a model
selection problem that can be easily solved with dynamical algorithms such as
the one presented here.

For our next example, we generate time-series data directly from a theoretical
model PSD. The simulated data are represented with the solid gray line in
figure 4(a). We then use the machinery of Eryn to find the optimal
dimensionality for the problem, together with the best-fit parameters for the
knots. To ease the computational complexity, we compute the PSD of the time
series using the methodology developed by Tröbs & Heinzel (2006). In more
detail, we begin by choosing a new frequency grid, on which we compute the PSD
using the optimal number of averaged segments for each given frequency.
Following this method, we essentially split the time-series data at the
maximum number of segments that the given choice of window and percentage of
data overlap permits, which will allow us to estimate the PSD at each
frequency bin with minimal variance. By carefully choosing the window function
and distance between the data points, one can compute a power spectrum with
minimal correlations between frequencies. The estimated spectrum, $D_{i}$, at
each frequency, $f_{i}$, is then used in the likelihood function given in Eq.
(35) below (the spectrum $D_{i}$ is represented by the red data-points in
figure 4). For more details about this method of computing the PSD, we refer
the reader to (Armano et al., 2018; Tröbs & Heinzel, 2006) and for a similar
application to the work of Baghi et al. (2023).

Finally, we also keep two knots fixed at the edges of the spectra, allowing
the sampler to estimate only their amplitude, while the rest of the knot
parameters (and their number) are left to be estimated from the data.

(a)

(b)

(c)

Figure 4:  Results for power spectra modelling with a shape agnostic model.
(a) The simulated data (gray), generated from the theoretical model (dashed
black line). The PSD computed on an equally-spaced logarithmic grid with the
method of (Tröbs & Heinzel, 2006; Armano et al., 2018), which is used for
inference, is represented with the red data points. The pink solid lines
represent models drawn randomly from the posterior chains. (b) The optimal
B-spline knots estimated by the dynamical parameter estimation procedure. As
shown from this histogram, the optimal interior knot count for this data
converges to six, corresponding to eight total knots including the two edge
knots. (c) Posterior slice for the knot parameters, $(\log f_{j,k},\,\log
S_{j,k})$, after stacking the MCMC chains across all model dimensions, $k$.
This illustrates where the model prefers to place spline knots, which clearly
corresponds to where the spectral density is changing most rapidly. It is also
evident from this figure that we essentially “scan” the true noise shape (pink
solid line), by placing knots across the frequency range (see text for more
details).

For the spline knot positions, $\\{\log f_{j,k}\\}$, and amplitudes, $\\{\log
S_{j,k}\\}$, we adopt uniform priors that cover the complete parameter space.
Here, the $j$ index corresponds to the knot number for the given model order
$k$. We also use a ladder of $10$ temperatures, with $10$ walkers each, while
maintaining the same settings for the adaptivity of the temperatures as in
section 3.3.1. Each walker is initialized at a random point on the parameter
space, after drawing the dimensionality $k$ of the model from
$k\sim\mathcal{U}[1,20]$. We adopt a Gaussian likelihood, with its logarithm
written as

$\log
p(D|\vec{\theta}_{k})\propto-\frac{1}{2}\sum_{i}n_{i}\left(\frac{D_{i}}{{\cal
N}_{i,k}(\vec{\theta}_{k})}+\log{\cal N}_{i,k}(\vec{\theta}_{k})\right),$ (35)

where $D_{i}$ is the PSD data value for the given frequency $f_{i}$, as
computed by the method presented in (Tröbs & Heinzel, 2006; Armano et al.,
2018), using $n_{i}$ averaged segments. The ${\cal N}_{i,k}(\vec{\theta}_{k})$
is the spline noise model of order $k$ evaluated at $f_{i}$, that depends on a
parameter set

$\vec{\theta}_{k}=\\{\log f_{1,k},\cdots,\log f_{k,k},\,\log S_{0},\cdots,\log
S_{k,k},\,\log S_{k+1}\\},$ (36)

in which the $\log S_{0}$ and $\log S_{k+1}$ parameters refer to the logarithm
of the PSD amplitude of the two fixed knots at the “edges” of the spectrum.
Those two parameters correspond to our zeroth model order ($k=0$), thus they
are always being explored by the walkers of Eryn.

The results are shown in figure 4. In particular, in figure 4(b) we show the
histogram of the recovered number of knots for the particular data-set. It is
clear that $8$ spline knots are preferred, two of them being fixed at the
edges of the spectrum, and the other six knots free to take any position in
the given frequency range. In figure 4(c) we show the 2D sliced posteriors for
the spline parameters, $\\{\log S_{j,k}\\}$ and $\\{\log f_{j,k}\\}$. In this
figure, we again stack all the MCMC samples across model orders. The true
spectrum is indicated by the orange solid line. There is an interesting
outcome of this toy investigation; while there is a preferable dimensionality
of the model, there is a weak constraint on the actual positions of the knots.
We find that the sampler is virtually “scanning” the PSD data, showing
slightly higher preference for locations between $-6$ and $-4$ in
$\log$-frequency, where the spectrum follows a more complicated shape.
Finally, in figure 4(a), the data (gray solid line and red data points), is
shown together with model evaluations drawn from the posterior samples (pink
solid lines).

##  4  Examples from Gravitational Wave astronomy

In recent years, we have witnessed the beginning of Gravitational Wave
Astronomy. Since the first detection (Abbott et al., 2016) dozens of waveform
signatures have been measured with the current network of observatories. At
the time of the writing of this paper, more than 90 events have been recorded
(Abbott et al., 2021a), the vast majority of them are black hole (BH) binary
mergers, with a few of them being binary neutron star (NS) and BH-NS mergers.
At the same time, detector networks are being improved (The LIGO Scientific
Collaboration, 2019; Abbott et al., 2020) and there are plans to expand them
with the addition of new observatories, such as the Einstein Telescope
(Maggiore et al., 2020; Punturo et al., 2010) or Cosmic Explorer (Evans et
al., 2021; Abbott et al., 2017). Those detectors will unlock the sky to larger
redshifts $z$, allowing access to a vast number of potential sources. In
addition, space missions, such as LISA (Amaro-Seoane et al., 2017, 2012), are
predicted to be signal-dominated observatories, with many types of sources
populating their data streams. In fact, we expect that source confusion will
be one of the primary challenges in future data analysis efforts in
gravitational wave astronomy. In a typical data-set, we expect an unknown
number of signals, originating from sources that generate waveforms with
different characteristics. Those range from the stellar-mass BH binaries now
frequently observed by ground-basd detectors, to the supermassive BH binaries,
extreme mass ratio inspirals, ultra compact Galactic binaries (UCB), and
stochastic GW signals from both astrophysical and possibly cosmological origin
(Amaro-Seoane et al., 2017, 2012; Auclair et al., 2023). For this final
example, we will focus on the LISA mission, and in particular on the case of
discriminating UCB signals.

###  4.1  Application to LISA data and the Ultra Compact Galactic Binaries

LISA is going to measure GW signals in the $\mathrm{mHz}$ regime, accessing
sources of all the aforementioned types. As already discussed, the most
numerous of them are going to be the UCBs, which will be almost monochromatic
in the LISA band. Out of the millions of sources, only
$\sim\mathcal{O}(10^{4})$ will be individually resolvable, and the rest will
generate a confusion signal. As a consequence, for the duration of the
mission, we will need to disentangle tens of thousands of sources which will
be overlapping in both time and frequency domains. This is no trivial task,
but various different strategies have already been proposed for analyzing such
challenging data-sets. For example, Gaussian Processes can be employed (Strub
et al., 2022), or Swarm Optimization techniques (Zhang et al., 2021), or
hybrid swarm-based algorithms (Bouffanais & Porter, 2016). Pipelines based on
MCMC methods have been tested extensively (Crowder & Cornish, 2007;
Littenberg, 2011; Littenberg et al., 2020; Littenberg & Cornish, 2023), and
have been demonstrated to be able to tackle complex cases where signals are
overlapping.

Here, we will focus on the same problem, employing Eryn to solve a down-scaled
version of the UCB challenge. It is down-scaled because we focus only on a
single narrow frequency band, containing several overlapping signals, in the
presence of instrumental noise only333No confusion signal from other
unresolved UCBs is considered in this investigation.. In addition, we focus
solely on demonstrating the capabilities of Eryn on dynamic parameter
estimation for UCB type sources and no other types of signals are contained in
the data (e.g., chirping signatures from supermassive BH binaries). At the
same time, we have access to the level of instrumental noise, which is shown
in both panels of figure 5. Searching for the UCB signals across the complete
LISA band requires a more elaborate implementation of this simplified
pipeline. This pipeline will be focusing on solving the complete second LISA
Data Challenge (LISA Data Challenges Working Group, 2022), and is going to be
presented in future work.

# |  $f_{\mathrm{gw}}~{}[\mathrm{mHz}]$ |  $\rho_{\mathrm{opt}}$
---|---|---
1 |  3.99780 |  9.98
2 |  3.99781 |  46.70
3 |  3.99784 |  4.55
4 |  3.99854 |  39.45
5 |  3.99873 |  13.02
6 |  3.99882 |  8.47
7 |  3.99919 |  10.88
8 |  3.99939 |  19.07
9 |  3.99964 |  20.00
10 |  3.99965 |  7.99
Table 1: The optimal SNR $\rho_{\mathrm{opt}}$ for each of the 10 injected
sources, computed for the given duration of the mission (see eq. (42)). The
dominant emission frequency $f_{\mathrm{gw}}$ is also given for reference.

We choose to work on the frequency segment between $3.997$ and
$4~{}\mathrm{mHz}$, which contains $10$ UCB objects, drawn directly from the
LDC2 catalogue (LISA Data Challenges Working Group, 2022). Those are shown in
the top panel of figure 5 which shows the power spectrum of the $A$ data
channel of LISA. We use the two noise-orthogonal $A$, $E$, and $T$ Time Delay
Interferometry variables (Tinto & Dhurandhar, 2005; Prince et al., 2002; Baghi
et al., 2021), which are linear combinations of the LISA relative frequency
TDI Michelson measurements $X$, $Y$, and $Z$ as:

$\displaystyle A$ $\displaystyle=\frac{1}{\sqrt{2}}(Z-X),\quad
E=\frac{1}{\sqrt{6}}(X-2Y+Z),$ (37) $\displaystyle T$
$\displaystyle=\frac{1}{\sqrt{3}}(X+Y+Z).$

In ideal conditions (equal noises across spacecrafts, and equal LISA arms),
the noise in $A$ and $E$ is independent, while the $T$ data stream can be used
as a signal-insensitive null channel, useful for instrument noise calibration.
Since we perform analysis on a noise-free injection, we will be neglecting the
$T$ channel altogether. We simulated the injection data for an observation
time of $\mathrm{T}_{\mathrm{obs}}=1~{}\mathrm{year}$.

The optimal SNR for each injected source, $\rho_{\mathrm{opt}}$, is given in
table 1. The $\rho_{\mathrm{opt}}$ quantity refers to the SNR of each source
in isolation, with respect only to the instrumental noise, and can be
calculated as

$\rho_{\mathrm{opt}}^{2}=\sum_{C}\left(h_{C}|h_{C}\right)_{C},$ (38)

with $C\in\\{A,E\\}$ the noise-orthogonal TDI channels of eq. (37), while the
$\left(\cdot|\cdot\right)$ notation represents the noise weighted inner
product expressed for two time series $a$ and $b$ as

$\left(a|b\right)=2\int\limits_{0}^{\infty}\mathrm{d}f\left[\tilde{a}^{\ast}(f)\tilde{b}(f)+\tilde{a}(f)\tilde{b}^{\ast}(f)\right]/\tilde{S}_{n}(f).$
(39)

The tilde represents the data in the Fourier frequency domain, and the
asterisk indicates complex conjugate. The $\tilde{S}_{n}(f)$ is the one-sided
PSD of the noise for a given TDI channel. Under our assumptions
$S_{n,A}(f)=S_{n,E}(f)$.

Figure 5:  Top panel: The simulated data used for demonstrating the
capabilities of Eryn in tackling a high-dimensional problem. A total of 10
Ultra Compact Binaries in the vicinity of our Galaxy emitting Gravitational
Wave signals at the $\mathrm{mHz}$ range. Here we plot the power spectral
density of the $A$ data channel of LISA. The catalogue of sources is taken
from the second LISA Data Challenge (LISA Data Challenges Working Group,
2022). Each signal is represented by a different colour. Bottom panel: The
same data-set, now comparing the injected signal against the solution yielded
by Eryn (see text for more details). We have plotted the shaded area by
sampling the joint posterior on model order $k$ and the corresponding
parameters.

For our investigation we chose to analyze noiseless data (no noise
realization), while in the likelihood we are using the PSD noise levels taken
from the LISA design studies (LISA Science Study Team, 2018). For the signals,
we utilize the fast frequency-domain UCB waveform model of Cornish &
Littenberg (2007). Then, the two polarizations of an emitting UCB can be
written as

$\displaystyle h_{+}(t)$ $\displaystyle=\frac{2\mathcal{M}}{D_{L}}\left(\pi
f_{\mathrm{gw}}(t)\right)^{2/3}\left(1+\cos^{2}\iota\right)\cos{\phi(t)},$
(40) $\displaystyle h_{\times}(t)$
$\displaystyle=-\frac{4\mathcal{M}}{D_{L}}\left(\pi
f_{\mathrm{gw}}(t)\right)^{2/3}\cos{\iota}\sin{\phi(t)},$

where $\mathcal{M}$ is the chirp mass, $f_{\mathrm{gw}}$ is the instantaneous
gravitational wave frequency, $D_{L}$ is the luminosity distance, $\iota$ is
the inclination of the binary orbit, and $\phi(t)$ is the gravitational wave
phase. The phase $\phi$ can be expressed as
$\phi=\phi_{0}+2\pi\int^{t}f_{\mathrm{gw}}(t^{\prime})dt^{\prime}$, with
$\phi_{0}$ being an initial arbitrary phase shift. The LISA constellation is
assumed to be rigid and with equal arms, while the spacecrafts are assumed to
follow analytic Keplerian orbits (Cornish & Rubbo, 2003). Under these
assumptions, it is straightforward to compute its response to the almost
monochromatic waveforms of eqs. (40) (Cornish & Rubbo, 2003; Babak et al.,
2021). For more details about the waveform model, the response of the
instrument, and the orbits of the constellation, we refer the reader to
Cornish & Littenberg (2007); Robson et al. (2018); Katz et al. (2022); Cornish
& Rubbo (2003); Babak et al. (2021).

In our simplified scenario, each binary signal in the Solar System Barycenter
is determined by a set of eight parameters. Those are the
$\vec{\theta}=\\{\mathcal{A},\,f_{\mathrm{gw}}~{}[\mathrm{mHz}],\,\dot{f}_{\mathrm{gw}}~{}[\mathrm{Hz/s}],\,\phi_{0},\,\cos\iota,\psi,\lambda,\sin\beta\\}$,
where $\mathcal{A}$ is the overall amplitude, $\dot{f}_{\mathrm{gw}}$ is the
first derivative of the gravitational-wave frequency, $\psi$ the polarization,
$\lambda$ is the ecliptic longitude, and $\beta$ the ecliptic latitude of the
binary. The amplitude of the signal is calculated as

$\mathcal{A}=\left(2\mathcal{M}^{5/3}\pi^{2/3}f_{\mathrm{gw}}^{2/3}\right)/D_{L},$
(41)

which can be used to obtain a rough SNR estimate, via (Littenberg et al.,
2020)

$\rho^{2}=\frac{\mathcal{A}^{2}\mathrm{T}_{\mathrm{obs}}\sin^{2}(f_{\mathrm{gw}}/f_{\ast})}{4S_{n}(f_{\mathrm{gw}})},$
(42)

with $S_{n}(f_{\mathrm{gw}})$ being the instrumental noise power spectral
density at frequency $f_{\mathrm{gw}}$, and $f_{\ast}=1/(2\pi L)$, where $L$
the LISA arm length. Given eq. (41) and (42), we find it convenient to
directly sample on $\rho$ instead of $\mathcal{A}$, which also yields a more
illustrative measure of the amplitude of each binary. Then, if $d$ is the
measured TDI data and $\mathrm{h}$ the given GW signal after applying the
response of the instrument, the logarithm of the likelihood for an arbitrary
number $k$ of UCB signals can be written as

$\log
p(d|\vec{\theta}_{k})\propto\left(d|\mathrm{h}_{k}\right)-\frac{1}{2}\left(\mathrm{h}_{k}|\mathrm{h}_{k}\right),$
(43)

where for the sake of convenience, we have defined
$\mathrm{h}_{k}=\sum_{k}\mathrm{h}(\vec{\theta}_{k})$.

We use wide uniform priors for all the rest of the binary parameters, covering
essentially the complete parameter space. The exception is again the amplitude
(SNR), $\rho$, where we adopt a prior which was first introduced in Cornish &
Littenberg (2015) and then adapted in Littenberg et al. (2020). The prior
density can be expressed as

$p(\rho)=\frac{3\rho}{4\rho_{\ast}\left(1+\rho/(4\rho_{\ast})\right)},$ (44)

where $\rho_{\ast}$ is a given constant that specifies the peak of the above
distribution. This distribution is designed to prevent the proposal of sources
with very small SNR in the model, as it drops sharply for $\rho\rightarrow 0$.
Those weak sources do not significantly affect the likelihood, and so their
inclusion must be penalised by the prior 444We remind that the prior is also
used as our “birth” proposal here.. This prior choice forces the sampler to
explore only potential sources with non-zero SNR, avoiding populating the
chains with numerous undetectable signals. This prior performs adequately in
this problem, but there are other solutions one could adopt in order to keep
control of the number of very weak sources. This discussion, which sets the
grounds for a global-fit analysis pipeline for the LISA data (Littenberg et
al., 2020), is out of the scope of this paper, but a more detailed description
will be presented in a future work.

(a)

(b)

Figure 6: Left panel: In this figure, we show the posterior on the number of
UCB sources in the data. The true injected number is shown with the red dashed
line. It is clear that, for the given measurement duration of the particular
data set, we manage to confidently resolve eight binaries out of a total of
ten. Right panel: Corner plot for two of the eight parameters characterising
each UCB source. These are the amplitude, expressed as an SNR $\rho$, and the
dominant emission (or initial) frequency, $f_{\mathrm{gw}}~{}[\mathrm{mHz}]$
(see text for more details). The violet crosses represent the injected
parameter values. A corner plot for more parameters is shown in figure A.7 in
the Appendix.

##### Search Phase:

Before parameter estimation, we initiate a search phase of our analysis, with
the aim of getting the walkers to a better starting point on the posterior
surface. This phase consists of an iterative brute force procedure, based on
drawing a very large number of proposals, then maximizing the likelihood over
the initial phase $\phi_{0}$, and finally perform a rapid MCMC sampling over
the parameter space, using only a one-source model (therefore there is no
dynamical parameter spaces). In particular, we draw $5\times 10^{6}$ points in
the parameter space, and after phase maximization, we use them as starting
points to a parallel tempered MCMC run with $N_{T}=10$ temperatures, each
running with $n_{w}=500$ walkers. When this step concludes, we keep the 100
best samples in terms of likelihood value and use their corresponding
parameter estimates as starting points for the parameter estimation portion of
the analysis. We also use the maximum Likelihood solution to subtract the
source found from the data. We then use the residual data to search for
another source, and this process repeats until there is no signal found with
SNR $\rho>5$. In between successive iterations of the single-binary search, we
run another MCMC over all sources found so far in order to adjust the
parameters to account for correlations and overlap between sources. After
convergence, we found eight sources in our data-set with an optimal SNR
greater than 5. We take these found sources and add them to all walkers in the
sampler at the beginning of the full MCMC run described below.

##### Parameter estimation:

During this step, we perform hybrid MCMC sampling, where we both update the
found sources (in-model) and dynamically search for new and weaker sources in
the data employing Reversible Jump sampling. For the number of signals $k$, we
adopt a uniform prior $k\sim\mathcal{U}[6,\,20]$. For the sake of convenience
in this simple application, we keep the six loudest binaries found during the
search phase as fixed. This means that we still sample their waveform
parameters, but they are not allowed to be removed by the Reversible Jump
process. We chose this setup in order to accelerate the convergence of the
algorithm, being confident that these sources are part of the solution. In
future work, this will be adjusted to deal with the much larger complexity of
the full problem.

Concerning the sampler settings, we use the adaptive parallel tempering scheme
of section 2.4, building a temperature ladder of $N_{T}=10$ temperatures, with
$100$ walkers for each temperature. For this run, we have also utilized the
Multiple Try Metropolis algorithm (see section 2.3) in order to improve the
acceptance rate in the reversible jump proposal. We have also tried the
Delayed-Rejection scheme which is implemented for Eryn, but we found that the
Multiple-Try strategy yields more efficient sampling. Finally, we have
utilized the basic stretch and group stretch proposals that were described in
section 2.1.

After convergence, the result is shown in figures 5 and 6. In figure 6(a), the
sampled posterior on the number of sources $k$ is presented. In this
histogram, we have added the six fixed binaries to the actual number of
signals being sampled via the Reversible Jump algorithm. It is fairly obvious
that we have managed to confidently resolve eight out of the ten injected
binary signals. This fact that we do not favour $10$ sources can be explained
partly by the low SNR of the signals (see table 1) and partly by confusion
from source overlap (also shown in figure 5). Additionally, the result of
figure 6(a) depends on the given observation duration. The greater the
$T_{\mathrm{obs}}$, the better our ability to resolve the confused sources.
Thus, in that case, we should expect more Reversible Jump iterations across
the higher dimensional models.

On the right panel, in figure 6(b), the ensemble 2D posterior slice is shown,
for two selected parameters. We call it ensemble because we are again
“stacking” all the chains for these two parameters for all sources for all
model orders $k$. We chose to show only the amplitude (the $\rho$ parameter
explained in eq. (42) above) and the dominant emission frequency
$f_{\mathrm{gw}}$, which illustrates the number of sources resolved, and how
they overlap in frequency. A corner plot for more parameters is shown in
figure A.7 in Appendix A. We also show the true injection values, marked as
crosses, on top of the 2D posterior. From this plot alone, one can see that
the sampler is exploring efficiently the parameter space, converging to the
true values of the resolvable binaries that were injected.

##  5 Discussion

We have implemented Eryn, a Bayesian sampling package capable of performing
efficient trans-dimensional inference, by employing different techniques that
improve its acceptance rate. These techniques are the affine invariant
sampling, the adaptive parallel tempering, the delayed rejection, and multiple
try metropolis, in combination with the construction of informative proposal
distributions for the parameters of the models. The structure of Eryn is based
on the widely used software emcee (Foreman-Mackey et al., 2013), enhanced with
the ability of performing Reversible Jumps Green (1995) between different
model spaces. The sampler capabilities have been demonstrated with toy models
that are commonly encountered in different data analysis problems. We have
begun with an application to signal detection, and in particular to searching
for simple signals in the form of Gaussian pulse signals in the presence of
Gaussian noise (see section 3.3.1).

In section 3.3.2, we applied our algorithm to a problem of modeling power
spectra with arbitrary shapes in frequency domain. In such cases, it is
convenient to define models based on B-splines, which are able to faithfully
capture the shape of any spectral data. However, in order to avoid over-
fitting situations, the optimal order of the model (i.e. the optimal number of
spline knots), needs to be estimated from the data. This can be done either
sequentially, by trying models of different dimensionality and then comparing
their performance, or dynamically, by using trans-dimensional algorithms such
as Eryn. This class of problems is often encountered in cosmology (Planck
Collaboration et al., 2020a, b), where the signal of interest is stochastic in
nature, and sometimes the prior knowledge on its shape is very limited. As
already discussed, this is especially true for future GW observatories, which
open the possibility of detecting such signals from both astrophysical and
cosmological origin (Amaro-Seoane et al., 2017; Auclair et al., 2023; Baghi et
al., 2023). The different theoretical models produce spectra with distinct
shapes, increasing the need for shape-agnostic spectral models, such as the
B-spline used here.

Finally, in section 4.1, we demonstrated Eryn in a more complicated problem,
that of the analysis of ultra compact binary signals measured by the future
LISA detector. These objects are going to produce the majority of the signals
in the LISA data, each emitting almost monochromatic radiation. Their vast
number will generate a confusion foreground, while only a few thousand of them
will be resolvable from the data. We employ our tools described in this work,
together with a search phase that is based on iteratively running the sampler
on “static models” (no trans-dimensional moves) with phase-maximized
likelihoods. We do these runs on the residuals of each iteration, with the aim
of extracting all bright sources. In order to achieve faster convergence of
our parameter estimation run, we choose to keep the brightest sources found
during search as fixed (minimum number of model order is $k=6$), while the
Reversible Jump algorithm is used to search for weaker signals in the data.
This is purely a choice that allows quick convergence in this fully-controlled
and simplified LISA data-set.

We perform this analysis for a mission duration of
$T_{\mathrm{obs}}=1~{}\mathrm{year}$ and only on a single narrow frequency
band around $4~{}\mathrm{mHz}$, which contains a total of $10$ binary signals.
It is worth noting here, that the synthetic data were produced assuming
idealized conditions. This means that we do not consider any data
irregularities, such as data gaps and glitches and spectral lines, or any
other contamination originating from the mixing of signals of different types
(such as supermassive BH binaries). In the end, as shown in figure 6(a), we
manage to recover eight out of ten injected signals. This result makes sense
given the relative strength of the injections, and their waveform overlap.
Many of the injected sources have an optimal SNR in isolation which is rather
low (see table 1), so these are more susceptible to deterioration when we
account for signal overlap.

The above investigations demonstrate that the dynamical parameter estimation
capabilities of Eryn are suitable for these types of problems. This feature is
missing from already existing libraries such as bilby Ashton & Talbot (2021),
which are used by the Gravitational Wave community. bilby offers a wide
selection of tools which are necessary for the Bayesian analysis of
Gravitational Wave data. These include implementations of likelihood and prior
functions, instrument response models and spectral densities, waveform models,
and a number of samplers to choose from (both MCMC and nested samplers). On
the other hand, Eryn is based on parallel tempering MCMC enhanced with
Reversible Jump, which allows Bayesian analyses for a wider set of problems,
in which sampling of a dynamical parameter space is needed.

Eryn has already been used in several works that have been already published
(Katz et al., 2022; Katz, 2022; Baghi et al., 2023; Sasli et al., 2023), or
are going to appear soon. The work presented in this paper is the initial part
of our efforts toward implementing a data analysis pipeline for LISA data.
This pipeline will be demonstrated on the LDC2 data-set (LISA Data Challenges
Working Group, 2022), which contains multiple types of signals overlapping in
both time and frequency domains. That being said, Eryn is a generic and
versatile sampler, which can be used in any investigation that requires
Reversible Jump sampling, and to our knowledge is one of the very few
statistical tools of this kind that is not specialized to a single type of
analysis (see discussion in section 2.6).

## Acknowledgements

We wish to thank S. Babak, M. Le Jeune, S. Marsat, T. Littenberg, and N.
Cornish for their useful comments and very fruitful discussions. N Stergioulas
and N Karnesis acknowledge support from the Gr-PRODEX 2019 funding program
(PEA 4000132310).

## Data Availability

The sampler code can be found at https://github.com/mikekatz04/Eryn. The code
for the toy examples of section 3.3 can be found in the documentation of Eryn
(https://mikekatz04.github.io/Eryn/html/index.html#). For the analysis of
section 4, we used the GBGPU package (https://github.com/mikekatz04/GBGPU),
which contains an implementation of the UCB waveforms of the LDC software
(https://lisa-ldc.lal.in2p3.fr/code).

## References

  * Abbott et al. (2016) Abbott B. P., et al., 2016, Phys. Rev. Lett., 116, 061102
  * Abbott et al. (2017) Abbott B. P., et al., 2017, Classical and Quantum Gravity, 34, 044001
  * Abbott et al. (2020) Abbott B. P., et al., 2020, Living reviews in relativity, 23, 1
  * Abbott et al. (2021a) Abbott R., et al., 2021a, preprint (arXiv:2111.03606)
  * Abbott et al. (2021b) Abbott R., et al., 2021b, SoftwareX, 13, 100658
  * Abbott et al. (2021c) Abbott R., et al., 2021c, Phys. Rev. D, 104, 022004
  * Amaro-Seoane et al. (2012) Amaro-Seoane P., et al., 2012, Classical and Quantum Gravity, 29, 124016
  * Amaro-Seoane et al. (2017) Amaro-Seoane P., et al., 2017, preprint,  (arXiv:1702.00786)
  * Amaro-Seoane et al. (2022) Amaro-Seoane P., et al., 2022, preprint (arXiv:2203.06016)
  * Andrieu & Thoms (2008) Andrieu C., Thoms J., 2008, Statistics and Computing, 18, issue 4, 343
  * Armano et al. (2018) Armano M., et al., 2018, Phys. Rev. Lett., 120, 061101
  * Ashton & Prix (2018) Ashton G., Prix R., 2018, Phys. Rev. D, 97, 103020
  * Ashton & Talbot (2021) Ashton G., Talbot C., 2021, Mon. Not. Roy. Astron. Soc., 507, 2037
  * Auclair et al. (2023) Auclair P., et al., 2023, Living Rev. Rel., 26, 5
  * Babak et al. (2021) Babak S., Petiteau A., Hewitson M., 2021, preprint (arXiv:2108.01167)
  * Baghi et al. (2021) Baghi Q., Thorpe J. I., Slutsky J., Baker J., 2021, Phys. Rev. D, 103, 042006
  * Baghi et al. (2023) Baghi Q., Karnesis N., Bayle J.-B., Besançon M., Inchauspé H., 2023, J. Cosmology Astropart. Phys., 2023, 066
  * Baio (2012) Baio G., 2012, Bayesian Methods in Health Economics. Chapman & Hall/CRC Biostatistics Series, CRC Press, https://books.google.gr/books?id=2RzSBQAAQBAJ
  * Bédard et al. (2012) Bédard M., Douc R., Moulines E., 2012, Stochastic Processes and their Applications, 122, 758
  * Betancourt & Girolami (2013) Betancourt M. J., Girolami M., 2013, preprint, p. arXiv:1312.0906 (arXiv:1312.0906)
  * Biwer et al. (2019) Biwer C. M., Capano C. D., De S., Cabero M., Brown D. A., Nitz A. H., Raymond V., 2019, Publ. Astron. Soc. Pac., 131, 024503
  * Bouffanais & Porter (2016) Bouffanais Y., Porter E. K., 2016, Phys. Rev. D, 93, 064020
  * Breivik et al. (2020) Breivik K., et al., 2020, ApJ, 898, 71
  * Brooks (1998) Brooks S. P., 1998, Journal of the Royal Statistical Society. Series D (The Statistician), 47, 69
  * Brooks & Gelman (1998) Brooks S., Gelman A., 1998, Journal of Computational and Graphical Statistics, 7, 434
  * Christensen & Meyer (2022) Christensen N., Meyer R., 2022, Rev. Mod. Phys., 94, 025001
  * Collaboration & the Virgo Collaboration (2020) Collaboration T. L. S., the Virgo Collaboration 2020, Classical and Quantum Gravity, 37, 055002
  * Cornish (2021) Cornish N. J., 2021, Phys. Rev. D, 104, 104054
  * Cornish & Littenberg (2007) Cornish N. J., Littenberg T. B., 2007, Phys. Rev. D, 76, 083006
  * Cornish & Littenberg (2015) Cornish N. J., Littenberg T. B., 2015, Classical and Quantum Gravity, 32, 135012
  * Cornish & Rubbo (2003) Cornish N. J., Rubbo L. J., 2003, Phys. Rev. D, 67, 022001
  * Crowder & Cornish (2007) Crowder J., Cornish N. J., 2007, Phys. Rev. D, 75, 043008
  * De et al. (2018) De S., Finstad D., Lattimer J. M., Brown D. A., Berger E., Biwer C. M., 2018, Phys. Rev. Lett., 121, 091102
  * Dellaportas et al. (2002) Dellaportas P., Forster J. J., Ntzoufras I., 2002, Statistics and Computing, pp 27–36
  * Dinh et al. (2016) Dinh L., Sohl-Dickstein J., Bengio S., 2016, Density estimation using Real NVP, doi:10.48550/ARXIV.1605.08803, https://arxiv.org/abs/1605.08803
  * (36) Documentation Eryn., , Eryn Documentation , mikekatz04.github.io/Eryn/
  * Durkan et al. (2019) Durkan C., Bekasov A., Murray I., Papamakarios G., 2019, Neural Spline Flows, doi:10.48550/ARXIV.1906.04032, https://arxiv.org/abs/1906.04032
  * Evans et al. (2021) Evans M., et al., 2021, preprint, p. arXiv:2109.09882 (arXiv:2109.09882)
  * Foreman-Mackey et al. (2013) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, Publ. Astron. Soc. Pac., 125, 306
  * Frenkel (2006) Frenkel D., 2006, Waste-Recycling Monte Carlo. Springer Berlin Heidelberg, Berlin, Heidelberg, pp 127–137, doi:10.1007/3-540-35273-2_4, https://doi.org/10.1007/3-540-35273-2_4
  * Gelman & Rubin (1992) Gelman A., Rubin D. B., 1992, Statistical Science, 7, 457
  * Gelman et al. (1996) Gelman A., Roberts G. O., Gilks W. R., 1996, in Bernardo J. M., Berger J. O., Dawid A. P., Smith A. F. M., eds, , Bayesian Statistics. Oxford University Press, Oxford, pp 599–608
  * Gelman et al. (2004) Gelman A., Carlin J. B., Stern H. S., Rubin D. B., 2004, Bayesian Data Analysis, 2nd ed. edn. Chapman and Hall/CRC
  * Gilks & Wild (1992) Gilks W. R., Wild P., 1992, Journal of the Royal Statistical Society. Series C (Applied Statistics), 41, 337
  * Gilks et al. (1995) Gilks W., Richardson S., Spiegelhalter D., 1995, Markov Chain Monte Carlo in Practice. Chapman & Hall/CRC Interdisciplinary Statistics, Taylor & Francis, %****␣eryn.tex␣Line␣900␣****http://books.google.com/books?id=TRXrMWY_i2IC
  * Godsill (2001) Godsill S. J., 2001, Journal of Computational and Graphical Statistics, 10, 230
  * Goggans & Chi (2004) Goggans P. M., Chi Y., 2004, in Bayesian Inference and Maximum Entropy Methods in Science and Engineering. pp 59–66, doi:10.1063/1.1751356
  * Goodman & Weare (2010) Goodman J., Weare J., 2010, Communications in Applied Mathematics and Computational Science, 5, 65
  * Green (1995) Green P. J., 1995, Biometrika, 82, 711
  * Green & Mira (2001) Green P. J., Mira A., 2001, Biometrika, 88, 1035
  * Haario et al. (2001) Haario H., Saksman E., Tamminen J., 2001, Bernoulli, 7, 223
  * Haario et al. (2006) Haario H., Laine M., Mira A., Saksman E., 2006, Stat Comp, pp 339–354
  * Hastings (1970) Hastings W. K., 1970, Biometrika, 57, 97
  * Hermans et al. (2019) Hermans J., Begy V., Louppe G., 2019, Likelihood-free MCMC with Amortized Approximate Ratio Estimators, doi:10.48550/ARXIV.1903.04057, https://arxiv.org/abs/1903.04057
  * Hitchcock (2003) Hitchcock D. B., 2003, The American Statistician, 57, 254
  * Hobson et al. (2009) Hobson M. P., Jaffe A. H., Liddle A. R., Mukherjee P., Parkinson D., 2009, Bayesian Methods in Cosmology
  * Hoffman & Gelman (2011) Hoffman M. D., Gelman A., 2011, preprint,  (arXiv:1111.4246)
  * Hogg & Foreman-Mackey (2018) Hogg D. W., Foreman-Mackey D., 2018, The Astrophysical Journal Supplement Series, 236, 11
  * Hukushima & Nemoto (1996) Hukushima K., Nemoto K., 1996, Journal of the Physical Society of Japan, 65, 1604
  * Isi et al. (2022) Isi M., Farr W. M., Chatziioannou K., 2022, Phys. Rev. D, 106, 024048
  * Johannes & Polson (2009) Johannes M., Polson N., 2009, MCMC Methods for Financial Econometrics. https://www8.gsb.columbia.edu/researcharchive/articles/564
  * Johannes & Polson (2010) Johannes M., Polson N., 2010, in AïT-SAHALIA Y., HANSEN L. P., eds, Handbooks in Finance, Vol. 2, Handbook of Financial Econometrics: Applications. Elsevier, San Diego, pp 1–72, doi:https://doi.org/10.1016/B978-0-444-53548-1.50003-9, https://www.sciencedirect.com/science/article/pii/B9780444535481500039
  * Joseph (2019) Joseph A., 2019. SpringerBriefs in Physics. Springer (arXiv:1912.10997), doi:10.1007/978-3-030-46044-0
  * Karnesis et al. (2014) Karnesis N., et al., 2014, Phys. Rev. D, 89, 062001
  * Karnesis et al. (2021) Karnesis N., Babak S., Pieroni M., Cornish N., Littenberg T., 2021, Phys. Rev. D, 104, 043019
  * Kass & Raftery (1995) Kass R. E., Raftery A. E., 1995, Journal of the American Statistical Association, 90, 773
  * Katz (2022) Katz M. L., 2022, Phys. Rev. D, 105, 044055
  * Katz et al. (2022) Katz M. L., Danielski C., Karnesis N., Korol V., Tamanini N., Cornish N. J., Littenberg T. B., 2022, Mon. Not. Roy. Astron. Soc., 517, 697
  * Kawamura et al. (2008) Kawamura S., et al., 2008, Journal of Physics: Conference Series, 122, 012006
  * Kendall et al. (2005) Kendall W., Liang F., Wang J., 2005, Markov Chain Monte Carlo: Innovations and Applications. Institute for Mathematical Sciences lecture notes series, World Scientific, https://books.google.gr/books?id=cB-ZZX2HcWQC
  * Khan et al. (2005) Khan Z., Balch T., Dellaert F., 2005, IEEE Transactions on Pattern Analysis and Machine Intelligence, 27, 1805
  * LISA Data Challenges Working Group (2022) LISA Data Challenges Working Group 2022, LISA Data Challenges, https://lisa-ldc.lal.in2p3.fr/
  * LISA Science Study Team (2018) LISA Science Study Team 2018, Technical Report ESA-L3-EST-SCI-RS-001, LISA Science Requirements Document, https://www.cosmos.esa.int/documents/678316/1700384/SciRD.pdf. ESA, https://www.cosmos.esa.int/documents/678316/1700384/SciRD.pdf
  * Lartillot & Philippe (2006) Lartillot N., Philippe H., 2006, Systematic Biology, 55, 195
  * Littenberg (2011) Littenberg T. B., 2011, Phys. Rev. D, 84, 063009
  * Littenberg & Cornish (2009) Littenberg T. B., Cornish N. J., 2009, Phys. Rev. D, 80, 063007
  * Littenberg & Cornish (2010) Littenberg T. B., Cornish N. J., 2010, Phys. Rev. D, 82, 103007
  * Littenberg & Cornish (2023) Littenberg T. B., Cornish N. J., 2023, Phys. Rev. D, 107, 063004
  * Littenberg et al. (2020) Littenberg T. B., Cornish N. J., Lackeos K., Robson T., 2020, Phys. Rev. D, 101, 123021
  * Liu et al. (2000) Liu J. S., Liang F., Wong W. H., 2000, Journal of the American Statistical Association, 95, 121
  * Lopes & West (2004) Lopes H. F., West M., 2004, Statistica Sinica, 14, 41
  * Luo et al. (2016) Luo J., et al., 2016, Classical and Quantum Gravity, 33, 035010
  * Maggiore et al. (2020) Maggiore M., et al., 2020, JCAP, 03, 050
  * Marrs (1997) Marrs A., 1997, in Jordan M., Kearns M., Solla S., eds, Vol. 10, Advances in Neural Information Processing Systems. MIT Press, https://proceedings.neurips.cc/paper/1997/file/0ed9422357395a0d4879191c66f4faa2-Paper.pdf
  * Martino (2018) Martino L., 2018, Digital Signal Processing, 75, 134
  * Martino et al. (2011) Martino L., Del Olmo V. P., Read J., 2011, preprint, p. arXiv:1112.4048 (arXiv:1112.4048)
  * Maturana-Russel et al. (2019) Maturana-Russel P., Meyer R., Veitch J., Christensen N., 2019, Phys. Rev. D, 99, 084006
  * McMillan (2017) McMillan P. J., 2017, MNRAS, 465, 76
  * Menéndez et al. (1997) Menéndez M., Pardo J., Pardo L., Pardo M., 1997, Journal of the Franklin Institute, 334, 307
  * Metropolis et al. (1953) Metropolis N., Rosenbluth A. W., Rosenbluth M. N., Teller A. H., Teller E., 1953, J. Chem. Phys., 21, 1087
  * Mira (2001) Mira A., 2001, Metron - International Journal of Statistics, LIX, 231
  * Muller (1991) Muller P., 1991, Technical Report 91-09, A generic approach to posterior integration and Gibbs sampling. Purdue University
  * Neal (2011) Neal R., 2011, in , Handbook of Markov Chain Monte Carlo. pp 113–162, doi:10.1201/b10905
  * Nelemans et al. (2001) Nelemans G., Yungelson L. R., Portegies Zwart S. F., 2001, A&A, 375, 890
  * Piegl & Tiller (1996) Piegl L., Tiller W., 1996, The NURBS Book, second edn. Springer-Verlag, New York, NY, USA
  * Planck Collaboration et al. (2020a) Planck Collaboration Aghanim, N. et al., 2020a, A&A, 641, A1
  * Planck Collaboration et al. (2020b) Planck Collaboration Aghanim N., et al., 2020b, A&A, 641, A5
  * Prince et al. (2002) Prince T. A., Tinto M., Larson S. L., Armstrong J. W., 2002, Phys. Rev. D, 66, 122002
  * Punturo et al. (2010) Punturo M., et al., 2010, Classical and Quantum Gravity, 27, 194002
  * Ren et al. (2023) Ren Z., et al., 2023
  * Ritter & Tanner (1992) Ritter C., Tanner M., 1992, Journal of the American Statistical Association, 87, 861
  * Robert & Casella (2011) Robert C., Casella G., 2011, Statistical Science, 26, 102
  * Roberts & Rosenthal (2009) Roberts G. O., Rosenthal J. S., 2009, Journal of Computational and Graphical Statistics, 18, 349
  * Robson et al. (2018) Robson T., Cornish N. J., Tamanini N., Toonen S., 2018, Phys. Rev. D, 98, 064012
  * Sasli et al. (2023) Sasli A., Karnesis N., Stergioulas N., 2023, preprint (arXiv:2305.04709)
  * Sharma (2017) Sharma S., 2017, ARA&A, 55, 213
  * Shawhan (2003) Shawhan P., 2003, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 502, 396
  * Sivia & Skilling (2006) Sivia D., Skilling J., 2006, Data Analysis: A Bayesian Tutorial. Oxford science publications, OUP Oxford, https://books.google.gr/books?id=lYMSDAAAQBAJ
  * Sorensen & Gianola (2007) Sorensen D., Gianola D., 2007, Likelihood, Bayesian, and MCMC Methods in Quantitative Genetics. Statistics for Biology and Health, Springer New York, https://books.google.gr/books?id=c3ezZ3rK1HsC
  * Strub et al. (2022) Strub S. H., Ferraioli L., Schmelzbach C., Stähler S. C., Giardini D., 2022
  * Swendsen & Wang (1986) Swendsen R. H., Wang J.-S., 1986, Phys. Rev. Lett., 57, 2607
  * The LIGO Scientific Collaboration (2019) The LIGO Scientific Collaboration 2019, Gravitational wave astronomy with LIGO and similar detectors in the next decade, doi:10.48550/ARXIV.1904.03187, https://arxiv.org/abs/1904.03187
  * Thrane & Talbot (2019) Thrane E., Talbot C., 2019, Publications of the Astronomical Society of Australia, 36, e010
  * Tinto & Dhurandhar (2005) Tinto M., Dhurandhar S. V., 2005, Living Reviews in Relativity, 8, 4
  * (115) TolkienGateway, Tolkien Gateway, http://tolkiengateway.net/wiki/Eryn_Vorn
  * Trias et al. (2009) Trias M., Vecchio A., Veitch J., 2009
  * Tröbs & Heinzel (2006) Tröbs M., Heinzel G., 2006, Measurement, 39, 120
  * Vallisneri (2008) Vallisneri M., 2008, Phys. Rev. D, 77, 042001
  * Virtanen et al. (2020) Virtanen P., Gommers R., Oliphant T. E., Haberland M., Reddy T., et al., 2020, Nature Methods, 17, 261
  * Vousden et al. (2016) Vousden W., Farr W., Mandel I., 2016, MNRAS, 455, 1919
  * Wong et al. (2020) Wong K. C., et al., 2020, MNRAS, 498, 1420
  * Wu et al. (2020) Wu H., Köhler J., Noé F., 2020, Stochastic Normalizing Flows, doi:10.48550/ARXIV.2002.06707, https://arxiv.org/abs/2002.06707
  * Xie et al. (2010) Xie W., Lewis P. O., Fan Y., Kuo L., Chen M.-H., 2010, Systematic Biology, 60, 150
  * Yu et al. (2021) Yu X., Wang X., Lu C., Wang Y., 2021, IOP Conference Series: Earth and Environmental Science, 660, 012124
  * Zhang et al. (2021) Zhang X.-H., et al., 2021, Phys. Rev. D, 104, 024023

## Appendix A

In figure A.7 we show the triangle plot of the stacked posterior points as
sampled by Eryn, for the investigation of section 4. The difference to figure
6(b) is that here we include more parameters of the sources, but we still do
not include all parameters for the sake of clarity.

Figure A.7:  A triangle plot showing the 2D posterior slices for the
application of section 4.1, but for a greater selection of parameters than
figure 6(b). The rest of the parameters, if plotted stacked in the same
manner, result in surfaces that cannot be so easily interpreted, and therefore
have been left out. The true injected values are marked with crosses.