# Modelling the Impact of Scandals: the case of the 2017 French Presidential
Election

Yassine Bouachrine and Carole Adam

(This is an ENSIMAG internship report originally written in June 2019 by
intern Yassine Bouachrine under the supervision of Carole Adam)

###### Abstract

This paper proposes an agent-based simulation of a presidential election,
inspired by the French 2017 presidential election. The simulation is based on
data extracted from polls, media coverage, and Twitter. The main contribution
is to consider the impact of scandals and media bashing on the result of the
election. In particular, it is shown that scandals can lead to higher
abstention at the election, as voters have no relevant candidate left to vote
for. The simulation is implemented in Unity 3D and is available to play
online. Keywords: agent-based simulation, computational social choice, voting
models

##  1 Introduction

During the 2017 French presidential election, the media had a very impactful
role in the shift of the opinion away from the election’s favorite François
Fillon. The seriousness of the accusations against the candidate led to Fillon
plummeting in the polls. We will try to model the impact of both conventional
and social media through scandal diffusion, in order to better understand the
dynamics underlying the voting process.

There are a variety of existing models for the voting process. However, most
of the models see voting as a discrete event, while it is actually a sample at
a given instant of a system in continuous evolution. This paper tries to shed
some light on improvements which could take into account the inherent dynamism
that comes with the interactions of the voters.

Another issue is that most computer voting models are adapted to the American
context. In France, there are multiple candidates participating in the first
round of the elections but in reality, few of them are actually considered by
the voters, despite being aligned with their ideals. This is due to the nature
of the voting process in France which takes place in two rounds. This leads
the voters to cast their vote strategically towards candidates who have a
chance of making it out of the first round.

The Voter Autrement111Voter Autrement : https://vote.imag.fr/ experiment
explored the effects of this strategic voting in 2017 by testing various
voting systems during the presidential election. The results showed that the
alternative voting methods yield vastly different results (in terms of who is
elected), especially for candidates such as J-L. Mélenchon and H. Hamon
(systematic improvement) or F. Fillon and M. Le Pen (systematic decline).
Other voting methods such as candidate ranking make the strategic vote
useless, since every candidate gets a chance to make it to the second turn.

Finally, the context for the 2017 presidential election is even more unique in
regard to the state of the French political scene. It is the first
presidential election since the split of the centre-right party (UMP), and it
is taking place amid the overall dissatisfaction of the population with the
French Socialist Party (PS, left wing). This will be further elaborated on in
the Voting models section as it is relevant to studying the impact of such a
context on the existing models.

Our goal here is to model the impact of scandals over the course of the
election and the change in the opinion of the voters. Last year, A. Soutif
[13] tried to model the election process using an agent-based model by feeding
the agents the results of the polls reported by the media. His goal in that
study was to model the impact of the polls on the votes and on the strategic
vote in particular, as polls give information on which candidates have chances
of making it out of the first round. Our approach aims to complement this work
by showing the additional impact of media through the diffusion of scandals
during the campaign.

The first part of the paper (Section 2) focuses on the data analysis upon
which the model is built. The second part of the paper (Section 3) addresses
voting models and our implementation of the suggested model.

##  2 Data analysis

Building a model requires the availability of sufficient and relevant data
about the phenomenon we are trying to model.

###  2.1 Comparing poll results and media trends

In our case, we took the aggregated poll results from various organisms [14]
shown in Figure 1, and interpolated linearly when we had missing data (mostly
B. Hamon in the early polls). We then compared the result of these polls with
the evolution of the presence of the candidates across traditional media and
Twitter.

Figure 1: Evolution of the polls during the elections (Source: Wikipedia [15])

Google Trends222https://trends.google.fr/trends/ is a Google website that
offers statistics on Google Search queries. We used it to evaluate the
searched queries associated with the the top five candidates in the News
section, the results are shown on Figure 2. An interesting point can be made
around how F. Fillon dominates the media presence in the early months of the
election, alongside B. Hamon which we can ignore since it is mostly due to the
French Socialist Party presidential primary.

Figure 2: Evolution of the search queries for the candidates

We can see that media coverage of this candidate spikes at multiple points in
time (Figure 3), which is likely due to the “Penelope gate”, which is a
scandal associated with the alleged fictitious employment of members of
Fillon’s family. Further, the media frenzy surrounding the “Penelope gate”
looks correlated with the evolution of the polls, as shown in Figure 3.

Figure 3: Evolution of the polls and search queries for Fillon

The takeaway from this comparison is that further media bashing does not seem
to impact polls as much. This can be attributed to the nature of the scandal
being a punctual event, and that people willing to cast away their vote
already did so with the initial media outlets. The biggest beneficiary of this
is definitely the closest candidate in the political spectrum: E. Macron.

Another interesting candidate is M. Le Pen who dropped significantly in the
polls (Figure 4) once the media picked up that she could actually finish first
at the first round.

Figure 4: Evolution of the polls and search queries for Le Pen

An hypothesis could be that there is actually a bidirectional relationship
between polls and media coverage, with each one affecting the other. For
example, in F. Fillon’s case it is the disclosure of the scandal by Le Canard
Enchaîné that led to him dropping significantly, whereas in M. Le Pen’s case,
it is the poll results that led to an increase in coverage, and then her poll
results dropped as a result of this increased coverage.

###  2.2 Processing Twitter data

The limits of the information we can leverage from Google Trends is that it
does not tell us about the nature or content of the coverage. We cannot know
for sure if the increasing number of articles are rather positive or negative
ones. We can deduce that a posteriori by observing the impact on the polls,
but it is what we are actually trying to model. Therefore, we used Twitter
data, and analyze the tweets during the primary round to try and visualize
opinion trends during the election.

####  2.2.1 Data sets

Two datasets were used in our work:

  * •

The first one is from Kaggle 333https://www.kaggle.com/jeanmidev/french-
presidential-election: Kaggle is an online exchange platform for
datascientist, users can publish datasets among other things. This Kaggle
dataset contains tweets sampled during the elections. It is very rich but the
data collection rate varies and some tweets appear to be truncated. More
details about this dataset can be found at the source.

  * •

The other dataset is a courtesy of E. Duble, research engineer at LIG. It
contains an anonymized collection, sampling only geotagged tweets, that are
mentioning the top hashtags during the election. The importance of hashtags
has been shown for instance by the Politoscope project [10]

####  2.2.2 Clustering

Building a vectorial representation of the tweets can be done in various ways.
M. Campr and K. Jezek [3] provided a performance evaluation of various methods
for paragraph vectorization. At first, we opted for Tweet2Vec [6], which
relies on character-based representations (as opposed to word representations
for Doc2Vec [8]) that perform better for content such as tweets. However,
Tweet2vec uses hashtag prediction to train the model, which is limited for our
use-case since we already have a restricted number of hashtags. It also takes
longer to train compared to word-based models.

We used Facebook fastText 444Facebook FastText: https://fasttext.cc/ to
generate embeddings, and enhanced them with their term frequency-inverse
document frequency (TF-IDF, [12]) in the corpus. We then computed the tweet
embeddings as the average of the word embeddings.

Performing Principal Component Analysis (PCA, [16]) on the 100-dimensional
tweet embeddings did not yield very good results as the embeddings are already
built to minimize colinearity.

##  3 Voting model

In this section, we provide an overview of existing models and their
limitations before presenting the model we built through the observation of
the data.

###  3.1 Existing models and their limitations

Doing a taxonomy of existing models is outside the scope of this paper, there
are good resources available for that [1, 9, 7]. We will focus on some
limitations of the existing models, starting with the psychological model. As
hinted to in the introduction, the context of the election makes partisan
identification hard to rely on: the schism of the centre-right party reshaped
the political scene entirely. Also, the overall dissatisfaction with F.
Hollande hurt the socialist party. Towards the end of his mandate, his
popularity rating was lower than Macron’s was during the Gilet Jaunes protests
[11]. On top of that, the appearance of new actors on the scene such as En
Marche further shook the scene. En Marche made retrospective voting irrelevant
as the party had never held responsibilities. The novelty Macron brings to the
table, and his ambition of uniting the political parties, gave him a
considerable advantage.

All of these circumstances made the elections very volatile. Even more
sophisticated models such as the funnel of causality have to be rethought.

Figure 5: Funnel of causality, source [4]

Media has to have a bigger role in this funnel, especially social media as it
has been shown to be a good indicator of standings [2], almost as good as
traditional polls. And that is, despite it being sensitive to social
engineering (cf. Cambridge Analytica’s impact on the American presidential
election).

###  3.2 Proposed model

Our model is an enhancement of proximity models, in order to take into account
the diffusion of scandals and the movement of neighboring agents.

####  3.2.1 Simulation environment and initialization

The environment is a 100 X 100 units 2D plan. There are two types of active
agents: candidates and voters.

For the simulation, we define:

  * •

