Normalizing Flows for Hierarchical Bayesian Analysis:
A Gravitational Wave Population StudyDavid Ruhe
AI4Science Lab, AMLab, Anton Pannekoek Institute
University of Amsterdam
d.ruhe@uva.nl
Miles Cranmer
Department of Astrophysical Sciences
Princeton University
mcranmer@princeton.edu

Kaze Wong
Center for Computational Astrophysics
Flatiron Institutekwong@flatironinstitute.org
Patrick Forré
AI4Science Lab, AMLab
University of Amsterdam
p.d.forre@uva.nl

Abstract
We propose parameterizing the population distribution of the gravitational wave
population modeling framework (Hierarchical Bayesian Analysis) with a normalizing flow. We first demonstrate the merit of this method on illustrative experiments
and then analyze four parameters of the latest LIGO data release: primary mass,
secondary mass, redshift, and effective spin. Our results show that despite the small
and notoriously noisy dataset, the posterior predictive distributions (assuming aprior over the parameters of the flow) of the observed gravitational wave population
recover structure that agrees with robust previous phenomenological modeling
results while being less susceptible to biases introduced by less-flexible distribution
models. Therefore, the method forms a promising flexible, reliable replacement for
population inference distributions, even when data is highly noisy. 

1Introduction

The number of detected gravitational wave events is rapidly increasing [1, 31, 51, 32, 45, 2, 3].
Consequently, datasets of inferred properties such as mass, spin, and redshift have become feasible,
the most recent being produced from GWTC-3 [3]. At this point, there is sufficient data for populationlevel analyses [26, 14, 50, 39, 40, 18, 21, 22, 38, 42, 5, 4]. Many previous methods assume a simplistic
phenomenological shape of the distribution such as a power law. However, as indicated in [42], it is
crucial that the proposed model class can accurately represent the data; otherwise, the reconstruction
may miss essential population-level properties and correlations. An imposed phenomenological shape
can suffer from human-induced bias, leading to potentially incorrect conclusions. In addition, thedetection rate increases steeply with improved detector sensitivity, and it will be difficult to extract all
the information from growing datasets. [42] move away from this approach by proposing a mixture
model of weighted Gaussians (and a power law) that is expected to be capable of modeling a varietyof complex distributions.
In this work, we continue this direction, but use normalizing flows [41, 37, 33] as our model class.
The main advantage is that normalizing flows have successfully been able to represent complex
data distributions in previous works and can be fit to such distributions in a straightforward manner.
We thus rid ourselves almost entirely of the (potentially biased) constraints put previously on theinference models.
Summarizing, we show that:
Machine Learning and the Physical Sciences workshop, NeurIPS 2022.

0.4

4

2

0

2

0.3

6

0.2

4

0.1

2

0.0

4

p( |w)
Ground Truth

Densityq( |xi))p( |X)

q( |xi)

7
6
5
4
3
2
1
0

Average

4

2

0

2

4

0

Figure 1: Left: density plots of the inferred parameters q(θ|xi ) (some bold for visualization purposes)
for each event xi in our simulated dataset. Right: despite the noise in the samples, the normalizing
flow is able to recover the true population Gaussian mixture. We include the heuristic of averaging
the noisy posterior samples, yielding an incorrect result.
• In the Hierarchical Bayesian Analysis (HBA) framework (a Bayesian graphical model usedin the gravitational wave community), a normalizing flow can, despite having access tominimal and noisy data, recover structure in the data that agrees with either ground truth orrobust results from previous studies. Further, it can reveal patterns that are easily missed bycustom phenomenological modeling.
• Using a prior over the free parameters of the normalizing flow, we can get a posterior
predictive distribution that generates interpretable confidence intervals.

2

MethodLet x ∈ X ⊂ Rd be a gravitational wave time series and X := {x1 , . . . , xN } be an observed datasetof such events. We understand the physical implication of an event x through a vector of associatedparameters θ ∈ T ⊂ Rt . For binary mergers, these parameters include both masses, the redshift,and effective spin. Parameter inference for individual events comes with a lot of uncertainty due tomeasurement noise intrinsic to the underlying time series. Further, these posterior distributions canhave a non-trivial shape. This is exemplified in fig. 1 and fig. 2, where we show on the left the density
estimates of θi associated with xi for every i = 1, . . . , N . It is clear that heuristics such as taking
the mean of the samples of the posterior distributions and fitting the population model to those is
not sufficient. Therefore, the gravitational wave community has resorted to Hierarchical Bayesian
Analysis (HBA) [28, 16, 47], a Bayesian graphical model (see appendix A) for population-level
inference. By taking into account measurement uncertainty, HBA distills the posterior samples into asingle population model, which can be seen as a form of deconvolution [9]. In this way, despite huge
uncertainties in the parameters of individual events, we can still draw scientifically sound downstream
conclusions about important characteristics of the Universe.
2.1