The appeasement delta $\Delta\alpha\in[0,1]$, rate at which the repulsion of
the candidate diminishes

  * •

The falloff rate for the potential of the scandals $\Delta\rho\in[0,1]$

  * •

The maximum openness for the voters $\sigma_{max}\in[0,100]$, defines how far
a voter considers his surroundings

  * •

The maximum tolerance for the voters $\theta_{max}\in[0,+\infty[$

####  3.2.2 Agents and their attributes

For candidates $C$, we define:

  * •

Position at time $t$ as $\psi_{t}\in[0,1]^{2}$ initialized manually

  * •

Repulsion at time $t$ as $\gamma_{t}\in[0,1]$ with $\gamma_{0}=0$

  * •

A list S of scandals with $S_{i}$ being the i-th one.

For voters $V$, we define:

  * •

Position at time $t$ as $\psi_{t}\in[0,1]^{2}$ with
$\psi_{0}\sim\mathcal{U}^{2}(0,1)$ (distance to a candidate inversely
proportional to agreement with this candidate)

  * •

Openness as $\sigma\in[0,\sigma_{max}]$ with $\sigma\sim
clamp_{[0,1]}(\mathcal{N}(0.5,0.2^{2}))\sigma_{max}$ (the radius is which a
voter considers agents around him)

  * •

Charisma as $\kappa\in[0,1]$ with $\kappa\sim
clamp_{[0,1]}(\mathcal{N}(0.5,0.2^{2}))$ (the chance to influence others
around)

  * •

Tolerance as $\theta\in[0,1]$ with $\theta\sim
clamp_{[0,1]}(\mathcal{N}(0.5,0.2^{2}))\theta_{max}$ (the threshold for
repulsion before dismissing a candidate completely)

  * •

Conformity as $\eta\in[0,1]$ with $\eta\sim
clamp_{[0,1]}(\mathcal{N}(0.5,0.2^{2}))$

For scandals, we define:

  * •

Potential at time $t$ as $\rho_{t}\in[0,1]$ with $\rho_{0}$ initialized
manually by the user (parameter in the simulation)

####  3.2.3 Simulation update

At each time-step of the simulation, we update the entities. To simplify the
equations, we assume that the values are clamped to their domain.

For scandals, the potential decreases with time, at the falloff rate:

$\rho_{t+1}(x)=\rho_{t}(x)-\Delta\rho$ (1)

For candidates, the position is static, and the repulsion increases with each
scandal, and decreases with time at the pace set by the appeasement delta:

$\gamma_{t+1}(x)=\gamma_{t}(x)-\Delta\alpha+\sum_{y\in S(x)}\rho_{t+1}(y)$ (2)

For voters, the position evolves as their opinion about the different
candidates evolves based on their surroundings:

$\begin{split}\psi_{t+1}(x)&=\psi_{t}(x)\\\
&+\operatorname*{argmin}_{y}\\{||\psi_{t}(x)-\psi_{t}(y)||_{2}\mid y\in
C\land\gamma_{t+1}(y)<\theta(x)\\}\frac{1}{1+\gamma_{t+1}(y)}\\\
&+\eta\sum_{y\in
V\land||\psi_{t}(x)-\psi_{t}(y)||_{2}<\sigma(x)}\kappa(y)(\psi_{t}(x)-\psi_{t}(y))\end{split}$
(3)

When the simulation stops, each voter votes for the closest candidate still
considered in the openness radius around him. If there are none, the voter
withholds his vote (chooses abstention).

####  3.2.4 Implementation details

The simulation is built in Unity 2018.3.5f1 555Unity 3D: https://unity.com/.
It allowed for faster prototyping and also supports a wide range of platforms
to run the simulation on. The simulation is available to play online at
http://lig-tdcge.imag.fr/votsim/ or to download as a WebGL export666WebGL
export of the simulator: https://ensiwiki.ensimag.fr/index.php?title=IRL_-
_Modélisation_de_la_dynamique_des_opinions_des_électeurs.

A first screen lets the user select the values of the global parameters of the
simulation (Figure 6): the number of voters and candidates, the appeasement
delta (rate at which the scandals decrease, which determines the duration of
their effect on opinions), and the maximum values for tolerance and openness
(individual values of all agents are then set randomly under this boundary).

Figure 6: Parameter selection before starting the simulation

On the next screen, the user can modify the simulation speed. Voters are
moving in the environment towards or away from the candidates. The user can
also trigger a scandal and choose its intensity and target candidate, by using
a button at the bottom of the window, in order to observe the influence on the
movements of the voters (see Figure 7). The intensity of the scandal then
decreases with time.

Figure 7: Screenshot of the simulation after triggering a scandal

####  3.2.5 Discussion and results

Regarding the motivations behind the model, we wanted to enhance the existing
models with the observations made from our data analysis. First, regarding the
initialization, the justification behind the uniform distribution for the
voters’ positions is the unique context behind the 2017 election that we
detailed earlier, with many voters not knowing which parties to consider. A
fine-tuned Gaussian mixture model could also be explored.

Most of the reasoning behind the model is based on the reactions to the
Penelopegate, with some voters completely turning their backs on F. Fillon
(which we model by a scandal being above their tolerance threshold) and some
only showing hesitation (tolerance threshold not reached). A temporary dip
followed by a partial recovery in the polls supports the model: regardless of
the further media bashing around the event, voters have a threshold over which
further coverage has no effect.

Over the simulated scenarios, one of the most interesting observations is that
scandals tend to be tied with an increased abstention rate. In our model this
is represented by the voter moving too far from all candidates (due to
repulsion generated by scandals about their favourite candidates, or to
diverging opinions with the others), so that no valid candidate is still
present in the openness radius when the election comes; in that case the voter
prefers to choose abstention. Our model can therefore reproduce and explain
the 2017 presidential election’s high abstention rate in the first round, at
22.23% [5].

##  4 Conclusion

We have seen that media coverage of the campaign scandals can have a big
impact on the election results. The simulation showed that scandals can
totally shape the result of the election and that scandals profit to the
closer candidates on the political spectrum. The more interesting finding was
how scandals impact the abstention rate, which is in agreement with the
observations made in the context of the French 2017 presidential election and
the high abstention rate recorded.

There is still much to do to reach a unified model, a first step in that
direction would be enhancing the simulation with the results from A. Soutif’s
experiments regarding the impact of the polls [13]. We could then initialize
the model to match the French political scene at the beginning of the first
round and test if it corresponds to the observed election results. If the
model is validated, we could explore alternative scenarios for the election:
how different scandals could have led to different results and particularly
what would have happened if there were no scandals involving the pre-campaign
favorite F. Fillon.

## References

  * [1] Rui Antunes.  Theoretical models of voting behaviour.  Exedra, 4(1):145–70, 2010.
  * [2] David Anuta, Josh Churchin, and Jiebo Luo.  Election bias: Comparing polls and twitter in the 2016 us election.  arXiv preprint arXiv:1701.06232, 2017.
  * [3] Michal Campr and Karel Ježek.  Comparing semantic models for evaluating automatic document summarization.  In International Conference on Text, Speech, and Dialogue, pages 252–260. Springer, 2015.
  * [4] Russell J Dalton et al.  Citizen politics in Western democracies: Public opinion and political parties in the United States, Great Britain, West Germany, and France.  Chatham House Chatham, NJ, 1988.
  * [5] Laurent de Boissieu.  Participation et abstention aux élections.  France politique, 2016.
  * [6] Bhuwan Dhingra, Zhong Zhou, Dylan Fitzpatrick, Michael Muehl, and William W Cohen.  Tweet2vec: Character-based distributed representations for social media.  Technical report, Arxiv, 2016.  https://arxiv.org/abs/1605.03481.
  * [7] Marco Giugni.  Theoretical models of voting behaviour.  https://baripedia.org/wiki/Theoretical_models_of_voting_behaviour, 2011\.
  * [8] Jey Han Lau and Timothy Baldwin.  An empirical evaluation of doc2vec with practical insights into document embedding generation.  arXiv preprint arXiv:1607.05368, 2016.
  * [9] Michael S Lewis-Beck and Mary Stegmaier.  Economic models of voting.  In The Oxford handbook of political behavior. Oxford, 2007.
  * [10] Maziyar Panahi, Noe Gaumont, and David Chavalarias.  How political communities are using hashtags against their opponents in upcoming french presidential election: #presidentielle2017.  Technical report, ISC-PIF, Politoscope, 2017.  https://politoscope.org/2017/04/hashtags-matter-presidentielle2017/.
  * [11] Opinion Publique.  Popularité de François Hollande : Chiffres et analyses de la popularité de François Hollande et de son gouvernement.  Opinion Publique: Sondages publiés et analyses sur les élections, la politique et les sujets de société, 2014.  https://opinionpublique.wordpress.com/category/popularite-de-francois-hollande/.
  * [12] Juan Ramos et al.  Using tf-idf to determine word relevance in document queries.  In Proceedings of the first instructional conference on machine learning, volume 242.1, pages 29–48, 2003.
  * [13] Albin Soutif.  Simulation multi-agent de systèmes de vote.  Technical report, ENSIMAG. Projet d’Introduction à la Recherche en Laboratoire, encadré par Carole Adam et Sylvain Bouveret, 2017.  https://ensiwiki.ensimag.fr/index.php?title=Albin_Soutif_:_Simulation_multi-agent_de_syst%C3%A8mes_de_vote.
  * [14] Wikipedia.  Opinion polling for the 2017 french presidential election.  https://en.wikipedia.org/wiki/Opinion_polling_for_the_2017_French_presidential_election, 2017\.
  * [15] Wikipedia.  Évolution des intentions de vote au premier tour de l’élection.  https://commons.wikimedia.org/wiki/File:Evolution_des_intentions_de_vote_à_l’élection_présidentielle_2017.png, 2017\.
  * [16] Svante Wold, Kim Esbensen, and Paul Geladi.  Principal component analysis.  Chemometrics and intelligent laboratory systems, 2(1-3):37–52, 1987\.


#  bilby in space: Bayesian inference for transient gravitational-wave signals
observed with LISA

C. Hoy1, L. K. Nuttall1
1University of Portsmouth, Portsmouth, PO1 3FX, United Kingdom
<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The Laser Interferometer Space Antenna (LISA) is scheduled to launch in the
mid 2030s, and is expected to observe gravitational-wave candidates from
massive black-hole binary mergers, extreme mass-ratio inspirals, and more.
Accurately inferring the source properties from the observed gravitational-
wave signals is crucial to maximise the scientific return of the LISA mission.
bilby , the user-friendly Bayesian inference library, is regularly used for
performing gravitational-wave inference on data from existing ground-based
gravitational-wave detectors. Given that Bayesian inference with LISA includes
additional subtitles and complexities beyond it’s ground-based counterpart, in
this work we modify bilby to perform parameter estimation with LISA. We show
that full nested sampling can be performed to accurately infer the properties
of LISA sources from transient gravitational-wave signals in a) zero-noise and
b) idealized instrumental noise. By focusing on massive black-hole binary
mergers, we demonstrate that higher order multipole waveform models can be
used to analyse a year’s worth of simulated LISA data, and discuss the
computational cost and performance of full nested sampling compared with
techniques for optimising likelihood calculations, such as the heterodyned
likelihood.

###### keywords:

gravitational waves – black hole binaries – methods: data analysis

††pubyear: 2023††pagerange: bilby in space: Bayesian inference for transient
gravitational-wave signals observed with LISA–References

##  1 Introduction

Gravitational-wave (GW) signals are regularly observed with existing ground-
based GW detectors (Aasi et al., 2015; Acernese et al., 2014; Akutsu et al.,
2021). By focusing on the Hertz to kiloHertz range of the GW spectrum, ground-
based GW detectors have observed $\sim 100$ signals from neutron star and
stellar-mass black hole coalescences (Abbott et al., 2021a; Nitz et al., 2023;
Venumadhav et al., 2020; Zackay et al., 2019b, a; Mehta et al., 2023), with
many more expected during the ongoing fourth GW observing run. To probe lower
frequencies, GW detectors must transition to space, as unavoidable terrestrial
noise means that ground-based GW detectors will never be sensitive to the
milliHertz frequency range (Danzmann, 1996).

The Laser Interferometer Space Antenna (LISA) is designed to explore the
milliHertz region of the GW spectrum, and is scheduled to launch in the mid
2030s (Amaro-Seoane et al., 2017). This frequency range is rich in sources,
with massive black hole binaries (MBHBs) (Klein et al., 2016), extreme mass
ratio inspirals (Glampedakis et al., 2002; Babak et al., 2007; Babak et al.,
2017) and galactic white dwarf binaries (Nelemans et al., 2001) all expected.
Inferring the properties of these systems from the observed gravitational-wave
signals will e.g. provide stringent tests of general relativity in the strong-
field regime, help constrain the Hubble’s constant, and investigate the effect
of dark matter on GW propagation (Amaro-Seoane et al., 2013). Unfortunately,
inferring the properties of LISA sources includes additional subtitles and
complexities beyond it’s ground-based counterpart. For instance, unlike
existing ground-based GW detectors that can be assumed static throughout the
observed GW signal, the LISA observatory is continuously orbiting the Sun.
This affects e.g. waveform model generation as well as the commonly used
assumption that the noise is stationary and uncorrelated. Accurately inferring
the properties of LISA sources is crucial to maximize the scientific output of
the LISA mission.

bilby , the user-friendly Bayesian inference library, is designed to infer the
source properties from the observed GW signal (Ashton et al., 2019; Romero-
Shaw et al., 2020). bilby has been extensively used to analyse data from
ground-based GW detectors, including individual GW observations (e.g. Abbott
et al., 2021a), population properties using hierarchical Bayesian inference
(e.g. Abbott et al., 2023; Saleem et al., 2022; Talbot et al., 2019), tests of
general relativity (e.g. Abbott et al., 2021b) and more (e.g. You et al.,
2021; Powell & Müller, 2022; Hübner et al., 2020; Dietrich et al., 2020; Hoy,
2022). Given that bilby is regularly used with a proven track record for
producing reliable results, it is natural to extend bilby to perform Bayesian
inference on LISA sources. It also has the advantage of trivially interfacing
with parallel-bilby  (Smith et al., 2020), which reduces wall-time through
massively parallel Bayesian inference, as well as numerous methods for
reducing overall cost (e.g. Cornish, 2010; Zackay et al., 2018; Cornish, 2021;
Morisaki et al., 2023; Morisaki, 2021; Williams et al., 2021; Krishna et al.,
2023).

In this paper, we outline our developments to bilby to perform Bayesian
inference on (simulated) LISA data. We demonstrate that full Bayesian
inference can be performed on the three time-delay-interferometry (TDI) LISA
observables (Tinto et al., 2002) for a range of waveform models. By focusing
on MBHB mergers, we show that bilby accurately infers the source properties in
zero-noise, as well as idealized instrumental noise. We compare the standard
likelihood commonly used in GW astronomy against optimised likelihoods
designed to reduce the computational cost (Cornish, 2010; Zackay et al., 2018;
Cornish, 2021; Krishna et al., 2023). We show that although the latter reduces
the CPU time by a factor of $\sim 10^{2}$, a small bias is observed. We
finally conclude by showing that bilby can perform Bayesian inference on LISA
data with higher order multipole waveform models, and discuss how the
inclusion of higher order multipoles can help break degeneracies in the
parameter space. Although the work presented here focuses on MBHBs, our
implementation is generic and can in principle perform Bayesian inference for
any transient GW signals in LISA data; our only requirement is GW models for
more exotic systems.

We note that there already exists methods to perform Bayesian inference on
LISA sources: balrog, a software package under development for LISA data
analysis (see e.g. Klein et al., 2022; Pratten et al., 2023), pycbc-inference
(Weaving et al., 2023), numerous algorithms that employ simplified fisher-
matrix approximations (see e.g. Vecchio, 2004; Berti et al., 2005; Lang &
Hughes, 2006; Arun et al., 2009) and more (see e.g. Cornish, 2022; Katz,
2022a; Littenberg & Cornish, 2023). However, our implementation is generic,
includes additional waveform models defined in LALSimulation (LIGO Scientific
Collaboration, 2018), has the ability to interface with a diverse range of
samplers, and builds upon the modularity and ease of using bilby , meaning
that it facilitates use by beginners.

The paper is organised as follows: in Section 2 we review Bayesian inference,
describe the likelihoods typically used in gravitational-wave astronomy, and
introduce the bilby package. In Section 2.1 we outline our modifications to
bilby to enable Bayesian inference on LISA data. In Section 3 we validate our
changes and confirm the effectiveness of bilby for analysing LISA data. We
infer the properties of MBHBs in zero-noise in Section 3.1, and idealised
instrumental noise in Section 3.2. Results of the analyses presented in this
paper are provided as accompaniments to this paper.

##  2 Bayesian inference

Bayesian inference is the process of estimating the parameters,
$\boldsymbol{\theta}$, of a model, $m$, given the observed data $d$. These
properties are represented by the model-dependent posterior probability
density function (PDF), which is calculated through Bayes’ theorem,

$p(\boldsymbol{\theta}|d,m)=\frac{p(\boldsymbol{\theta}|m)\,p(d|\boldsymbol{\theta},m)}{\mathcal{Z}}.$
(1)

Here, $p(\boldsymbol{\theta}|m)$ is the probability that the model has
parameters $\boldsymbol{\theta}$ given the model $m$, otherwise known as
prior, $p(d|\boldsymbol{\theta},m)$ is the probability of observing the data
$d$ given the model parameters $\boldsymbol{\theta}$ and model $m$, otherwise
known as the likelihood, and $\mathcal{Z}$ is the probability of observing the
data given the model $m$, otherwise known as the evidence.