Hierarchical Bayesian Analysis

We are interested in obtaining a population model for θ after observing X. First, we introduce amodel parameter w ∈ W ⊂ Rw with a fixed prior distribution p(w). For a fully specified generativemodel, we further need to introduce distributions p(θi |w) and p(xi |θi ), i = 1, . . . , N . We are
QN Rinterested in the posterior distribution p(w|X) ∝ p(w) · i=1 dθi p(θi |w) p(xi |θi )and the posterior predictive distributionZp(θN +1 |X) = dw p(w|X) p(θN +1 |w) = Ep(w|X) [p(θN +1 |w)] ,(1)

which expresses the probability of observing a new event with parameters θN +1 given the (old)observations X while taking into account the uncertainty in w.
While we usually do not have access to the marginal p(θi ) we can make use of a (usually physics-informed prior) distribution q(θi ), and then receive samples from the conditional2

101

p(m1|X)q(m1|xi)

100
10 1
10 2
10 3

101

GW150914
GW170817
GW190521

0

20

40

60

m1 [M ]

80

100
10 1
10 2
10 3

100

0

20

40

60

m1 [M ]

80

100

Figure 2: Left: the GWTC-3 posterior density plots. Right: ourRRRinferred population model using
normalizing flows for the inspiral primary mass p(m1 |X) =dm1 dz dXeff p(θ|X) where weused a Monte Carlo estimate of eq. (1) to generate a 90% percentile interval.
R|θi ) q(θi )θi ∼ q(θi |xi ) := p(xiq(x, where q(xi ) :=dθi p(xi |θi ) q(θi ). I.e., we assume thati)
we have access to samples from the individual posterior distributions. As such, our dataset is
D := {{θ11 , . . . , θ1M1 } , . . . , {θN 1 , . . . , θN MN }} , where for every event i ∈ (1, . . . , N ) we have
Mi samples from q(θi |xi ). We are first interested in maximizing an approximate posterior using D:



Mi 
N
Y
X
p(w)
1
p(θ
|w)
ij
,

(2)
w⋆ := arg max p̃(w|X),p̃(w|X) :=w∈W
Z(X, D) i=1 Mi j=1 q(θij )with normalizing constant Z(X, D) (derivation appendix B). A similar Monte-Carlo-based estimatoris presented in [44]. After obtaining w⋆ through, e.g., stochastic gradient descent, we run Markov
chain Monte-Carlo around the MAP solution to obtain samples w ∼ p̃(w|X). These samples and ourmodel p(θ|w) can be used to approximate the posterior predictive distribution eq. (1). It is important
to model the whole posterior distribution, as the uncertainty in the parameters tells us where we candraw robust conclusions about the underlying population.
2.2

Normalizing Flows

So far, we have not discussed how we parameterize p(θ|w). To facilitate increasingly larger datasetsand to be able to model higher-dimensional distributions, the goal of this research is to explore more
flexible and scalable models. As such, we propose normalizing flows [41, 37, 33, appendix C] as asuitable model class. Normalizing flows transform a density using a chain of a smooth, invertible
mappings f : Rt → Rt such that θK = fK ◦ · · · ◦ f2 ◦ f1 (θ0 ), which constructs an arbitrarily
complex density. We parameterize each fk using a set of weights w, and construct a log-likelihoodPK∂f w
log p(θK |w) = log p0 (θ0 ) − k=1 log det ∂θkk , which we can optimize for w using, e.g., stochastic gradient descent. In this work, we use planar flows [37] and block neural autoregressive flows [12],as we experimentally found that they are sufficiently flexible while not using too many parameters.
Furthermore, we found that the resulting distributions are smooth, which makes them physically moreplausible and allows for better automated discovery of interpretable models [49].

3Results

3.1

Illustrative Examples

In an illustrative example, we use a mixture of Gaussians as the ground truth marginal distribution overθ. The observable distribution is a Gaussian with mean |θ|, from which we draw X. Using a uniform
prior q(θ), we obtain posterior draws θ ∼ q(θ|xi ) for all events using rejection sampling. Due tothe noninvertible forward model, the posteriors are bimodal and noisy. Since the true underlying
population distribution and the posteriors are nontrivial, standard techniques will not be effective here.
Still, the normalizing flow can, using the Hierarchical Bayesian Analysis framework, convert the
noisy posterior draws shown in fig. 1 (left) into a faithful population model as shown in fig. 1 (right).
3

Further, we perform a more realistic (though less clearly visualizable) experiment with synthetic
data in appendix D.2 where we sample data from previously suggested gravitational wave population
models (e.g., by [4]) and recover the proposed distributions.
3.2

GWTC-3: Primary Mass Marginal

3.3

0.5
1.0
1.5
2.0
0.8
0.4
0.0
0.4
0.8

0.2
0.4
0.6
0.8

20
40
60
80

0.8 0.4 0.0 eff0.4 0.8

z

0.5 1.0 1.5 2.0

q

0.2 0.4 0.6 0.8

In our experiments on the GWTC-3 catalog, we
consider four gravitional wave parameters: m1(primary mass), m2 (secondary mass), z (redshift) and Xeff (effective inspiral spin). I.e., θ :=(m1 , m2 , z, Xeff ) and t = 4. Note that the number of recorded gravitational wave events is currently rather small (N = 76). Still, we show that
we can retrieve structure from the data that agrees
with previous studies. After getting p(θ|w⋆ ) by
maximizing eq. (2), we obtain, using Hamiltonian Monte Carlo [30], L = 10, 000 samples
w ∼ p(w|θ), starting the chain at w∗ . Next, we
obtain the posterior predictive distribution by theapproximation p(θ|X) ≈ Ep(w|X) [p(θ|w)] ≈PL1qzm1 [M ]eff
l=1 p(θ|wl ). We show the resulting obL
served populationRRR
posterior predictive marginalp(mdm2 dz dXeff p(θ|X) = Figure 3: Pair plot of joint marginals of our inRRR 1 |X) =dm2 dz dXeff p(m1, m2, z, Xeff |X) in fig. 2. ferred p(θ|X) Note that we did not take in toNote, again, that we did not impose any strong account selection bias effects which is especiallyrestrictions on the model class; this result was noticeable in the redshift marginal. Hence, this is
fully recovered from just the data. Overall, we an observed population analysis.
see (apart from selection bias effects [28]) similarity with, e.g., the P OWER L AW + P EAK and M ULTI -P EAK model of [5, 4]. Specifically, weobserve distinct peaks around m1 ∼ 10 M⊙ , m1 ∼ 40 M⊙ , and m1 ∼ 60 M⊙ . Finally, we alsoreassuringly observe that the uncertainty contracts or expands in the parts where we expect it to. I.e.,
it expands at regions where we do not have many observations, and contracts in the regions where wedo (e.g., around the high-mass modes).
GWTC-3: Observed Population Distribution

m2
We display the marginal joint distributions in fig. 3. Here, we define q := mas the ratio of secondary
1
to primary mass of the binary. Note this figure shows the observed population results. This isespecially noticeable in the redshift marginal z, which is heavily biased towards low-redshift events.
Further, we expect more lower primary mass events m1 . Xeff and q are expected to be less affectedby not including the selection function. To draw conclusions about the true population, one should
include selection bias effects. Optimization was highly unstable when including the selection functionin the objective. Hence, including selection bias effects was left for future investigation. Also, q is
expected to peak at q = 1, with an infinitely steep cutoff at q > 1. We did not observe this due to thesmoothness of the flow. Future work could look into reparameterizations to enable this behavior.

We can tentatively draw a few conclusions from the observed population figure. Interestingly, itsuggests that m1 and q do not simply correlate through a power law. Instead, the distribution shows
an inverted V shape. The redshift to mass correlations seem plausible and consistent. The effective
spin marginal also aligns with previous literature. In [29] is was noted that q does not correlate
positively with Xeff , which is also confirmed here.

4

Discussion

We parameterized the population model in the Hierarchical Bayesian Analysis framework with anormalizing flow. In illustrative experiments, we confirmed that the model can retrieve the actual
population density under nontrivial posterior samples. On gravitational wave data from GWTC-3, weconsidered primary mass, secondary mass, redshift, and effective spin and recovered an observed
4population model that agrees with robust previous phenomenological modeling results despite thedataset being small and highly noisy. This paves the way for less constrained, tractable, and flexible
population-level inference in noisy settings, especially once more higher-dimensional data become
available, allowing for automated discovery of structure in the data.

5

Broader Impact

This research proposed to combine normalizing flows with the HBA framework. The current work
found an application to gravitational wave event data, but in principle, the method can be implemented
in many fields of science. While this work purely aims to aid scientific discovery, the approach canalso be used in non-scientific settings dealing with noisy observations, which can pose privacy-related
issues depending on the data type and source. As such, we encourage more data privacy awarenessand establishing data protection policies.

6

Acknowledgements

This work was largely performed as part of the machine learning summer program of the Flatiron
Institute. We would like to thank the scientific software development community, without whom thiswork would not be possible [20, 43, 36, 7, 35, 46, 11, 48].

References[1] BP Abbott, Richard Abbott, TDea Abbott, S Abraham, F Acernese, K Ackley, C Adams,
RX Adhikari, VB Adya, Christoph Affeldt, et al. Gwtc-1: a gravitational-wave transient catalog
of compact binary mergers observed by ligo and virgo during the first and second observing
runs. Physical Review X, 9(3):031040, 2019.
[2] R Abbott, TD Abbott, S Abraham, F Acernese, K Ackley, A Adams, C Adams, RX Adhikari,
VB Adya, Christoph Affeldt, et al. Gwtc-2: compact binary coalescences observed by ligo andvirgo during the first half of the third observing run. Physical Review X, 11(2):021053, 2021.
[3] R Abbott, TD Abbott, F Acernese, K Ackley, C Adams, N Adhikari, RX Adhikari, VB Adya,C Affeldt, D Agarwal, et al. Gwtc-3: compact binary coalescences observed by ligo and virgo
during the second part of the third observing run. arXiv preprint arXiv:2111.03606, 2021.
[4] R Abbott, TD Abbott, F Acernese, K Ackley, C Adams, N Adhikari, RX Adhikari, VB Adya,C Affeldt, D Agarwal, et al. The population of merging compact binaries inferred using
gravitational waves through gwtc-3. arXiv preprint arXiv:2111.03634, 2021.
[5] Rich Abbott, TD Abbott, S Abraham, F Acernese, K Ackley, A Adams, C Adams, RX Adhikari,
VB Adya, Christoph Affeldt, et al. Population properties of compact objects from the second
ligo–virgo gravitational-wave transient catalog. The Astrophysical journal letters, 913(1):L7,
2021.
[6] Rianne van den Berg, Leonard Hasenclever, Jakub M Tomczak, and Max Welling. Sylvester
normalizing flows for variational inference. arXiv preprint arXiv:1803.05649, 2018.
[7] Eli Bingham, Jonathan P Chen, Martin Jankowiak, Fritz Obermeyer, Neeraj Pradhan, Theofanis
Karaletsos, Rohit Singh, Paul Szerlip, Paul Horsfall, and Noah D Goodman. Pyro: Deep
universal probabilistic programming. The Journal of Machine Learning Research, 20(1):973–
978, 2019.
[8] Charles Blundell, Julien Cornebise, Koray Kavukcuoglu, and Daan Wierstra. Weight uncertainty
in neural network. In International conference on machine learning, pages 1613–1622. PMLR,
2015.
[9] Jo Bovy, David W Hogg, and Sam T Roweis. Extreme deconvolution: Inferring complete
distribution functions from noisy, heterogeneous and incomplete observations. The Annals ofApplied Statistics, 5(2B):1657–1677, 2011.
[10] Thomas A Callister, Carl-Johan Haster, Ken KY Ng, Salvatore Vitale, and Will M Farr. Whoordered that? unequal-mass binary black hole mergers have larger effective spins. The Astrophysical Journal Letters, 922(1):L5, 2021.
5

[11] Adam D Cobb and Brian Jalaian. Scaling hamiltonian monte carlo inference for bayesian neuralnetworks with symmetric splitting. Uncertainty in Artificial Intelligence, 2021.
[12] Nicola De Cao, Wilker Aziz, and Ivan Titov. Block neural autoregressive flow. In Uncertainty
in artificial intelligence, pages 1263–1273. PMLR, 2020.
[13] Laurent Dinh, Jascha Sohl-Dickstein, and Samy Bengio. Density estimation using real nvp.
arXiv preprint arXiv:1605.08803, 2016.
[14] Maya Fishbach, Daniel E Holz, and Will M Farr. Does the black hole merger rate evolve with
redshift? The Astrophysical Journal Letters, 863(2):L41, 2018.
[15] Maya Fishbach, Daniel E. Holz, and Will M. Farr. Does the Black Hole Merger Rate Evolvewith Redshift? Astrophys. J. Lett., 863(2):L41, 2018.
[16] Sebastian M Gaebel, John Veitch, Thomas Dent, and Will M Farr. Digging the population ofcompact binary mergers out of the noise. Monthly Notices of the Royal Astronomical Society,484(3):4008–4023, 2019.
[17] Yarin Gal and Zoubin Ghahramani. Dropout as a bayesian approximation: Representing model
uncertainty in deep learning. In international conference on machine learning, pages 1050–1059.
PMLR, 2016.
[18] Shanika Galaudage, Colm Talbot, and Eric Thrane. Gravitational-wave inference in the catalog
era: Evolving priors and marginal events. Physical Review D, 102(8):083026, 2020.
[19] Emiel Hoogeboom, Victor Garcia Satorras, Jakub Tomczak, and Max Welling. The convolution
exponential and generalized sylvester flows. Advances in Neural Information Processing
Systems, 33:18249–18260, 2020.
[20] John D Hunter. Matplotlib: A 2d graphics environment. Computing in science & engineering,9(03):90–95, 2007.
[21] Shasvath J Kapadia, Sarah Caudill, Jolien DE Creighton, Will M Farr, Gregory Mendell, Alan
Weinstein, Kipp Cannon, Heather Fong, Patrick Godwin, Rico KL Lo, et al. A self-consistentmethod to estimate the rate of compact binary coalescences with a poisson mixture model.
Classical and Quantum Gravity, 37(4):045007, 2020.
[22] Chase Kimball, Colm Talbot, Christopher PL Berry, Matthew Carney, Michael Zevin, Eric
Thrane, and Vicky Kalogera. Black hole genealogy: Identifying hierarchical mergers withgravitational waves. The Astrophysical Journal, 900(2):177, 2020.
[23] Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint
arXiv:1412.6980, 2014.
[24] Durk P Kingma and Prafulla Dhariwal. Glow: Generative flow with invertible 1x1 convolutions.
Advances in neural information processing systems, 31, 2018.
[25] Durk P Kingma, Tim Salimans, Rafal Jozefowicz, Xi Chen, Ilya Sutskever, and Max Welling.
Improved variational inference with inverse autoregressive flow. Advances in neural information
processing systems, 29, 2016.
[26] Ely D Kovetz, Ilias Cholis, Patrick C Breysse, and Marc Kamionkowski. Black hole massfunction from gravitational wave measurements. Physical Review D, 95(10):103010, 2017.
[27] Balaji Lakshminarayanan, Alexander Pritzel, and Charles Blundell. Simple and scalable
predictive uncertainty estimation using deep ensembles. Advances in neural information
processing systems, 30, 2017.
[28] Ilya Mandel, Will M Farr, and Jonathan R Gair. Extracting distribution parameters from multiple
uncertain observations with selection biases. Monthly Notices of the Royal Astronomical Society,486(1):1086–1093, 2019.
[29] B McKernan, KES Ford, T Callister, WM Farr, R O’Shaughnessy, R Smith, E Thrane, andA Vajpeyi. Ligo–virgo correlations between mass ratio and effective inspiral spin: testing theactive galactic nuclei channel. Monthly Notices of the Royal Astronomical Society, 514(3):3886–
3893, 2022.
[30] Radford M Neal et al. Mcmc using hamiltonian dynamics. Handbook of markov chain montecarlo, 2(11):2, 2011.
6

[31] Alexander H Nitz, Collin Capano, Alex B Nielsen, Steven Reyes, Rebecca White, Duncan A
Brown, and Badri Krishnan. 1-ogc: The first open gravitational-wave catalog of binary mergers
from analysis of public advanced ligo data. The Astrophysical Journal, 872(2):195, 2019.
[32] Alexander H Nitz, Thomas Dent, Gareth S Davies, Sumit Kumar, Collin D Capano, Ian
Harry, Simone Mozzon, Laura Nuttall, Andrew Lundgren, and Márton Tápai. 2-ogc: Opengravitational-wave catalog of binary mergers from analysis of public advanced ligo and virgo
data. The Astrophysical Journal, 891(2):123, 2020.
[33] George Papamakarios, Eric T Nalisnick, Danilo Jimenez Rezende, Shakir Mohamed, and Balaji
Lakshminarayanan. Normalizing flows for probabilistic modeling and inference. J. Mach.
Learn. Res., 22(57):1–64, 2021.
[34] George Papamakarios, Theo Pavlakou, and Iain Murray. Masked autoregressive flow for density
estimation. Advances in neural information processing systems, 30, 2017.
[35] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan,
Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al. Pytorch: An imperative
style, high-performance deep learning library. Advances in neural information processing
systems, 32, 2019.
[36] Fabian Pedregosa, Gaël Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion,
Olivier Grisel, Mathieu Blondel, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, et al. Scikitlearn: Machine learning in python. the Journal of machine Learning research, 12:2825–2830,
2011.
[37] Danilo Rezende and Shakir Mohamed. Variational inference with normalizing flows. In
International conference on machine learning, pages 1530–1538. PMLR, 2015.
[38] Javier Roulet, Tejaswi Venumadhav, Barak Zackay, Liang Dai, and Matias Zaldarriaga. Binary
black hole mergers from ligo/virgo o1 and o2: Population inference combining confident andmarginal events. Physical Review D, 102(12):123022, 2020.
[39] Javier Roulet and Matias Zaldarriaga. Constraints on binary black hole populations fromligo–virgo detections. Monthly Notices of the Royal Astronomical Society, 484(3):4216–4229,
2019.
[40] Rory JE Smith, Colm Talbot, Francisco Hernandez Vivanco, and Eric Thrane. Inferring the
population properties of binary black holes from unresolved gravitational waves. Monthly
Notices of the Royal Astronomical Society, 496(3):3281–3290, 2020.
[41] Esteban G Tabak and Eric Vanden-Eijnden. Density estimation by dual ascent of the loglikelihood. Communications in Mathematical Sciences, 8(1):217–233, 2010.
[42] Vaibhav Tiwari. Vamana: modeling binary black hole population with minimal assumptions.
Classical and Quantum Gravity, 38(15):155007, 2021.
[43] Stefan Van Der Walt, S Chris Colbert, and Gael Varoquaux. The numpy array: a structure for
efficient numerical computation. Computing in science & engineering, 13(2):22–30, 2011.
[44] Maxime Vandegar, Michael Kagan, Antoine Wehenkel, and Gilles Louppe. Neural empiricalbayes: Source distribution estimation and its applications to simulation-based inference. In
International Conference on Artificial Intelligence and Statistics, pages 2107–2115. PMLR,
2021.
[45] Tejaswi Venumadhav, Barak Zackay, Javier Roulet, Liang Dai, and Matias Zaldarriaga. New
binary black hole mergers in the second observing run of advanced ligo and advanced virgo.
Physical Review D, 101(8):083030, 2020.
[46] Pauli Virtanen, Ralf Gommers, Travis E Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, et al. Scipy 1.0:
fundamental algorithms for scientific computing in python. Nature methods, 17(3):261–272,
2020.
[47] Salvatore Vitale, Davide Gerosa, Will M Farr, and Stephen R Taylor. Inferring the properties ofa population of compact binaries in presence of selection effects. In Handbook of Gravitational
Wave Astronomy, pages 1–60. Springer, 2022.
[48] Michael L Waskom. Seaborn: statistical data visualization. Journal of Open Source Software,
6(60):3021, 2021.
7

[49] Kaze W. K Wong and Miles Cranmer. Automated discovery of interpretable gravitational-wave
population models, 2022.
[50] Daniel Wysocki, Jacob Lange, and Richard O’Shaughnessy. Reconstructing phenomenologicaldistributions of compact binaries via gravitational wave observations. Physical Review D,100(4):043012, 2019.
[51] Barak Zackay, Tejaswi Venumadhav, Liang Dai, Javier Roulet, and Matias Zaldarriaga. Highly
spinning and aligned binary black hole merger in the advanced ligo first observing run. PhysicalReview D, 100(2):023007, 2019.

Checklist
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper’s
contributions and scope? [Yes](b) Did you describe the limitations of your work? [Yes] E.g., the exclusion of selection
bias effects.
(c) Did you discuss any potential negative societal impacts of your work? [Yes] In the
broader impact statement.
(d) Have you read the ethics review guidelines and ensured that your paper conforms tothem? [Yes]2. If you are including theoretical results...
(a) Did you state the full set of assumptions of all theoretical results? [Yes](b) Did you include complete proofs of all theoretical results? [Yes]3. If you ran experiments...
(a) Did you include the code, data, and instructions needed to reproduce the main experimental results (either in the supplemental material or as a URL)? [Yes] Code will be
released upon acceptance.
(b) Did you specify all the training details (e.g., data splits, hyperparameters, how they
were chosen)? [Yes] In the supplementary material.
(c) Did you report error bars (e.g., with respect to the random seed after running experiments multiple times)? [Yes](d) Did you include the total amount of compute and the type of resources used (e.g., typeof GPUs, internal cluster, or cloud provider)? [Yes]4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets...
(a) If your work uses existing assets, did you cite the creators? [Yes] Will be done inacknowledgements statement.
(b) Did you mention the license of the assets? [N/A](c) Did you include any new assets either in the supplemental material or as a URL? [N/A](d) Did you discuss whether and how consent was obtained from people whose data you’reusing/curating? [N/A](e) Did you discuss whether the data you are using/curating contains personally identifiable
information or offensive content? [N/A]
5. If you used crowdsourcing or conducted research with human subjects...
(a) Did you include the full text of instructions given to participants and screenshots, if
applicable? [N/A](b) Did you describe any potential participant risks, with links to Institutional Review
Board (IRB) approvals, if applicable? [N/A](c) Did you include the estimated hourly wage paid to participants and the total amountspent on participant compensation? [N/A]8

xN +1
xi

θN +1

θi
N

w
Figure 4: Graphical model of Hierarchical Bayesian Analysis. The red arrows are the ones that ournormalizing flow parameterizes.

A

Graphical Model

The graphical we use is shown in fig. 4. We colored the arrows that our normalizing flow parameterizes
red. To estimate p(θN +1 | X), samples θi ∼ q(θi | xi ) are inferred for all xi in the first stage. In
this work we assume we have those samples and use them to estimate the remaining densities in
eq. (1) through the objective eq. (2).

B

Hierarchical Bayesian Likelihood|θi ) q(θi )Using q(θi |xi ) := p(xiq(x, we set
i)



q(x)
p̂(x|w) := dθ q(θ|x) ·p(θ|w)q(θ)
p(θ|w).
= q(x) · Eq(θ|x)q(θ)
Note that q(xi ) is a constant with respect to our model parameters w. Hence,Np(w) Yp(θ|w)
p̂(w|X) :=Eq(θ|xi )Z(X, D) i=1
q(θ)



M 
N
p(w) Y  1 Xi p(θj |w) =: p̃(w|X)≈
Z(X, D) i=1 Mi j=1 q(θj )
Z

C

(3)
(4)

(5)

(6)

Normalizing Flows

To facilitate increasingly larger datasets and to be able to model higher-dimensional distributions, thegoal of this research is to explore more flexible and scalable models. As such, we propose normalizing
flows[41, 37, 33] as a suitable model class. Normalizing flows transform a density using a smooth,
invertible mapping f : Rt → Rt with inverse g := f −1 . As such, if we transform a random variable
θ, then θ ′ has∂f −1p(θ ′ ) = p(θ) det
.
(7)
∂θ ′
If we use a chain of these transformations, such thatθK = fK ◦ · · · ◦ f2 ◦ f1 (θ0 ),(8)
we can construct an arbitrarily complex density. We parameterize each fk using a set of weights w,
and construct a log-likelihood
KX∂fkw
log p(θK |w) = log p0 (θ0 ) −log det
,
(9)
∂θk
k=1

9

which we can optimize for w using, e.g., stochastic gradient descent. The downside is a potentiallycomputationally expensive Jacobian determinant. As such, much work has been conducted on findingtransformations that have computationally cheap Jacobian determinants [13, 25, 34, 6, 24, 19]. In
this work, we use planar flows [37] and block neural autoregressive flows [12], as we experimentallyfound that they are sufficiently flexible while not using too many parameters. Furthermore, we foundthat the resulting distributions are (relative to some other parameterizations) smooth, which makes
them physically more plausible.
Planar Flow Planar flow uses the transformation
fkw (θ) := θ + uk h(vk⊤ θ + bk ),(10)

where (uk , zk , bk ) ∈ w are free parameters, and h is the hyperbolic tangent. These transformations
contract and expand the inputs perpendicular to the plane defined vk⊤ θ − bk = 0. It can be shownthat the Jacobian determinant of this mapping can be computed in O(t) [37]. As such, we can use
a sequence of these transformations to model p(θ|w), our gravitational wave population model(eq. (2)).
Block Neural Autoregressive Flow Neural autoregressive flows decomposes a joint distribution
over θ ∈ Rt into t conditional distributions. the transformations fkw then yield lower triangular
Jacobians and hence cheaply computable determinants. Block neural autoregressive flows thenparameterize fkw using a neural network that uses block lower-triangular matrices. Their diagonal
elements are strictly positive, ensuring monotonicity. The off-diagonal elements do not requiremonotonicity as they only play a role in the conditioning part of the conditional distributions.

D

Experimental Details

All experiments were conducted locally.
D.1

Illustrative Example

The ground-truth population is defined by a mixture θ ∼ N (−0.75, 0.52 )N (0.75, 0.52 ). The forward
model is given by x|θ ∼ N (|θ|, 0.12 ). This noninvertible forward mapping results in noisy, bimodal
posteriors. We sample N = 1024 events. Using a uniform prior q(θ) := U (−5, 5) and the true
forward model, we draw Mi = 128 posterior samples using rejection sampling from q(θ|xi ) for
every xi , i ∈ 1, . . . , 1024. We parameterize p(θ|w) as a 4-layer planar flow and a standard Gaussianbase distribution. We use the Adam optimizer [23] with learning rate 0.01 so obtain w⋆ . We foundempirically that optimization runs better when mini-batches are taken for both the events and also theposterior draws. In this experiment, we batch sizes of N ′ := 1024 events and Mi′ := 4.
D.2

Synthetic Data Experiment

Here, we simulate a population model using previously published models and investigate whether
our flow-based approach can faithfully recover those distributions. We start with the P OWER L AW +P EAK model for p(m1 , m2 ) from [4] using the means of the posterior distributions for ground-truth
population parameters. For a conditional effective spin distribution we use [10], again with the meanof the inferred distributions as ground-truth parameter values. Finally, we model the distributionof redshift using a power law with spectral index κ = 1.7 following [15]. Note that there is an
extra factor in the merger coming from the comoving volume factor. The resulting samples areshown in fig. 5 (left). We simulate Gaussian posteriors in our Hierarchical Bayesian framework. The
distribution that the normalizing flow recovers is shown on the right in fig. 5. We see that the flow
recovers most of the correlations and shapes of the joint distribution faithfully.
We reparameterize the m1 samples to log-space and q samples to logit-space. We sample N = 1024
events and use 32 hierarchical samples per batch. We used the Adam optimizer [23] with learning
rate of 0.001 and 1 layer of block neural autoregressive flow with 2 8-dimensional hidden layers.
100.2 0.0 0.2eff0.5
5 0 5 0q

0.2 0.4 0.6 0.8

qeff
0.0 0.2 0.4 0.6 0.2 0.4 0.6 0.8

q

eff

2
4
6
8

m1 [M ]

0.2
0.4
0.6
0.8

2
4
6
8

z

0.2
0.05
0
0.2
5
0.5
0

eff

15
30
45

q

0.0
0.2
0.4
0.6

0.2
0.4
0.6
0.8

20
40
60
80

z

z

2 4 6 82 4 6 8m1 [M ]z

Figure 5: Left: simulated gravitational wave data. Right: recovered distributions by the normalizing
flow.
D.3

GWTC-3We have N = 76 gravitational wave events available. The number of posterior samples per event
ranges from Mi = 3194 to Mi = 268806. We use 10 layers of planar flow, resulting in a total of
w = 90 free parameters, limiting the risk of severe over-fitting resulting in unnatural population
distribution shapes. We use mini-batches of N ′ := 76 events and Mi′ := 32 posterior samples per
event. We use the Adam optimizer [23] with learning rate 0.001. to find w⋆ .

E

Future WorkIn this study, we did not include selection bias effects. Future work should adjust for this to obtain
more realistic models of population parameters of gravitational waves. We also expect that we did
not fully explore all parameterizations (and hyperparameters) that recover the population model
as faithfully as possible, both in terms of the normalizing flow and how the input data should betransformed for the flow to optimally learn the densities.
Further, [49] have presented a promising approach to distill black-box models into symbolic equations.
Such an approach can also be applied to our normalizing flows, yielding more interpretable versionsof the presented population densities.
Finally, we used Hamiltonian Monte-Carlo to get samples w ∼ p(w|θ). This is not scalable when thenumber of weights increases in order to be able to fit more complex data. To resolve this, follow-up
works could explore approaches based on the Bayesian Neural Network literature, e.g., [8, 17, 27].

11