Although an analytic form of the prior and likelihood is known, it is often
not possible to trivially evaluate Eq. 1. Consequently, stochastic sampling
techniques, such as Nested sampling (Skilling, 2004, 2006) or Markov-Chain-
Monte-Carlo (MCMC) (Metropolis & Ulam, 1949), are often used to draw samples
from the unknown posterior PDF (although see e.g. Pankow et al., 2015; Lange
et al., 2018; Delaunoy et al., 2020; Green et al., 2020; Chua & Vallisneri,
2020; Green & Gair, 2021; Dax et al., 2021; Gabbard et al., 2022; Tiwari et
al., 2023; Fairhurst et al., 2023, for alternative methods).

In GW astronomy, the noise in the GW detector is assumed stationary and
Gaussian, and the Whittle likelihood (Whittle, 1951) is often used (Romano &
Cornish, 2017; Thrane & Talbot, 2019). Here a parameterised GW model, $m$, is
evaluated for a set of parameters, $\boldsymbol{\theta}$, and compared with
the GW strain data, $d$. For a network of GW detectors, a Whittle likelihood
is assumed for each case, and the full likelihood is obtained by taking the
product of the single-detector likelihoods (see e.g. Veitch et al., 2015;
Ashton et al., 2019, for details). This inherently assumes that the noise is
uncorrelated in each detector. Of course, the GW observed in different
detectors will have small differences owing to their relative position on
Earth or in space. An additional step is therefore required to calculate the
expected GW as seen in each of the detectors in the network.

As the Whittle likelihood is typically computationally expensive to evaluate,
the heterodyned likelihood was introduced (Cornish, 2010; Zackay et al., 2018;
Cornish, 2021). By assuming that the model evaluated at two high likelihood
points is similar, and hence the ratio is a smoothly varying function, summary
data for a fiducial point can be pre-computed, and likelihoods for surrounding
points can be evaluated in a fraction of the time. It has been shown that the
heterodyned likelihood can reduce the computational cost of GW Bayesian
inference by a factor of $10^{4}$ compared to evaluating the Whittle
likelihood (Zackay et al., 2018). A downside of the heterodyned likelihood is
that a fiducial point of high likelihood must be known _a priori_.

###  2.1 bilby and modifications for LISA

Numerous packages have been developed to stochastically sample the unknown
posterior PDF for GW signals (e.g. Veitch et al., 2015; Ashton et al., 2019;
Biwer et al., 2019). Most notably, LALInference is a Bayesian inference
library that employs custom built Nested and MCMC algorithms, and was the main
workhorse of the LIGO–Virgo–KAGRA (LVK) collaboration for many years (Veitch
et al., 2015). bilby was later introduced as an alternative tool, which
prioritised modularity and ease of accessibility (Ashton et al., 2019). Unlike
LALInference , bilby interfaces with existing ‘off-the-shelf’ samplers (e.g.
Foreman-Mackey et al., 2013; Handley et al., 2015; Speagle, 2020; Williams et
al., 2021; Ashton & Talbot, 2021). For the case of the nested sampler dynesty
(Speagle, 2020), bilby incorporates additional modifications to enhance it’s
performance for GW problems. Since the third LVK GW observing run, bilby has
primarily been used to analyse data from ground-based GW detectors (Abbott et
al., 2021b, 2023).

Although bilby has a proven track record of analysing data from ground-based
GW detectors, there are several additional features that are required in order
to infer the properties of LISA sources. These include implementing the LISA
observatory into the bilby.gw.detector module, described in Section 2.1.1, and
incorporating additional models for calculating GWs as seen by LISA, see
Section 2.1.2.

####  2.1.1 The LISA observatory

LISA will consist of three spacecraft, each containing a free falling test
mass. Lasers will be relayed between each of the spacecraft, and the distance
between free falling test masses will be monitored over time. A major
challenge is that the laser frequency noise will be several orders of
magnitude larger than any change in the arm length due to expected GW signals
(Amaro-Seoane et al., 2017). For this reason, Time-delayed interferometery
(TDI) was introduced as an additional post-processing step. TDI is a technique
that combines data from three spacecraft into virtual interferometric
observables, where the laser frequency noise is suppressed by several orders
of magnitude, improving the sensitivity of LISA to GW signals. Numerous TDI
observables have been proposed with different underlying assumptions. The
second-generation TDI observables, $X$, $Y$ and $Z$, allowed for linearly
changing arm lengths, making them the most applicable for realistic LISA
scenarios (Tinto & Armstrong, 1999; Armstrong et al., 1999; Estabrook et al.,
2000; Vallisneri, 2005; Tinto & Dhurandhar, 2021; Tinto et al., 2023).
However, the second-generation TDI observables are correlated in their noise
properties. For this reason, a linear transformation is applied to the $X$,
$Y$ and $Z$ observables, forming a set of approximately uncorrelated TDI
observables: $A$, $E$ and $T$ (Prince et al., 2002).

We provide the LISA class to handle the LISA observatory. We treat LISA as a
network of distinct virtual interferometers, each corresponding to separate
TDI observables. Although our implementation is generic to any TDI
observables, allowing for possible improvements beyond second-generation TDI,
the likelihood for a network of virtual interferometers assumes the noise is
uncorrelated in each detector. This paper therefore uses the $A$, $E$ and $T$
channels.

####  2.1.2 Alternative GW models

In order to evaluate the GW likelihood, a parameterised model must be
evaluated for a given set of parameters, and compared with the GW strain data.
Since it is sensible to perform Bayesian inference on the TDI observables, see
Section 2.1.1, our parameterised GW model must calculate the expected GW in
each of these channels.

An additional complication is that the GWs from expected LISA sources will
likely remain observable for weeks, if not months. During this time, the LISA
observatory will continuously change it’s orientation with respect to the
source, as a result of it’s orbit around the Sun. This implies that when
calculating the expected GW as seen by LISA, we must apply a time and
frequency dependent transformation, $\mathcal{T}(f,t)$, to the source-frame GW
(Marsat et al., 2021). This is in contrast to ground-based GW detectors, which
are assumed to be static during the whole GW signal.

The lisa_binary_black_hole function returns a dictionary of expected GWs
observed in each TDI observable for a provided set of parameters
$\boldsymbol{\theta}$. Since we argued that the $A$, $E$ and $T$ channels are
the most logical for calculating the likelihood in Section 2.1.1, we currently
only consider these TDI observables. Of course, this can be extended in future
work. Our implementation calculates GWs for parameters provided in the solar
system baricenter (SSB) or LISA frame. When the LISA frame is used,
transformations are applied to map the parameters into the SSB frame. This is
required in order to calculate $\mathcal{T}(f,t)$, which depends on the source
location and polarization defined in the SSB frame (Katz, 2022a).

We interface with the BBHx (Katz et al., 2020; Katz, 2022b, 2021) and
LALSimulation (LIGO Scientific Collaboration, 2018) libraries. BBHx is an open
source software that generates the GW polarizations, calculates the time and
frequency dependent transformations, $\mathcal{T}(f,t)$, and returns the
expected GWs in the $A$, $E$ and $T$ channels. Currently BBHx implements the
IMRPhenomD (Khan et al., 2016; Husa et al., 2016) and IMRPhenomHM (London et
al., 2018) waveform models; both models assume black hole spins aligned with
the orbital angular momentum, with the former considering only the leading
order quadrupole contribution, while the latter allows for additional higher-
order multipole content. The user can call the BBHx models through bilby by
passing a waveform_approximant of BBHx_IMRPhenomD or BBHx_IMRPhenomHM
respectively. In order to reduce computational cost, interpolation can be
enabled in the BBHx waveform models by specifying the waveform argument
direct=False. Here, a reduced set of frequencies are used to evaluate the
waveform, controlled with the waveform argument length, and interpolation is
used to provide the waveform at the desired frequency points (Katz et al.,
2020; Katz, 2022b, 2021); by default length = 1024. In this work, we do not
use interpolation when using the BBHx waveform models, in order to reduce the
possibility of interpolation errors. If the user would like to use additional
models beyond those included in BBHx, we use the LALSimulation library to
generate GW polarizations, and use BBHx functions to project the GW
polarizations into the $A$, $E$ and $T$ channels. This allows us to interface
with e.g. the cutting edge binary black hole model IMRPhenomXPHM (Pratten et
al., 2021).

Figure 1: Mismatches between $10^{5}$ GW signals calculated using the BBHx and
LALSimulation implementations of IMRPhenomD, when projected onto the A TDI
channel.

To validate our LALSimulation implementation, we randomly drew $10^{5}$ MBHBs
with chirp masses $1<\mathcal{M}\,[10^{6}\,M_{\odot}]<5.0$, mass ratios
$0.25<q<1.0$, spins $0<\chi<0.99$ and luminosity distances
$1<d_{L}[\text{Gpc}]<100$ (all other parameters randomly chosen from an
unconstrained boundary), and calculated the expected GW in the $A$ TDI channel
when using the BBHx IMRPhenomD implementation, and the LALSimulation
IMRPhenomD implementation. We quantify the difference between GWs by the
mismatch, which is defined as,

$\mathfrak{M}=1-\max_{\phi,t}\left[4\,\mathrm{Re}\int_{0}^{\infty}df\frac{\tilde{a}(f)\tilde{b}^{\ast}(f)}{S(f)}\right]$
(2)

where $a$ and $b$ are the GWs calculated using BBHx and LALSimulation
respectively, $S$ is the one-sided power spectral density (PSD), and the
maximisation is over the phase $\phi$ and time $t$111By maximising over time,
we are implicitly assuming that the LISA response function is not varying in
time across the duration of this signal. However, because the majority of the
signal-to-noise ratio from MBHBs occurs in the final day before merger (Marsat
et al., 2021), it is safe to assume that LISA is approximately stationary
during this time period.. The mismatch ranges between $0$ and $1$, where a
mismatch of $0$ implies the two GWs are identical. As shown in Figure 1, we
find a median mismatch of $4.0\times 10^{-6}$, with a large tail extending to
smaller mismatches. Interestingly, we find that in general larger mismatches
are correlated with smaller chirp masses, where the signal durations are
longer. Ideally all mismatches would be close to numerical precision, but we
point out that the average mismatch remains two orders of magnitude lower than
the mismatch between IMRPhenomD and numerical solutions to general relativity
(Khan et al., 2016). It is well known that the mismatch can be mapped to an
approximate signal-to-noise ratio (SNR) at which two signals are
indistinguishable at 90% confidence (Baird et al., 2013). Assuming eleven
degrees of freedom (for example the two masses and spins of both black holes,
sky location, luminosity distance, inclination, phase, polarization and merger
time), we find that on average, the two approaches are indistinguishable for
SNRs less than $\sim 1400$, and $70\%$ of signals are indistinguishable for
SNRs less than $\sim 1000$. Therefore although there may be convention
differences between the two approaches, which become more evident at larger
durations, it is unlikely to cause significant differences when performing
Bayesian inference on most MBHB systems. We are therefore confident that our
LALSimulation implementation is reasonable.

##  3 Results

To test our developments to bilby to perform Bayesian inference on LISA data,
we analyse MBHB signals. We choose signals that are similar to those expected
when LISA is online. For cases where the true values are known, systematic
biases can be studied.

We perform Bayesian inference with the dynesty Nested sampler (Speagle, 2020),
and compare the computational cost and performance of evaluating the full
Whittle likelihood with the heterodyned likelihood. Since evaluating the full
Whittle likelihood is likely computationally expensive, we perform Bayesian
inference with parallel-bilby . We use standard bilby when evaluating the
heterodyned likelihood, and set the true values as the fiducial point. Of
course, in reality the true values are unknown, meaning that this is an
idealised case. For all analyses we use 1000 live points, the bilby-
implemented rwalk sampling algorithm with an average of 60 accepted steps per
MCMC, and evaluate the likelihood for frequencies $10^{-4}\,\mathrm{Hz}\leq
f\leq 10^{-1}\,\mathrm{Hz}$. Sampling is performed in the LISA frame, as it is
the most natural for LISA data analysis tasks, and we coherently analyse the
$A$, $E$ and $T$ channels. All analyses use uninformative and wide priors, and
we sample in the standard parameters: chirp mass, mass ratio, spin magnitudes,
inclination angle, phase, ecliptic latitude, ecliptic longitude, polarization
and merger time. We reduce the computational cost by analytically
marginalizing over the luminosity distance of the source, and reconstructing
the posterior in post-processing, see Thrane & Talbot (2019) for details.

To compare posterior distributions obtained with the heterodyned and full
Whittle likelihoods, we use the Jensen-Shannon divergence (JSD) (Lin, 1991).
The JSD ranges between $0\,\mathrm{bits}$ and $1\,\mathrm{bit}$, where a
JSD=$0\,\mathrm{bits}$ (JSD=$1\,\mathrm{bit}$) implies statistically identical
(distinct) distributions. The JSD has been regularly used in GW astronomy for
comparing posterior distributions (see e.g. Abbott et al., 2019, 2021c). A
general rule of thumb is that a JSD $<0.05\,\mathrm{bits}$ implies that the
distributions are in good agreement (Abbott et al., 2019).

We analyse a zero-noise injection in Section 3.1, and an injection in mock
LISA data with idealised Gaussian instrumental noise in Section 3.2. We
obtained mock LISA data through the Sangria dataset (Le Jeune & Babak, 2022)
released as part of the LISA data challenge (LDC). The LDC contains multiple
challenges that will eventually cover all of the complexities associated with
LISA data analysis. For this paper, we only consider the Sangria dataset as a
proof of principle, and leave further challenges to future work. The PSDs used
in this work were generated from the Sangria dataset using gwpy (Macleod et
al., 2021) and the Welch method. Unlike Weaving et al. (2023), we assume
Galactic binaries form part of the noise, and therefore we include them in our
PSD generation. By constructing PSDs in this way we are inherently removing
the time-dependence, both from the foreground Galactic binaries and from
LISA’s orbit around the Sun. This is therefore an idealised case, but we note
that it will unlikely cause a significant difference for MBHBs, where the SNRs
are expected to be large, and rapidly accumulating within the final day before
merger (Marsat et al., 2021). This assumption will need to be relaxed for
signals of longer duration. We note that the foreground Galactic binary
problem may be navigated by using a Global Fit approach (see e.g. Littenberg &
Cornish, 2023), and the PSD may be able to include a time-frequency dependence
by using the methods outlined in e.g. Cornish (2020).

Figure 2: Inferred posterior distribution for the total mass $M$, mass ratio
$q$, primary spin $\chi_{1}$ and secondary spin $\chi_{2}$ when analysing a
zero noise injection with the BBHx_IMRPhenomD model. Results were obtained
when analysing the A, E and T TDI channels. The dotted and solid lines show
the results when using the heterodyned and full Whittle likelihood
respectively. The grey cross hairs show the true values, and the 2-dimensional
contours show the 90% confidence interval.

###  3.1 Zero-noise Injections

We first simulate a GW produced by a MBHB with total mass
$M=m_{1}+m_{2}=3.1\times 10^{6}\,M_{\odot}$, mass ratio $q=m_{2}/m_{1}=0.74$,
primary spin $\chi_{1}=0.51$ and secondary spin $\chi_{2}=0.14$. We assume the
source has spins aligned with the orbital angular momentum and is located at
an ecliptic latitude $\beta=2.8$ and ecliptic longitude $\lambda=2.3$, viewed
at an inclination angle of $\iota=1.3\,\mathrm{rad}$ and luminosity distance
$d_{L}=17\,\mathrm{Gpc}$. The phase and polarization were chosen to be
$\phi=4.0$, $\psi=1.1$ respectively, and the merger time was set to be
$t_{c}=6.7\times 10^{6}\,\mathrm{s}$. The injection is defined in the LISA
frame, and has SNR $\sim 550$222The parameters were chosen to be similar to
the zero-noise injection performed in Weaving et al. (2023).. We generate the
simulated signal with BBHx_IMRPhenomD, inject it into a month’s worth of zero-
noise data sampled at 5 second intervals, and perform Bayesian inference with
the BBHx_IMRPhenomD model for consistency.

Figure 3: _Top_ : Comparison between the reconstructed GW signal (black) and
the time-domain Sangria blind dataset (grey). The reconstructed GW signal was
obtained by coherently analysing the A, E and T TDI channels using the
aligned-spin limit of IMRPhenomXPHM, and the full Whittle likelihood. We shift
the reconstructed GW signal and data by the inferred merger time as defined in
the LISA frame $t_{c}$. The reconstructed GW signal is plotted as a band
representing the 90% credible region. _Bottom_ : Comparison between simulated
gaussian noise which has been coloured by the power spectral density (black),
and the residual after subtracting the reconstructed GW signal from the time-
domain Sangria blind dataset (grey). For simplicity, we only show the A TDI
channel in both panels.

In Figure 2, we compare the marginalized posterior distributions for the
intrinsic parameters (total mass, mass ratio, primary spin and secondary spin)
with the true values. In general we see that the full Whittle and heterodyned
likelihoods recover the true values within the 90% confidence interval.
Although all Jensen-Shannon divergences are less 0.05 bits, we see that the
Whittle likelihood results more accurately recover the true values, with a
small shift relative to the heterodyned likelihood in all cases.

When considering extrinsic parameters (merger time, sky location, luminosity
distance, inclination, phase, polarization), we expect to observe a bias in
the sky location owing to well known degeneracies; for aligned-spin leading
order quadrupole models, we expect to infer one of eight possibilities (Marsat
et al., 2021). Both likelihoods find support at the true ecliptic latitude,
but only the heterodyned likelihood identifies the true ecliptic longitude
mode, while the full Whittle likelihood favours the opposite sky mode. We see
that both likelihoods recover a quad-modal polarization, with the full Whittle
posterior shifted by $\pi/2$ compared to the heterodyned posterior – this is
consistent with favouring the opposite sky mode. Both likelihoods accurately
recover the merger time and phase. Interestingly, only the Whittle likelihood
observes a bias in the inclination angle and luminosity distance. Since we are
using a model with spins aligned with the orbital angular momentum, and
ignoring the effects of higher-order multipole content, we are unable to break
the well known distance–inclination degeneracy (see e.g. Marsat et al., 2021).
It is therefore expected to see a bias in both the distance and inclination
angle, especially at such large SNRs. Given that the heterodyned likelihood
uses the true values as the fiducial point, it is not surprising that it more
accurately recovers the distance and inclination angle of the source.

We found that, although both analyses obtained $\sim 5,500$ posterior samples
and evaluated the likelihood $\sim 200$ million times, the analysis that took
advantage of the heterodyned likelihood was $\sim 10^{2}$ times faster and
required $\sim 10$ times less memory than the full Whittle likelihood
analysis; the full Whittle likelihood analysis completed in $O(10^{4})$ CPU
hours and required $\sim 200\,\mathrm{GB}$ memory. We note, however, that the
heterodyned likelihood requires prior knowledge of a high likelihood point,
whilst the Whittle likelihood is completely general. To navigate this,
previous works either used the true injected values, or they performed
template-based matched-filter searches of the data prior to performing
Bayesian inference, and they used the best fit template as the high likelihood
point. We note that in some cases, the latter can lead to relatively poor
fitting templates (Weaving et al., 2023). Standard implementations of the
heterodyned likelihood may also not be appropriate for models that include the
effects of higher-order multipoles (Leslie et al., 2021). The full Whittle
likelihood may therefore be more appropriate for analyses where a high
likelihood fiducial point cannot be confidently found, or when models that
include higher order effects are required. We leave an investigation that
studies the impact of fiducial points, as well as potential algorithms for
using the heterodyned posteriors as an initial guess for massively restricting
the prior volume for full Whittle analyses, to future work.

###  3.2 Sangria Blind Dataset

We next analyse the Sangria Blind Dataset (Le Jeune & Babak, 2022) with
IMRPhenomXPHM (Pratten et al., 2021). IMRPhenomXPHM is the cutting edge binary
black hole model, which includes spin-induced orbital precession (Apostolatos
et al., 1994) and higher order multipoles. In this analysis, we restrict black
holes spins to be aligned with the orbital angular momentum, meaning that we
ignore the effects of spin-induced orbital precession. Of course, we could
have alternatively used IMRPhenomXHM (García-Quirós et al., 2020), which also
includes higher order multipole contributions and similarly neglects the
effects of spin-induced orbital precession. However, we chose to use the more
general IMRPhenomXPHM model, as it is more likely to be used in future
analyses; for instance investigating the impact of including mis-aligned black
hole spins for breaking possible degeneracies. This dataset contains idealized
Gaussian instrumental noise and simulated waveforms from Galactic white dwarf
binaries and MBHBs with parameters derived from astrophysical models. The
level of instrumental noise, the number of sources, and their parameters are
not disclosed. Since the heterodyned likelihood requires a fiducial point, and
the true injected parameters are unknown, we only use the full Whittle
likelihood. We analyse a year’s worth of data, sampled at 5 second intervals,
and generate PSDs using all available data. The PSDs differ from those used in
Section 3.1, but we similarly assume that Galactic binaries form part of the
noise, and therefore we include them in our PSD generation.

As with the current bilby implementations, we require an estimate for the
merger time of the binary in order to define a suitable prior distribution for
the merger time. This is typically provided from the matched-filter search
pipelines. However, given that signal 0 can clearly be seen in the Sangria
blind dataset, see Figure 3, it is trivial to estimate a suitable prior
distribution that encompasses the merger. In future analyses of blind data, we
would base our prior distributions on the best matching template from a LISA
search pipeline, see e.g. Weaving et al. (2023).

We see remarkable agreement between the reconstructed GW obtained with
IMRPhenomXPHM and the data, see Figure 3. When comparing the marginalized
posterior distributions with those obtained by Weaving et al. (2023) and
BBHx_IMRPhenomD (see Figs. 9 and 10 in Weaving et al. (2023)), we see that our
analysis a) more tightly constrains the intrinsic source properties, and b)
obtains smoother (more gaussian-like) posterior distributions. This is
expected as significant SNRs, as well as a lack of observed SNRs, in higher
order multipoles can help break degeneracies (see e.g. Marsat et al., 2021;
Pratten et al., 2023; Gong et al., 2023). We can directly calculate the SNRs
in higher order multipoles given the inferred posterior with simple-pe
(Fairhurst et al., 2023). For this signal, we infer SNRs in the
$(\ell,m)=(3,3)$ and $(4,4)$ multipoles to be $<1$ ($\sim 0.6$ and $\sim 0.3$
respectively). The near-zero SNRs are consistent with the expectations from
gaussian noise (Mills & Fairhurst, 2021), and therefore we find no significant
evidence for higher order multipole content in this signal. This is consistent
with the fact that the GW signals in the Sangria Blind Dataset were
constructed with IMRPhenomD (Le Jeune & Babak, 2022).

When subtracting the reconstructed GW from the data, the residual signal is
comparable with simulated gaussian noise coloured by the PSD, see Figure 3.
This implies that the MBHB signal has been appropriately subtracted from the
data, meaning that the residual can be used for further Bayesian inference
analyses as part of a Global Fit style approach (see e.g. Littenberg &
Cornish, 2023). We note that do not expect perfect agreement between the
residual and the simulated gaussian noise, since the simulated gaussian noise
is dependent on the specific noise realisation chosen.

The IMRPhenomXPHM analysis completed in $O(10^{5})$ CPU hours. Owing to the
larger data duration considered, compared with the zero-noise injection
presented in Section 3.1, more memory was required, $\sim 300\,\mathrm{GB}$.
For comparison a typical analysis of 8 seconds worth of ground-based GW
detector data sampled at $1/4096$ second intervals completes in $O(10^{3})$
CPU hours (Colleoni et al., 2021; Ramos-Buades et al., 2023). The increase in
computational cost is primarily driven by the waveform evaluation time and the
long data durations considered. There are numerous methods available for
reducing the computational cost without taking advantage of the heterodyned
likelihood. For example, reduced order quadrature models (Canizares et al.,
2015; Qi & Raymond, 2021; Morisaki et al., 2023), adaptive frequency
resolutions (Vinciguerra et al., 2017; Morisaki, 2021) and other techniques
(Williams et al., 2021; Lee et al., 2022; Islam et al., 2022; Wong et al.,
2023; Pathak et al., 2022; Tiwari et al., 2023) have all been used to good
effect in the past. Given that most of these methods are already implemented
in the bilby infrastructure, it is trivial to take advantage of these in
future analyses. As discussed throughout this work, the SNR from MBHBs rapidly
accumulates within the final day before merger (Marsat et al., 2021). This
implies that long data duration’s may not be necessary into order accurately
infer the properties of MBHBs. This will inherently reduce both the CPU and
memory costs.

Although the scalability of parallel-bilby has been investigated in detail in
Smith et al. (2020), we provide a further analysis specific to the LISA
context. We analysed the Sangria Blind Dataset with the aligned-spin limit of
IMRPhenomXPHM using 128, 256, 512, and 768 CPUs. Owing to computational cost,
we stopped the analyses after 12 hours wall time and compared the samplers
progress. In general, we find that more CPUs increased a) the number of times
the Whittle likelihood was called, b) the recorded iteration, and c) the
number of proposals used. Specifically, when comparing the 768 and 128 CPU
analyses, we find that the former called the likelihood $\sim 34$ times more
than the latter, and it was on the $\sim 5000$th iteration compared to $\sim
3500$th. This suggests that running with more CPUs will reduce the overall
wall time of the analysis.

##  4 Discussion

Inferring the properties of LISA sources is crucial for maximising the
scientific output of the LISA mission. Although numerous algorithms exist for
performing Bayesian inference, few are applicable for GW astronomy, and even
fewer are capable of analysing LISA data. In this work, we introduce a
modified version of bilby , the user-friendly Bayesian inference library, to
allow users to easily analyse LISA data with a suite of waveform models. We
use bilby as our starting point as it is already extensively used for
analysing data from ground-based GW detectors, and has a proven track record
for producing reliable results.

By focusing on MBHBs, we show that bilby can perform Bayesian inference on
LISA data and is capable of accurately inferring the source properties in
zero-noise as well as idealized instrumental noise. We compare the performance
of the standard likelihood used in GW astronomy, with an optimised likelihood
designed to reduce the computational cost: the heterodyned likelihood.
Although we find that analyses which use the full likelihood more accurately
infer the source properties, analyses using the heterodyned likelihood are
$\sim 10^{2}$ times faster. We note, however, that the heterodyned likelihood
is dependent on knowing a point of high likelihood _a priori_ , and standard
implementations may not be suitable for models with higher order multipole
content. The full Whittle likelihood may therefore be preferred for many
cases. We finally demonstrate that higher order multipole waveform models can
be performed on LISA data. We highlight that the computational cost can likely
be reduced by using techniques other than the heterodyned likelihood, already
implemented into the bilby infrastructure.

Our current implementation includes waveform models that are applicable for
MBHBs. As we have maintained the modularity of bilby , we anticipate future
extensions of this work to include alternative waveform models (Afshordi et
al., 2023).

## Acknowledgements

We would like to thank Jonathan Thompson for valuable comments on this
manuscript, as well as Gregory Ashton, Ian Harry and Connor Weaving for
continued discussions throughout this project. We are also grateful to Gregory
Ashton and Colm Talbot for their advice regarding merging our modifications
into the main bilby code-base. We thank the UKRI Future Leaders Fellowship for
support through the grant MR/T01881X/1. This work used the computational
resources provided by the ICG, SEPNet and the University of Portsmouth,
supported by STFC grant ST/N000064, and the DiRAC@Durham facility managed by
the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC
Facility (www.dirac.ac.uk). The equipment for the latter was funded by BEIS
capital funding via STFC capital grants ST/P002293/1, ST/R002371/1 and
ST/S002502/1, Durham University and STFC operations grant ST/R000832/1. DiRAC
is part of the National e-Infrastructure. This work used bilby $=$ v2.1.0,
bilby_pipe $=$ v1.0.11, parallel-bilby  $=$ v2.0.2, and dynesty $=$ v2.1.1.
PESummary (Hoy & Raymond, 2021) and matplotlib (Hunter, 2007) were used for
plotting.

_Data availability statement_ : In this work we analysed simulated and
publicly available data (Le Jeune & Babak, 2022). The results presented in
this paper are publicly available online (Hoy & Nuttall, 2023).

## References

  * Aasi et al. (2015) Aasi J., et al., 2015, Class. Quant. Grav., 32, 074001
  * Abbott et al. (2019) Abbott B. P., et al., 2019, Phys. Rev. X, 9, 031040
  * Abbott et al. (2021a) Abbott R., et al., 2021a, arXiv e-prints, p. arXiv:2111.03606
  * Abbott et al. (2021b) Abbott R., et al., 2021b, arXiv e-prints, p. arXiv:2112.06861
  * Abbott et al. (2021c) Abbott R., et al., 2021c, Phys. Rev. X, 11, 021053
  * Abbott et al. (2023) Abbott R., et al., 2023, Phys. Rev. X, 13, 011048
  * Acernese et al. (2014) Acernese F., et al., 2014, Classical and Quantum Gravity, 32, 024001
  * Afshordi et al. (2023) Afshordi N., et al., 2023, arXiv preprint arXiv:2311.01300
  * Akutsu et al. (2021) Akutsu T., et al., 2021, PTEP, 2021, 05A101
  * Amaro-Seoane et al. (2013) Amaro-Seoane P., et al., 2013, GW Notes, 6, 4
  * Amaro-Seoane et al. (2017) Amaro-Seoane P., et al., 2017, arXiv e-prints, p. arXiv:1702.00786
  * Apostolatos et al. (1994) Apostolatos T. A., Cutler C., Sussman G. J., Thorne K. S., 1994, Phys. Rev. D, 49, 6274
  * Armstrong et al. (1999) Armstrong J., Estabrook F., Tinto M., 1999, The Astrophysical Journal, 527, 814
  * Arun et al. (2009) Arun K. G., et al., 2009, Class. Quant. Grav., 26, 094027
  * Ashton & Talbot (2021) Ashton G., Talbot C., 2021, Mon. Not. Roy. Astron. Soc., 507, 2037
  * Ashton et al. (2019) Ashton G., et al., 2019, Astrophys. J. Suppl., 241, 27
  * Babak et al. (2007) Babak S., Fang H., Gair J. R., Glampedakis K., Hughes S. A., 2007, Phys. Rev. D, 75, 024005
  * Babak et al. (2017) Babak S., et al., 2017, Phys. Rev. D, 95, 103012
  * Baird et al. (2013) Baird E., Fairhurst S., Hannam M., Murphy P., 2013, Phys. Rev. D, 87, 024035
  * Berti et al. (2005) Berti E., Buonanno A., Will C. M., 2005, Phys. Rev. D, 71, 084025
  * 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
  * Canizares et al. (2015) Canizares P., Field S. E., Gair J., Raymond V., Smith R., Tiglio M., 2015, Phys. Rev. Lett., 114, 071104
  * Chua & Vallisneri (2020) Chua A. J. K., Vallisneri M., 2020, Phys. Rev. Lett., 124, 041102
  * Colleoni et al. (2021) Colleoni M., Mateu-Lucena M., Estellés H., García-Quirós C., Keitel D., Pratten G., Ramos-Buades A., Husa S., 2021, Phys. Rev. D, 103, 024029
  * Cornish (2010) Cornish N. J., 2010, arXiv e-prints, p. arXiv:1007.4820
  * Cornish (2020) Cornish N. J., 2020, arXiv e-prints, p. arXiv:2009.00043
  * Cornish (2021) Cornish N. J., 2021, Phys. Rev. D, 104, 104054
  * Cornish (2022) Cornish N. J., 2022, Phys. Rev. D, 105, 044007
  * Danzmann (1996) Danzmann K., 1996, Class. Quant. Grav., 13, A247
  * Dax et al. (2021) Dax M., Green S. R., Gair J., Macke J. H., Buonanno A., Schölkopf B., 2021, Phys. Rev. Lett., 127, 241103
  * Delaunoy et al. (2020) Delaunoy A., Wehenkel A., Hinderer T., Nissanke S., Weniger C., Williamson A. R., Louppe G., 2020, arXiv e-prints, p. arXiv:2010.12931
  * Dietrich et al. (2020) Dietrich T., Coughlin M. W., Pang P. T. H., Bulla M., Heinzel J., Issa L., Tews I., Antier S., 2020, Science, 370, 1450
  * Estabrook et al. (2000) Estabrook F. B., Tinto M., Armstrong J. W., 2000, Phys. Rev. D, 62, 042002
  * Fairhurst et al. (2023) Fairhurst S., Hoy C., Green R., Mills C., Usman S. A., 2023, Phys. Rev. D, 108, 082006
  * Foreman-Mackey et al. (2013) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, Publ. Astron. Soc. Pac., 125, 306
  * Gabbard et al. (2022) Gabbard H., Messenger C., Heng I. S., Tonolini F., Murray-Smith R., 2022, Nature Phys., 18, 112
  * García-Quirós et al. (2020) García-Quirós C., Colleoni M., Husa S., Estellés H., Pratten G., Ramos-Buades A., Mateu-Lucena M., Jaume R., 2020, Phys. Rev. D, 102, 064002
  * Glampedakis et al. (2002) Glampedakis K., Hughes S. A., Kennefick D., 2002, Phys. Rev. D, 66, 064005
  * Gong et al. (2023) Gong Y., Cao Z., Zhao J., Shao L., 2023, Phys. Rev. D, 108, 064046
  * Green & Gair (2021) Green S. R., Gair J., 2021, Mach. Learn. Sci. Tech., 2, 03LT01
  * Green et al. (2020) Green S. R., Simpson C., Gair J., 2020, Phys. Rev. D, 102, 104057
  * Handley et al. (2015) Handley W. J., Hobson M. P., Lasenby A. N., 2015, Mon. Not. Roy. Astron. Soc., 453, 4385
  * Hoy (2022) Hoy C., 2022, Phys. Rev. D, 106, 083003
  * Hoy & Nuttall (2023) Hoy C., Nuttall L., 2023, _in preparation_
  * Hoy & Raymond (2021) Hoy C., Raymond V., 2021, SoftwareX, 15, 100765
  * Hübner et al. (2020) Hübner M., Talbot C., Lasky P. D., Thrane E., 2020, Phys. Rev. D, 101, 023011
  * Hunter (2007) Hunter J. D., 2007, CSE, 9, 90
  * Husa et al. (2016) Husa S., Khan S., Hannam M., Pürrer M., Ohme F., Forteza X. J., Bohé A., 2016, Phys. Rev. D, 93, 044006
  * Islam et al. (2022) Islam T., Roulet J., Venumadhav T., 2022, arXiv e-prints, p. arXiv:2210.16278
  * Katz (2021) Katz M., 2021, mikekatz04/BBHx: First official public release, doi:10.5281/zenodo.5730688, https://doi.org/10.5281/zenodo.5730688
  * Katz (2022a) Katz M. L., 2022a, Physical Review D, 105, 044055
  * Katz (2022b) Katz M. L., 2022b, Phys. Rev. D, 105, 044055
  * Katz et al. (2020) Katz M. L., Marsat S., Chua A. J. K., Babak S., Larson S. L., 2020, Phys. Rev. D, 102, 023033
  * Khan et al. (2016) Khan S., Husa S., Hannam M., Ohme F., Pürrer M., Forteza X. J., Bohé A., 2016, Phys. Rev. D, 93, 044007
  * Klein et al. (2016) Klein A., et al., 2016, Phys. Rev. D, 93, 024003
  * Klein et al. (2022) Klein A., et al., 2022, arXiv e-prints, p. arXiv:2204.03423
  * Krishna et al. (2023) Krishna K., Vijaykumar A., Ganguly A., 2023, in preparation
  * LIGO Scientific Collaboration (2018) LIGO Scientific Collaboration 2018, LIGO Algorithm Library - LALSuite, free software (GPL), doi:10.7935/GT1W-FZ16
  * Lang & Hughes (2006) Lang R. N., Hughes S. A., 2006, Phys. Rev. D, 74, 122001
  * Lange et al. (2018) Lange J., O’Shaughnessy R., Rizzo M., 2018, arXiv e-prints, p. arXiv:1805.10457
  * Le Jeune & Babak (2022) Le Jeune M., Babak S., 2022, LISA Data Challenge Sangria (LDC2a), doi:10.5281/zenodo.7132178, https://doi.org/10.5281/zenodo.7132178
  * Lee et al. (2022) Lee E., Morisaki S., Tagoshi H., 2022, Phys. Rev. D, 105, 124057
  * Leslie et al. (2021) Leslie N., Dai L., Pratten G., 2021, Phys. Rev. D, 104, 123030
  * Lin (1991) Lin J., 1991, IEEE Transactions on Information Theory, 37, 145
  * Littenberg & Cornish (2023) Littenberg T. B., Cornish N. J., 2023, Phys. Rev. D, 107, 063004
  * London et al. (2018) London L., et al., 2018, PhRvL., 120, 161102
  * Macleod et al. (2021) Macleod D. M., Areeda J. S., Coughlin S. B., Massinger T. J., Urban A. L., 2021, SoftwareX, 13, 100657
  * Marsat et al. (2021) Marsat S., Baker J. G., Dal Canton T., 2021, Phys. Rev. D, 103, 083011
  * Mehta et al. (2023) Mehta A. K., Olsen S., Wadekar D., Roulet J., Venumadhav T., Mushkin J., Zackay B., Zaldarriaga M., 2023, arXiv e-prints, p. arXiv:2311.06061
  * Metropolis & Ulam (1949) Metropolis N., Ulam S., 1949, Journal of the American statistical association, 44, 335
  * Mills & Fairhurst (2021) Mills C., Fairhurst S., 2021, Phys. Rev. D, 103, 024042
  * Morisaki (2021) Morisaki S., 2021, Phys. Rev. D, 104, 044062
  * Morisaki et al. (2023) Morisaki S., Smith R., Tsukada L., Sachdev S., Stevenson S., Talbot C., Zimmerman A., 2023, arXiv e-prints, p. arXiv:2307.13380
  * Nelemans et al. (2001) Nelemans G., Yungelson L. R., Portegies Zwart S. F., 2001, Astron. Astrophys., 375, 890
  * Nitz et al. (2023) Nitz A. H., Kumar S., Wang Y.-F., Kastha S., Wu S., Schäfer M., Dhurkunde R., Capano C. D., 2023, Astrophys. J., 946, 59
  * Pankow et al. (2015) Pankow C., Brady P., Ochsner E., O’Shaughnessy R., 2015, Phys. Rev. D, 92, 023002
  * Pathak et al. (2022) Pathak L., Reza A., Sengupta A. S., 2022, arXiv e-prints, p. arXiv:2210.02706
  * Powell & Müller (2022) Powell J., Müller B., 2022, Phys. Rev. D, 105, 063018
  * Pratten et al. (2021) Pratten G., et al., 2021, Phys. Rev. D, 103, 104056
  * Pratten et al. (2023) Pratten G., Klein A., Moore C. J., Middleton H., Steinle N., Schmidt P., Vecchio A., 2023, Phys. Rev. D, 107, 123026
  * Prince et al. (2002) Prince T. A., Tinto M., Larson S. L., Armstrong J. W., 2002, Phys. Rev. D, 66, 122002
  * Qi & Raymond (2021) Qi H., Raymond V., 2021, Phys. Rev. D, 104, 063031
  * Ramos-Buades et al. (2023) Ramos-Buades A., Buonanno A., Estellés H., Khalil M., Mihaylov D. P., Ossokine S., Pompili L., Shiferaw M., 2023, arXiv preprint arXiv:2303.18046
  * Romano & Cornish (2017) Romano J. D., Cornish N. J., 2017, Living Rev. Rel., 20, 2
  * Romero-Shaw et al. (2020) Romero-Shaw I. M., et al., 2020, Mon. Not. Roy. Astron. Soc., 499, 3295
  * Saleem et al. (2022) Saleem M., Krishnendu N. V., Ghosh A., Gupta A., Del Pozzo W., Ghosh A., Arun K. G., 2022, Phys. Rev. D, 105, 104066
  * Skilling (2004) Skilling J., 2004, in AIP Conference Proceedings. AIP, doi:10.1063/1.1835238, https://doi.org/10.1063/1.1835238
  * Skilling (2006) Skilling J., 2006, Bayesian Analysis, 1, 833
  * Smith et al. (2020) Smith R. J. E., Ashton G., Vajpeyi A., Talbot C., 2020, Mon. Not. Roy. Astron. Soc., 498, 4492
  * Speagle (2020) Speagle J. S., 2020, Mon. Not. Roy. Astron. Soc., 493, 3132
  * Talbot et al. (2019) Talbot C., Smith R., Thrane E., Poole G. B., 2019, Phys. Rev. D, 100, 043030
  * Thrane & Talbot (2019) Thrane E., Talbot C., 2019, Publ. Astron. Soc. Austral., 36, e010
  * Tinto & Armstrong (1999) Tinto M., Armstrong J. W., 1999, Phys. Rev. D, 59, 102003
  * Tinto & Dhurandhar (2021) Tinto M., Dhurandhar S. V., 2021, Living Rev. Rel., 24, 1
  * Tinto et al. (2002) Tinto M., Estabrook F. B., Armstrong J. W., 2002, Phys. Rev. D, 65, 082003
  * Tinto et al. (2023) Tinto M., Dhurandhar S., Malakar D., 2023, Physical Review D, 107, 082001
  * Tiwari et al. (2023) Tiwari V., Hoy C., Fairhurst S., MacLeod D., 2023, Phys. Rev. D, 108, 023001
  * Vallisneri (2005) Vallisneri M., 2005, Phys. Rev. D, 72, 042003
  * Vecchio (2004) Vecchio A., 2004, Phys. Rev. D, 70, 042001
  * Veitch et al. (2015) Veitch J., et al., 2015, Phys. Rev. D, 91, 042003
  * Venumadhav et al. (2020) Venumadhav T., Zackay B., Roulet J., Dai L., Zaldarriaga M., 2020, Physical Review D, 101, 083030
  * Vinciguerra et al. (2017) Vinciguerra S., Veitch J., Mandel I., 2017, Class. Quant. Grav., 34, 115006
  * Weaving et al. (2023) Weaving C. R., Nuttall L. K., Harry I. W., Wu S., Nitz A., 2023, arXiv e-prints, p. arXiv:2306.16429
  * Whittle (1951) Whittle P., 1951, Uppsala: Almqvist & Wiksells Boktryckeri AB, Oxford, England
  * Williams et al. (2021) Williams M. J., Veitch J., Messenger C., 2021, Phys. Rev. D, 103, 103006
  * Wong et al. (2023) Wong K. W. K., Isi M., Edwards T. D. P., 2023, arXiv e-prints, p. arXiv:2302.05333
  * You et al. (2021) You Z.-Q., Zhu X.-J., Ashton G., Thrane E., Zhu Z.-H., 2021, Astrophys. J., 908, 215
  * Zackay et al. (2018) Zackay B., Dai L., Venumadhav T., 2018, arXiv e-prints, p. arXiv:1806.08792
  * Zackay et al. (2019a) Zackay B., Dai L., Venumadhav T., Roulet J., Zaldarriaga M., 2019a, arXiv preprint arXiv:1910.09528
  * Zackay et al. (2019b) Zackay B., Venumadhav T., Dai L., Roulet J., Zaldarriaga M., 2019b, Physical Review D, 100, 023007