# DFTS2: Simulating Deep Feature Transmission Over Packet Loss Channels
††thanks: This work was supported in part by the NSERC grant RGPIN-2021-02485.

Ashiv Dhondea, Robert A. Cohen, and Ivan V. Bajić  School of Engineering
Science, Simon Fraser University, Burnaby, BC, V5A 1S6, Canada

###### Abstract

In edge-cloud collaborative intelligence (CI), an unreliable transmission
channel exists in the information path of the AI model performing the
inference. It is important to be able to simulate the performance of the CI
system across an imperfect channel in order to understand system behavior and
develop appropriate error control strategies. In this paper we present a
simulation framework called DFTS2, which enables researchers to define the
components of the CI system in TensorFlow 2, select a packet-based channel
model with various parameters, and simulate system behavior under various
channel conditions and error/loss control strategies. Using DFTS2, we also
present the most comprehensive study to date of the packet loss concealment
methods for collaborative image classification models.

###### Index Terms:

Collaborative intelligence, distributed AI, feature transmission, error
resilience, packet loss

##  I Introduction

Collaborative intelligence (CI) [1, 2] is a “hybrid” AI deployment strategy.
To complete an AI task, CI typically involves splitting the Deep Neural
Network (DNN) computations between an edge device and a cloud service. The
resource-constrained edge device runs the initial few layers of the DNN on the
input signal (image, video, …) and produces a deep feature tensor. The latter
is transmitted over a communication channel to the cloud, where the remaining
layers of the DNN complete the inference task. CI has shown potential for
improved inference latency and energy use as compared to purely cloud-based or
edge-based AI deployment approaches in certain scenarios [1, 3].

Potential applications of CI, which include IoT systems, traffic monitoring,
autonomous vehicles or drones and point-of-care systems, may involve imperfect
communication channels. This motivates the need to investigate the performance
of the CI system across such channels in order to understand system behavior
and develop error/loss resilience strategies. To this end, DFTS – Deep Feature
Transmission Simulator [4] – was developed to study the effects of missing
feature data on inference accuracy in CI. Here we present a more sophisticated
simulation framework called DFTS2, with several new capabilities: additional
channel models and missing feature recovery methods from the recent
literature.

This paper is organized as follows. Section II discusses related work on
packet loss concealment in CI. Section III describes the operation of the
proposed simulation framework. Experiments are discussed in Section IV
followed by conclusions in Section V. An implementation of DFTS2 is available
online, along with packet traces and simulation files for results presented in
this work.111https://github.com/AshivDhondea/DFTS2

##  II Related work

Two approaches to error resilience in CI have been studied in the literature:
1) joint source-channel coding of deep features [5, 6] and 2) packet-based
feature transmission with simple error concealment, such as bilinear and
nearest-neighbor interpolation in [4], tensor completion [7, 8], and
inpainting [9].

Generic tensor completion methods such as Simple Low Rank Tensor Completion
(SiLRTC) [10] and High Accuracy Low Rank Tensor Completion (HaLRTC) [10] were
adapted to recover corrupted deep feature data packets in [7, 8, 9].
Furthermore, [7] proposed a method called ALTeC - Adaptive Linear Tensor
Completion - specifically designed to recover missing intermediate features by
pre-training on deep feature tensors. With a single-feature-row-per-packet
scheme over a random loss channel, ALTeC proved to be as least as accurate as
SiLRTC and HaLRTC while bringing a significant reduction in inference latency
[7]. In [8], a Content Adaptive Linear Tensor Completion (CALTeC) method was
proposed to exploit existing redundancies in deep feature tensors. With a
multiple-feature-rows-per-packet scheme over a Gilbert-Elliott channel, CALTeC
and HaLRTC produced similar Top-1 classofication accuracy in [8]. ALTeC
slightly outperformed CALTeC and HaLRTC in one scenario and offered slightly
worse performance than the other two in another scenario when the split point
of the DNN was closer to its input and the number of feature rows per packet
was increased. SiLRTC-completed tensors led to the poorest Top-1
classification accuracy in both scenarios in [8]. In [9], PDE-based image
inpainting methods, “Navier-Stokes” [11] and “Telea” [12], performed lost
packet concealment in collaborative object detection. Both inpainting methods
led to higher mean average precision than SiLRTC with 250 iterations, while
being more than 100$\times$ faster.

This paper presents an extended version of DFTS [4], called DFTS2, with more
advanced capabilities for simulating packet loss and loss recovery.
Specifically, the Gilbert-Elliot channel model [13] and packet traces are
introduced. Moreover, the recent feature recovery mechanisms mentioned above -
SiLRTC, HaLRTC, ALTeC, CALTeC, and feature inpainting - have been implemented.
Using DFTS2, we present a comprehensive comparison of the mentioned feature
recovery methods on collaborative image classification.

##  III DFTS2 simulator description

DFTS2 splits TensorFlow 2 DNN models into a mobile (or edge) sub-model and a
cloud (or remote) sub-model at a user-defined split point. After processing an
input image, the mobile sub-model produces an intermediate feature tensor,
which is transmitted over an imperfect communication channel to the remote
service. Prior to transmission over a packet network, the feature tensor is
packetized and quantized. Fig. 1 gives a system-level overview of a DFTS2
experiment. A lossy communication channel between the mobile device and the
remote service leads to data packets being damaged or lost due to errors in
data transmission and network congestion, for example. The different
components involved in a DFTS2 experiment are summarized in Table I and
discussed in detail in the following subsections.

Figure 1: DFTS2 system overview (Adapted from [14, Fig. 1])

###  III-A Intermediate feature tensor packetization

A feature tensor with height $h$, width $w$, and the number of channels $c$
will be denoted $\mathcal{X}\in\mathbb{R}^{h\times w\times c}$. An
illustration of a feature tensor extracted from layer add_1 of ResNet-18 [15]
is shown in Fig. 2. Layer add_1 corresponds to the output of the shortcut
connection and element-wise addition applied to the final output of the
conv2_x layers, as denoted in [15, Table 1]. Currently, there is no “standard”
approach to packetizing deep feature tensor data. DFTS2 employs a
packetization scheme inspired by a common practice in video streaming, where
video frames are encoded into packets of consecutive rows of frame elements
[16]. Similarly, each channel in a DFTS2 tensor is partitioned into groups of
$r_{p}$ rows, as indicated by the dashed green lines in Fig. 2. Zero padding
is done on the last packet of a channel if the number of channel rows ($h$) is
not an integer multiple of the number of packet rows ($r_{p}$). The user-
defined $r_{p}$ was set to 8 rows in our experiments. Packetized tensors were
min-max quantized [17] to 8 bits/element. While it is possible to further
compress the data by additional prediction, quantization, entropy coding [18,
19, 20, 21], and/or bit allocation [22], this was not done in our experiments
because the focus is on error resilience, not on compression.

Figure 2: A tensor $\mathcal{X}\in\mathbb{R}^{h\times w\times c}$ from layer
add_1 of ResNet-18. Several consecutive rows in each channel form a feature
data packet.

###  III-B Channel models

Four options for the communication channel model are offered in DFTS2
experiments: perfect channel (no loss), independent and identically
distributed (iid) random loss, Gilbert-Elliott (GE) channel [13], and user-
provided packet traces. The GE model, also known as the burst loss model, has
been shown to adequately capture observed packet loss patterns over the
Internet [13, 23]. This motivates our choice of GE channel models for our
experiments. This packet loss model is a two-state Markov model with a good
$(G)$ and a bad $(B)$ state. It can be parameterized by the burst loss
probability $P_{B}$ and the average burst length $L_{B}$. The conversion from
$P_{B}$ and $L_{B}$ to the four state transition probabilities is as follows:

$\displaystyle p_{B\to G}$ $\displaystyle=\frac{1}{L_{B}},\quad$
$\displaystyle p_{G\to B}=\frac{P_{B}}{L_{B}(1-P_{B})},$ (1) $\displaystyle
p_{B\to B}$ $\displaystyle=1-p_{B\to G},\quad$ $\displaystyle p_{G\to
G}=1-p_{G\to B}.$

Since DFTS2 can also load packet traces from imported files to create
corrupted deep feature tensor packets, it allows the user to provide packet
traces from an external source, such as from actual measurements over a
network or from a communication network simulation package.

When a feature tensor is transmitted, some of its packets will be lost,
depending on the channel realization. The remote client identifies missing
packets from packet sequence numbers provided by a transport-layer protocol
and reconstitutes the available part of the tensor. The available data is
inverse quantized, and the missing data (identified by a Boolean loss map) is
recovered using an error concealment method.

###  III-C Packet loss concealment

Missing value recovery in corrupted deep feature tensors means that more
accurate tensor data can be exploited by the back-end sub-model, resulting in
significant gains in classification accuracy [8].

SiLRTC [10] and HaLRTC [10] assume that the original, un-corrupted tensors are
low-rank. They iteratively reduce the rank of the damaged tensor via singular
value shrinkage to try to recover the original tensor values. Due to their
iterative operation, these general tensor completion methods introduce
significant processing delays, which is a problem in applications sensitive to
inference latency. SiLRTC and HaLRTC make no assumption about the spatial
distribution of missing values in a corrupted tensor, or how the tensor was
created in the first place. In [7], the Frobenius norm of the difference
between completed tensors over successive iterations of tensor completion with
SiLRTC and HaLRTC was found to have reached a plateau within 50 iterations.
Therefore, we ran SiLRTC and HaLRTC for 50 iterations in our experiments.

Besides SiLRTC and HaLRTC, DFTS2 also includes ALTeC [7], CALTeC [8], and
Navier-Stokes inpainting [11]. All these methods exploit spatial structure of
feature tensor channels to recover missing data. ALTeC assumes a linear
relationship between a packet of feature data and its spatial neighbors, as
well as collocated data in other feature channels. The coefficients of this
linear relationship are obtained by pre-training. Meanwhile, for each missing
feature packet, CALTeC finds the most correlated feature channel where
collocated data is available, and then derives an affine transformation to
fill in the missing data. Navier-Stokes inpainting is based on a “surface
flow” of available spatially neighboring data into the missing region.

###  III-D Simulation mode

The simulator can run in a “single-shot” mode or a Monte Carlo simulation
mode. In the single-shot mode, a single image is input to the edge sub-model,
the transmission of the corresponding feature packets is simulated, and the
received features are processed by the cloud sub-model. Monte Carlo simulation
operates like a single-shot mode, but in batches of images, and over multiple
channel realizations, so that average performance can be measured.

TABLE I: Main simulator components Component | Options
---|---
Model | Any TensorFlow 2 DNN model
Split layer | Any layer
Quantization | No quantization
Uniform $n$-bit
Packetization |  $r_{p}$ feature rows per packet
Channel model | Independent & identically distributed (iid)
Gilbert-Elliott
Packet traces
Feature loss concealment | None
SiLRTC [10]
HaLRTC [10]
ALTeC [7]
CALTeC [8]
Navier-Stokes [11]
Simulation mode | Single-shot
Monte Carlo

##  IV Experiments

Here we use DFTS2 to conduct what we believe are the most comprehensive
experiments on packet loss resilient collaborative image classification.
Experiments were run on the same 882 images used in [4, 8], which are a subset
of 10 classes of the ImageNet [24] test set. Pre-
trained222https://github.com/qubvel/classification_models DNNs, namely
ResNet-18 [15], ResNet-34 [15], DenseNet-121 [25], and EfficientNet-B0 [26],
were each split at the output of an early layer with a single connection to
downstream layers. This simulation scenario corresponds to a mobile/edge
device having very limited computational capabilities, hence the split at an
early layer. Of course, DFTS2 also supports split layers deeper in the DNN.
Layer add_1 for ResNet-18 and layer add_3 for ResNet-34 correspond to the
output of the shortcut connection and element-wise addition applied to the
final output of the conv2_x layers [15, Table 1]. This means that a single
deep feature tensor $\mathcal{X}\in\mathbb{R}^{56\times 56\times 64}$ is to be
transmitted to the cloud. For DenseNet-121, pool2_conv is the output of
convolution layer in the first transition layer block, as denoted in [25,
Table 1]. In this case, $\mathcal{X}\in\mathbb{R}^{56\times 56\times 128}$. In
the case of EfficientNet-B0 block2b_add tensors,
$\mathcal{X}\in\mathbb{R}^{56\times 56\times 24}$. With $r_{p}=8$ rows of
features/packet, we end up with 7 packets per tensor channel in the case of
ResNet-18 add_1, ResNet-34 add_3, DenseNet-121 pool2_conv, and EfficientNet-B0
block2b_add. When ResNet-18 is split at the output of the add_3 layer, the
transmitted feature tensor has dimensions $28\times 28\times 64$. With
$r_{p}=4$ rows of features/packet, we end up with 7 packets per tensor
channel.

Tensors are min-max quantized to $n=8$ bits before transmission over a
communication channel.

DFTS2 experiments were run using Python 3.6 and TensorFlow 2 in this study.
All packet concealment methods employed are executed on CPU only, though it
may be possible to develop GPU-compatible implementations.

Figure 3: Top-1 cloud sub-model classification accuracy averaged over Monte
Carlo simulations under a range of GE burst lengths on EfficientNet-B0
block2b_add tensors.

(a) Original tensor

(b) Corrupted tensor

(c) CALTeC-repaired tensor

(d) ALTeC-repaired tensor

(h) Meerkat test image

(g) Navier-Stokes inpainting

(f) SiLRTC-50

(e) HaLRTC-50

Figure 4: Tiled ResNet-34 add_3 tensor produced from a (h) Meerkat image: (a)
original; (b) corrupted; repaired with (c) CALTeC, (d) ALTeC, (e) HaLRTC with
50 iterations, (f) SiLRTC with 50 iterations and (g) Navier-Stokes inpainting.
Each tile delineated by green lines represents a channel in a feature tensor.
Images were mapped to grayscale for enhanced visualization.

Packets were transmitted over GE channels with burst loss probabilities
$P_{B}\in\\{0.01,0.10,0.20,0.30\\}$ and average burst lengths of
$L_{B}\in\\{1,2,\dotsc,7\\}$. The range of burst loss probabilities chosen for
Monte Carlo experiments is in accordance with a conventional practice in video
streaming research [16]. All tensor channels in our experiments are packetized
into 7 packets: $56\times 56$ tensor channels with 8 rows of tensor elements
per packet ($7\cdot 8=56$) and $28\times 28$ tensor channels with 4 rows of
features per packet ($7\cdot 4=28$). The range of $L_{B}$ was chosen up to 7,
to cover the case where the entire tensor channel is likely to be lost. In our
Monte Carlo experiments, for each pair of $(P_{B},L_{B})$ values, 20 channels
realizations were simulated for each batch of the 882 image test set. Boolean
loss maps for each test batch over 20 channel realizations of each GE channel
$(P_{B},L_{B})$ were saved. These loss maps were then used to re-create packet
loss patterns so that the same corrupted tensors were presented to each loss
concealment method.

In the results that follow, ‘No Loss’ or ‘NL’ is the performance of the system
without packet loss or quantization, i.e., the performance of the baseline DNN
model without splitting into edge and cloud sub-models. ‘No Completion’or ‘NC’
means no tensor completion, that is, tensor elements from missing packets are
simply set to zero. For each loss probability
$P_{B}\in\\{0.01,0.10,0.20,0.30\\}$, the cloud sub-model Top-1 accuracy
performance was accumulated over all average burst lengths and averaged.
Tables II, III, IV, V, and VI summarize Monte Carlo results for all simulation
scenarios considered. Fig. 3 shows the Top-1 cloud sub-model classification
accuracy over Monte Carlo experiments for EfficientNet-B0 block2b_add tensors.
Each marker in Fig. 3 thus represents 140 Monte Carlo runs through the 882
image test set. For all DNNs, SiLRTC and HaLRTC were run for 50 iterations,
which were found to be sufficient for convergence in [7]. Feature recovery
done with all packet loss concealment methods leads to significant
improvements in Top-1 classification accuracy over the ‘No Completion case in
Fig. 3.

TABLE II: ResNet-18 add_1 tensors: classification accuracy averaged over 140 Monte Carlo runs per burst loss probability $P_{B}$. ‘NL’ refers to the ‘No Loss’ case; ‘NC’ refers to the ‘No Completion’ case; other entries refer to the cloud prediction accuracy on tensors repaired with the indicated method. The highest classification accuracy values in each scenario are indicated in bold.  $\boldsymbol{P_{B}}$ | 0.3 | 0.2 | 0.1 | 0.01
---|---|---|---|---
Top-1 NL | 85.37% | 85.37% | 85.37% | 85.37%
Top-5 NL | 94.44% | 94.44% | 94.44% | 94.44%
Top-1 NC | 34.07% | 55.42% | 74.85% | 84.70%
Top-5 NC | 49.70% | 72.21% | 88.26% | 94.17%
Top-1 ALTeC | 58.57% | 75.75% | 83.05% | 85.31%
Top-5 ALTeC | 74.41% | 88.17% | 93.15% | 94.43%
Top-1 CALTeC | 69.71% | 78.81% | 83.30% | 85.33%
Top-5 CALTeC | 84.54% | 91.04% | 93.64% | 94.43%
Top-1 HaLRTC | 68.77% | 78.82% | 83.51% | 85.33%
Top-5 HaLRTC | 83.54% | 90.94% | 93.70% | 94.43%
Top-1 SiLRTC | 57.43% | 67.21% | 73.55% | 82.66%
Top-5 SiLRTC | 70.50% | 78.74% | 84.00% | 92.47%
Top-1 Navier-Stokes | 58.81% | 74.12% | 81.86% | 85.20%
Top-5 Navier-Stokes | 75.85% | 88.09% | 92.97% | 94.41%
TABLE III: ResNet-18 add_3 tensors: classification accuracy averaged over 140 Monte Carlo runs per burst loss probability $P_{B}$ $\boldsymbol{P_{B}}$ | 0.3 | 0.2 | 0.1 | 0.01
---|---|---|---|---
Top-1 NL | 85.37% | 85.37% | 85.37% | 85.37%
Top-5 NL | 94.44% | 94.44% | 94.44% | 94.44%
Top-1 NC | 65.67% | 75.46% | 81.62% | 85.11%
Top-5 NC | 82.82% | 89.32% | 92.68% | 94.36%
Top-1 ALTeC | 81.62% | 84.01% | 85.16% | 85.31%
Top-5 ALTeC | 92.66% | 93.87% | 94.29% | 94.49%
Top-1 CALTeC | 76.22% | 80.53% | 83.50% | 85.23%
Top-5 CALTeC | 89.47% | 92.00% | 93.63% | 94.45%
Top-1 HaLRTC | 77.00% | 81.75% | 84.31% | 85.28%
Top-5 HaLRTC | 89.84% | 92.71% | 93.94% | 94.47%
Top-1 SiLRTC | 25.97% | 32.66% | 38.62% | 66.33%
Top-5 SiLRTC | 34.32% | 41.31% | 47.93% | 78.07%
Top-1 Navier-Stokes | 77.99% | 81.37% | 83.88% | 85.24%
Top-5 Navier-Stokes | 90.20% | 92.48% | 93.75% | 94.44%
TABLE IV: ResNet-34 add_3 tensors: classification accuracy averaged over 140 Monte Carlo runs per burst loss probability $P_{B}$ $\boldsymbol{P_{B}}$ | 0.3 | 0.2 | 0.1 | 0.01
---|---|---|---|---
Top-1 NL | 85.83% | 85.83% | 85.83% | 85.83%
Top-5 NL | 94.78% | 94.78% | 94.78% | 94.78%
Top-1 NC | 26.27% | 50.94% | 74.93% | 85.58%
Top-5 NC | 41.77% | 68.09% | 87.81% | 94.62%
Top-1 ALTeC | 63.24% | 79.12% | 84.88% | 85.96%
Top-5 ALTeC | 77.98% | 90.20% | 94.04% | 94.82%
Top-1 CALTeC | 77.31% | 82.86% | 85.35% | 85.99%
Top-5 CALTeC | 90.03% | 93.31% | 94.59% | 94.87%
Top-1 HaLRTC | 72.84% | 81.53% | 85.09% | 85.95%
Top-5 HaLRTC | 86.01% | 92.29% | 94.49% | 94.85%
Top-1 SiLRTC | 57.15% | 66.21% | 72.64% | 83.41%
Top-5 SiLRTC | 68.41% | 76.57% | 82.34% | 92.79%
Top-1 Navier-Stokes | 71.45% | 81.08% | 84.97% | 85.98%
Top-5 Navier-Stokes | 85.62% | 92.04% | 94.34% | 94.85%
TABLE V: DenseNet-121 pool2_conv tensors: classification accuracy averaged over 140 Monte Carlo runs per burst loss probability $P_{B}$ $\boldsymbol{P_{B}}$ | 0.3 | 0.2 | 0.1 | 0.01
---|---|---|---|---
Top-1 NL | 89.00% | 89.00% | 89.00% | 89.00%
Top-5 NL | 96.03% | 96.03% | 96.03% | 96.03%
Top-1 NC | 42.15% | 66.84% | 82.50% | 88.52%
Top-5 NC | 56.21% | 80.60% | 92.68% | 95.98%
Top-1 ALTeC | 78.93% | 85.19% | 87.98% | 88.93%
Top-5 ALTeC | 91.15% | 95.08% | 96.02% | 96.06%
Top-1 CALTeC | 81.94% | 85.61% | 87.80% | 88.89%
Top-5 CALTeC | 92.79% | 95.11% | 96.01% | 96.05%
Top-1 HaLRTC | 83.38% | 86.72% | 88.32% | 88.86%
Top-5 HaLRTC | 93.78% | 95.62% | 96.09% | 96.06%
Top-1 SiLRTC | 13.90% | 19.68% | 30.16% | 72.71%
Top-5 SiLRTC | 18.68% | 26.04% | 38.66% | 82.88%
Top-1 Navier-Stokes | 82.09% | 85.87% | 87.85% | 88.83%
Top-5 Navier-Stokes | 92.73% | 95.15% | 95.98% | 96.05%
TABLE VI: EfficientNet-B0 block2b_add tensors: classification accuracy averaged over 140 Monte Carlo runs per burst loss probability $P_{B}$ $\boldsymbol{P_{B}}$ | 0.3 | 0.2 | 0.1 | 0.01
---|---|---|---|---
Top-1 NL | 86.28% | 86.28% | 86.28% | 86.28%
Top-5 NL | 96.15% | 96.15% | 96.15% | 96.15%
Top-1 NC | 0.30% | 2.80% | 17.38% | 75.51%
Top-5 NC | 0.69% | 4.96% | 25.53% | 87.98%
Top-1 ALTeC | 71.64% | 79.64% | 83.89% | 86.05%
Top-5 ALTeC | 85.06% | 91.31% | 94.80% | 96.08%
Top-1 CALTeC | 65.11% | 76.06% | 82.62% | 85.93%
Top-5 CALTeC | 83.55% | 90.89% | 94.65% | 96.10%
Top-1 HaLRTC | 65.95% | 76.69% | 82.96% | 85.98%
Top-5 HaLRTC | 83.40% | 91.11% | 94.82% | 96.09%
Top-1 SiLRTC | 65.71% | 76.51% | 82.83% | 85.97%
Top-5 SiLRTC | 83.07% | 90.88% | 94.71% | 96.08%
Top-1 Navier-Stokes | 61.14% | 74.44% | 82.32% | 85.91%
Top-5 Navier-Stokes | 80.34% | 90.19% | 94.57% | 96.12%

For all DNNs and split layers considered, ALTeC, CALTeC, HaLRTC and Navier-
Stokes inpainting bring clear improvements in cloud Top-1 and Top-5 accuracy
performance at all loss probabilities $P_{B}$ over the ‘No Completion’ case.
Overall, SilRTC has the worst performance. While it shows improvements at
$P_{B}=\\{0.2,0.3\\}$ for ResNet-18 add_1 and ResNet-34 add_3 and at all burst
loss probabilities for EfficientNet-B0 block2b_add, SiLRTC leads to
significantly worse performance than the ‘No Completion/NC’ case in all other
scenarios. With additional iterations, it may be possible to obtain improved
SiLRTC results, at the cost of poorer inference latency.

We present an example of a test image using ResNet-34 split at the output of
add_3. Fig. 4 shows the test image and tiled versions of the front-end sub-
model’s original output tensor, the same tensor corrupted after transmission
over a GE channel with $(P_{B},L_{B})=(0.3,7)$ and after repairs by various
methods in DFTS2. When no packet loss concealment was done on the corrupted
tensor, the back-end sub-model gave an incorrect Top-1 predicted class with
$29\%$ confidence. Meanwhile, all five packet concealment methods produced
repaired tensors which led to the correct class (‘meerkat’) with high
confidence ($>93\%$ with all five methods). Repaired tensors in Fig. 4 are
visually close to the original tensor.

##  V Conclusions

We presented a simulation framework called DFTS2 to study deep feature
transmission in collaborative intelligence. We discussed the main components
of this simulator, namely, the packetization, communication channel models,
and the packet loss concealment schemes. We also presented a comprehensive
study of packet loss concealment methods for collaborative image
classification. Among these, HaLRTC and CALTeC produced the best Top-1
classification accuracy performance overall in Monte Carlo experiments on
three deep neural networks. The proposed simulation framework is flexible and
extensible. It can be readily adapted to investigate deep feature compression
in collaborative intelligence, as well as the trade-off between compression
and error resilience. It can also be easily extended to study other
collaborative intelligence and video coding for machines (VCM) scenarios, such
as collaborative object detection, segmentation, tracking, and so on.

## References

  * [1] Y. Kang, J. Hauswald, C. Gao, A. Rovinski, T. N. Mudge, J. Mars, and L. Tang,  “Neurosurgeon: Collaborative intelligence between the cloud and mobile edge,”  in Proc. ASPLOS, Apr. 2017.
  * [2] I. V. Bajić, W. Lin, and Y. Tian,  “Collaborative intelligence: Challenges and opportunities,”  in Proc. IEEE ICASSP, 2021, pp. 8493–8497.
  * [3] A. E. Eshratifar, M. S. Abrishami, and M. Pedram,  “JointDNN: An efficient training and inference engine for intelligent mobile cloud computing services,”  IEEE Trans. Mobile Computing, vol. 20, no. 2, pp. 565–576, Feb. 2021.
  * [4] H. Unnibhavi, H. Choi, S. R. Alvar, and I. V. Bajić,  “DFTS: Deep feature transmission simulator,” IEEE MMSP’18 demo, 2018\.
  * [5] K. Choi, K. Tatwawadi, A. Grover, T. Weissman, and S. Ermon,  “Neural joint source-channel coding,”  in Proc. ICML, Jun. 2019, pp. 1182–1192.
  * [6] J. Shao and J. Zhang,  “Bottlenet++: An end-to-end approach for feature compression in device-edge co-inference systems,”  in Proc. IEEE ICC Workshops, Jun. 2020, pp. 1–6.
  * [7] L. Bragilevsky and I. V. Bajić,  “Tensor completion methods for collaborative intelligence,”  IEEE Access, vol. 8, pp. 41162–41174, 2020.
  * [8] A. Dhondea, R. A. Cohen, and I. V. Bajić,  “CALTeC: Content-adaptive linear tensor completion for collaborative intelligence,”  in Proc. IEEE ICIP, Anchorage, AK, Sep. 2021, pp. 2179–2183.
  * [9] I. V. Bajić,  “Latent space inpainting for loss-resilient collaborative object detection,”  in Proc. IEEE ICC, Jun. 2021, pp. 1–6.
  * [10] J. Liu, P. Musialski, P. Wonka, and J. Ye,  “Tensor completion for estimating missing values in visual data,”  IEEE Trans. Pattern Analysis and Machine Intelligence, vol. 35, no. 1, pp. 208–220, Jan. 2013.
  * [11] M. Bertalmio, A. L. Bertozzi, and G. Sapiro,  “Navier-Stokes, fluid dynamics, and image and video inpainting,”  in Proc. IEEE/CVF CVPR, 2001, vol. 1, pp. I–355–I–362.
  * [12] A. Telea,  “An image inpainting technique based on the fast marching method,”  J. Graphics Tools, vol. 9, no. 1, pp. 23–34, 2004.
  * [13] G. Hasslinger and O. Hohlfeld,  “The Gilbert-Elliott model for packet loss in real time services on the Internet,”  in Proc. 14th GI/ITG Conference - Measurement, Modelling and Evaluation of Computer and Communication Systems, 2008, pp. 1–15.
  * [14] R. A. Cohen, H. Choi, and I. V. Bajić,  “Lightweight compression of intermediate neural network features for collaborative intelligence,”  IEEE Open Journal of Circuits and Systems, vol. 2, pp. 350–362, 2021.
  * [15] K. He, X. Zhang, S. Ren, and J. Sun,  “Deep residual learning for image recognition,”  in Proc. IEEE/CVF CVPR, Jun. 2016, pp. 770–778.
  * [16] Y. Wang, J. Ostermann, and Y.-Q. Zhang,  Video Processing and Communications,  Prentice Hall, 2002.
  * [17] H. Choi and I. V. Bajić,  “Deep feature compression for collaborative object detection,”  in Proc. IEEE ICIP, Oct. 2018, pp. 3743–3747.
  * [18] H. Choi and I. V. Bajić,  “Near-lossless deep feature compression for collaborative intelligence,”  in Proc. IEEE MMSP, 2018, pp. 1–6.
  * [19] S. R. Alvar and I. V. Bajić,  “Multi-task learning with compressible features for collaborative intelligence,”  in Proc. IEEE ICIP, 2019, pp. 1705–1709.
  * [20] H. Choi, R. A. Cohen, and I. V. Bajić,  “Back-and-forth prediction for deep tensor compression,”  in Proc. IEEE ICASSP, 2020, pp. 4467–4471.
  * [21] R. A. Cohen, H. Choi, and I. V. Bajić,  “Lightweight compression of neural network feature tensors for collaborative intelligence,”  in Proc. IEEE ICME, Jul. 2020, pp. 1–6.
  * [22] S. R. Alvar and I. V. Bajić,  “Pareto-optimal bit allocation for collaborative intelligence,”  IEEE Transactions on Image Processing, vol. 30, pp. 3348–3361, 2021\.
  * [23] K. Stuhlmuller, N. Farber, M. Link, and B. Girod,  “Analysis of video transmission over lossy channels,”  IEEE Journal on Selected Areas in Communications, vol. 18, no. 6, pp. 1012–1032, 2000.
  * [24] O. Russakovsky et al.,  “Imagenet large scale visual recognition challenge,”  Int. J. Comput. Vision, vol. 115, no. 3, pp. 211–252, Dec. 2015\.
  * [25] G. Huang, Z. Liu, L. Van Der Maaten, and K. Q. Weinberger,  “Densely connected convolutional networks,”  in Proc. IEEE/CVF CVPR, 2017, pp. 2261–2269.
  * [26] M. Tan and Q. V. Le,  “EfficientNet: rethinking model scaling for convolutional neural networks,”  in Proc. ICML, Jun. 2019, pp. 6105–6114.


* [163] D. Dudal and T. G. Mertens, “Melting of charmonium in a magnetic field from an effective AdS/QCD model,” Phys. Rev. D 91, 086002 (2015) [arXiv:1410.3297 [hep-th]].
  * [164] D. Dudal, A. Hajilou and S. Mahapatra, “A quenched 2-flavour Einstein–Maxwell–Dilaton gauge-gravity model,” Eur. Phys. J. A 57, no.4, 142 (2021) [arXiv:2103.01185 [hep-th]].


11institutetext: Department of Astrophysics/IMAPP, Radboud University, P.O.
Box 9010, NL-6500 GL Nijmegen, The Netherlands
22institutetext: Department of Mathematics/IMAPP, Radboud University, P.O. Box
9010, NL-6500 GL Nijmegen, The Netherlands
33institutetext: SRON, Netherlands Institute for Space Research, Sorbonnelaan
2, NL-3584 CA Utrecht, The Netherlands 44institutetext: Institute of
Astronomy, KU Leuven, Celestijnenlaan 200D, B-3001 Leuven, Belgium
55institutetext: Key Laboratory of Particle Astrophysics, Institute of High
Energy Physics, Chinese Academy of Sciences, 19B Yuquan Road, Beijing 100049,
People’s Republic of China

# The GW-Universe Toolbox II: constraining the binary black hole population
with second and third generation detectors

Shu-Xu Yi<EMAIL_ADDRESS>Fiorenzo Stoppa 1122 Gijs Nelemans 113344 Eric
Cator 22

Context:  The GW-Universe Toolbox is a software package that simulates
observations of the Gravitational Wave (GW) Universe with different types of
GW detectors, including Earth-based/Space-borne laser interferometers and
Pulsar timing arrays. It is accessible as a website, and can also be imported
and run locally as a Python package.
Methods: We employ the method used by the GW-Universe Toolbox to generate a
synthetic catalogue of detection of stellar mass binary black hole (BBH)
mergers. As an example of its scientific application, we study how GW
observations of BBHs can be used to constrain the merger rate as function of
redshift and masses. We study advanced LIGO (aLIGO) and Einstein Telescope
(ET) as _two representatives for the 2nd and 3rd generation GW observatories_.
We also simulate the observations from a detector that is half as sensitive as
the ET at design which represents an early phase of ET. Two methods are used
to obtain the constraints on the source population properties from the
catalogues: the first uses a parameteric differential merger rate model and
applies a Bayesian inference on the parameters; The other one is non-
parameteric and uses weighted Kernel density estimators.
Results: The results show the overwhelming advantages of the 3rd generation
detector over the 2nd generation for the study of BBH population properties,
especially at a redshifts higher than $\sim 2$, where the merger rate is
believed to peak. With the simulated aLIGO catalogue, the parameteric Bayesian
method can still give some constraints on the merger rate density and mass
function beyond its detecting horizon, while the non-parametric method lose
the constraining ability completely there. The difference is due to the extra
information placed by assuming a specific parameterisation of the population
model in the Bayesian method. In the non-parameteric method, no assumption of
the general shape of the merger rate density and mass function are placed, not
even the assumption of its smoothness. These two methods represent the two
extreme situations of general population reconstruction. We also find that,
despite the numbers of detection of the half-ET can be easily compatible with
full ET after a longer observation duration, the catalogue from the full ET
can still give much better constraints on the population properties, due to
its smaller uncertainties on the physical parameters of the GW events.

###### Key Words.:

gravitational waves; method:statistical; stars: black holes

##  1 Introduction

The era of gravitational wave (GW) astronomy has begun with the first
detection of a binary black hole (BBH) merger in 2015 (Abbott et al., 2016).
There have been $\sim 90$ BBH systems discovered with $\sim 3$ years
observation in the LIGO/Virgo O1/O2/O3 runs (GWTC-1&2 catalogues: Abbott et
al. (2019b, 2020)). The number is catching up on that of the discovered
Galactic black holes (BH) in X-ray binaries during the last half century
(Corral-Santana et al., 2016). LIGO and Virgo are being upgraded towards their
designed sensitivity, and new GW observatories such as The Kamioka
Gravitational Wave Detector (KAGRA, Kagra Collaboration et al. 2019) and LIGO-
India (Iyer et al., 2011) will join the network in the near future.
Furthermore, there are plans for the next generation of GW observatories such
as the Einstein Telescope (ET, Punturo et al. 2010) and the Cosmic Explore
(CE, (Reitze et al., 2019)), which will push forward the detecting horizon
significantly. With the expectation of the exploding size of the BBH
population, the community has been actively studying what science can be
extracted from GW observations of BBHs. Applications include e.g. the cosmic
merger rate density of BBH (e.g. Vitale et al., 2019), the mass distribution
of stellar mass black hole (e.g. Kovetz, 2017; Kovetz et al., 2017; Abbott et
al., 2019a, 2021a), the rate of stellar nucleosynthesis (e.g. Farmer et al.,
2020) and constraints on the expansion history of the universe (Farr et al.,
2019). Simulating a catalogue of observations is the crucial step in the above
mentioned work. In all the previous literature, the simulation is performed
for specific detectors and observation duration and their catalogues are
difficult to extend to other usages. Therefore, we are building a set of
simulating tools which has the flexibility to different detectors, observation
duration, source populations and cosmological models. The GW-Universe Toolbox
is such a tool set that simulates observations of the GW Universe (including
populations of stellar mass compact binary mergers, supermassive black hole
binaries inspiral and mergers, extreme mass ratio inspirals, etc.) with
different types of GW detectors including Earth-based/Space-borne laser
interferometers and Pulsar timing arrays. It is accessible as a website111gw-
universe.org, or run locally as a Python package. For a detailed overview of
the functionalities and methods of GW-Universe Toolbox, please see Yi et al.
2021, referred to as the Paper-I hereafter.

Here we exhibit one application of the GW-Universe Toolbox by showing how the
cosmic merger rate and mass function of BH can be constrained by GW
observations using advanced LIGO (aLIGO, at design sensitivity) and ET 222The
current GW-Universe Toolbox can only simulate observations with single
detectors. A network of detectors is expected to be more efficient in the
following aspects: 1. the total SNR of a joint-observed GW source would be
larger (but in the same order of magnitude) than that from a single detector;
2. the uncertainties of source parameters would be smaller. In the current GW-
Universe Toolbox, the duty cycle of detectors are assumed to be 100%. In a
more realized simulation, where the duty cycles of single detectors are not
100%, a detector network also benefits from a larger time coverage. . The
cosmic merger rate density, or the volumetric merger rate as function of
redshift $\mathcal{R}(z)$, provides valuable information on the history of the
star formation rate (SFR), which is crucial for understanding the evolution of
the Universe, e.g, the chemical enrichment history, and the formation and
evolution of galaxies (Madau & Dickinson, 2014). The SFR as function of
redshift has been studied for decades mainly with UV/optical/infrared
observation of galaxies in different redshift (Serjeant et al., 2002; Doherty
et al., 2006; Thompson et al., 2006; Carilli et al., 2008; Magdis et al.,
2010; Rujopakarn et al., 2010; Ly et al., 2011; Sedgwick et al., 2011;
Kennicutt & Evans, 2012; Curtis-Lake et al., 2013; Madau & Dickinson, 2014).
Such electromagnetic (EM) wave probes are limited by dust extinction and the
completeness of samples at high redshifts is difficult to estimate (e.g.
Maiolino & Mannucci, 2019). Moreover, the connection between the direct
observable and the SFR history is model dependent on multiple layers (see e.g.
Chruślińska et al., 2021). GW provide a unique new probe which suffer less
from the above mentioned limitations.

Moreover, $\mathcal{R}(z)$ depends also on the distribution of delay times
$\mathbf{\tau}$, i.e., the time lag between the formation of the binary
massive stars and the merger of the two compact objects (e.g. Belczynski et
al., 2010; Kruckow et al., 2018). Therefore, a combination of the information
of the SFR from EM observations and $\mathcal{R}(z)$ from GW can improve our
knowledge of the delay time distribution. $P(\tau)$ together with primary mass
function $p(m_{1})$ (where $m_{1}$ is the mass of the primary BH) of BBH helps
us to better understand the physics in massive binary evolution, e.g., in the
common envelope scenario (e.g. Belczynski et al., 2016; Kruckow et al., 2018).

There are several previous works that studied the prospect of using simulated
GW observations to study $\mathcal{R}(z)$ and $p(m_{1})$ of BBH. Vitale et al.
(2019) shows that with three months observation with a third-generation
detector network, $\mathcal{R}(z)$ can be constrained to high precision up to
$z\sim 10$. In their work the authors did not take into account the signal-to-
noise (SNR) limits of detectors and instead made the assumption that all BBH
mergers in the Universe are in the detected catalogue. Kovetz et al. (2017)
studied how the observation with aLIGO can be used to give a constraint on
$p(m_{1})$. They assumed a certain parameterization of $p(m_{1})$, and use an
approximated method to obtain the covariance of the parameters based on a
simulated catalogue of aLIGO observations. Recently, Singh et al. (2021)
simulated observation with ET on compact object binary mergers, taking into
account the effect of the rotation of the Earth. They found that
$\mathcal{R}(z)$ and mass distribution can be so well constrained from ET’s
observation, that mergers from different populations can be distinguished.

In this work, we simulate the observed catalogues of BBH mergers in ten years
of aLIGO observation with its design sensitivity and one month ET observation
using the GW-Universe Toolbox. We also simulate the observation from a
detector that is half as sensitive as the ET in design which represents the
early phase of ET. The paper is organised as follows: In section 2, we
summarize the method used by the GW-Universe Toolbox to simulate the catalogue
of BBH mergers; In section 3, we use the synthetic catalogues to set
constraints on the cosmic merger rate and mass function of the population. We
employ two ways to constrain $\mathcal{R}(z)$ along with the mass function of
the primary BH. The first assumes parameteric forms of $\mathcal{R}(z)$ and
the mass function and uses Bayesian inference on their parameters. The other
method does not assume a parameteric formula and uses weighted Kernel density
estimators. The results are presented and compared. We conclude and discuss in
section 4.

##  2 Generating synthetic catalogues

###  2.1 method

We use the GW-Universe Toolbox to simulate observations on BBH with three
different detectors, namely, aLIGO (designed sensitivity), ET and “half-ET”
333Here we are referring to the ET-D design. For the key parameters of this ET
configuration, see Hild et al. (2011).. The noise power spectrum $S_{\rm{n}}$
for aLIGO and ET are obtained from
444https://dcc.ligo.org/LIGO-T1800044/public
http://www.et-gw.eu/index.php/etsensitivities. The “half-ET” represents our
expectation on the performance of the ET in its early phase. In figure 1, we
plot the $\sqrt{S_{\rm{n}}}$ as function of frequency for aLIGO, ET at their
designed sensitivity, along with that of the half-ET. For its $S_{\rm{n}}$, we
multiply $S_{\rm{n}}$ of the ET with a factor 10.89 to obtain a noise curve
with is in the half way between the designed aLIGO and ET.

Figure 1: The squared-root of the noise power spectrum density of the aLIGO
(blue) and ET (orange) at their designed sensitivity. The data of the plots
are obtained from https://dcc.ligo.org/LIGO-T1800044/public and http://www.et-
gw.eu/index.php/etsensitivities

Here we will give a brief summary of the method. For details, we refer to Yi
et al. (2021). We first calculate a function $\mathcal{D}(z,m_{1},m_{2})$,
which represents the detectability of a BBH with redshift $z$, primary and
secondary masses $m_{1}$, $m_{2}$. $\mathcal{D}(z,m_{1},m_{2})$ depends on the
antenna pattern and $S_{\rm{n}}$ of the detector, and the user designated
signal-to-noise (SNR) threshold. The other ingredient is the merger rate
density distribution of BBH in the Universe, $\dot{n}(z,m_{1},m_{2})$. In the
GW-Universe Toolbox, we integrate two parameterized function forms for
$\dot{n}(z,m_{1},m_{2})$, namely pop-A and pop-B (see Yi et al. 2021 for the
description on the population models). For more sophisticated
$\dot{n}(z,m_{1},m_{2})$ from population synthesis simulations, see Belczynski
et al. (2002); Mapelli et al. (2017); Mapelli & Giacobbo (2018); Arca Sedda et
al. (2021); Banerjee (2021); Mapelli et al. (2021); van Son et al. (2021);
Singh et al. (2021).

The distribution of detectable sources, which we denote as $N(z,m_{1},m_{2})$
can be calculated with:

$N_{\rm{D}}(\mathbf{\Theta})=\frac{T}{1+z}\frac{dV_{\rm{c}}}{dz}\dot{n}(\mathbf{\Theta})\mathcal{D}(\mathbf{\Theta}),$
(1)

where $T$ is the observation duration, and $\mathbf{\Theta}=(z,m_{1},m_{2})$.
An integration of $N_{\rm{D}}(\mathbf{\Theta})$ gives the expectation of total
detection number:

$N_{\rm{tot,exp}}=\int d\mathbf{\Theta}N_{D}(\mathbf{\Theta}).$ (2)

The actual total number in a realization of simulated catalogue is a Poisson
random around this expectation. The synthetic catalogue of GW events is drawn
from $N_{D}(\mathbf{\Theta})$ with a Markov-Chain Monte-Carlo (MCMC)
algorithm.

We can also give a rough estimation of the uncertainties on the intrinsic
parameters of each event in the synthetic catalogue using the Fisher
Information Matrix (FIM). Previous studies (e.g., Rodriguez et al. (2013))
found that the Fisher information matrix will sometimes severely overestimated
the uncertainties, especially when the SNR is close to the threshold. We
implement a correction that if the relative uncertainty of masses calculated
with FIM is larger than $20\%$, $\delta m_{1,2}=0.2m_{1,2}$ is applied
instead.

The uncertainty of the redshift is in correlation with those of other external
parameters, such as the localization error. The actual localization
uncertainties are largely determined by the triangulation error of the
detector network, which can not be estimated with the current version of GW-
Universe Toolbox yet. As a result, we do not have a reliable way that
accurately estimates the errors on redshift. As a ad hoc method, we assign
$\delta z=0.5z$ as a typical representative for the uncertainties of aLIGO,
and $\delta z=0.017z+0.012$ as a typical representative for the uncertainties
of ET (Vitale & Evans, 2017). For half-ET, we assign $\delta z=0.1z+0.01$.

###  2.2 Catalogues

Using the above-mentioned method, we generate synthetic catalogues
corresponding to 10-years aLIGO, 1-month ET observations and 1-month
observation with half-ET. In figures 2, 3 and 4, we plot events in those
catalogues, along with the corresponding marginalized
$N_{D}(\mathbf{\Theta})$. The numbers in the catalogues are 2072, 1830 and 889
for the 10-year aLIGO, 1-month ET and 1-month half-ET respectively. The
catalogues agree with the theoretical number density, which proves the
validity of the MCMC sampling process.

Figure 2: aLIGO 10 years’ catalogue:  Black dashed lines are contours of
$N_{D}(\mathbf{\Theta})$; Blue points with error bars are 2072 events in the
simulated catalogues. Figure 3: ET 1 month’s catalogue:  Black dashed lines
are contours of $N_{D}(\mathbf{\Theta})$; Blue points with error bars are 1830
events in the simulated catalogues. Figure 4: half-ET 1 month’s catalogue:
Black dashed lines are contours of $N_{D}(\mathbf{\Theta})$; Blue points with
error bars are 889 events in the simulated catalogues.

The underlying population model for $\dot{n}(z,m_{1},m_{2})$ is pop-A, with
parameters: $R_{\rm{n}}=13$ Gpc-3yr-1, $\tau=3$ Gyr, $\mu=3$, $c=6$,
$\gamma=2.5$, $m_{\rm cut}=95\,M_{\odot}$, $q_{\rm cut}=0.4$, $\sigma=0.1$.
Those parameters are chosen because they result in simulated catalogues which
are in compatible with the real observations (see Yi et al. 2021).

##  3 Constraints on $\mathcal{R}(z)$ and $p(m_{1})$ from BBH catalogues

In this section, we study how well $\mathcal{R}$ and $p(m_{1})$ can be
constrained with BBH catalogues from observations with 2nd and 3rd generations
GW detectors. Two distinct methods will be used.

###  3.1 Bayesian method

In the first method, we use a parameterised form of
$\dot{n}(\mathbf{\Theta}|\mathcal{B})$. $\mathcal{B}$ are the parameters
describing the population rate density, which are
$\mathcal{B}=(R_{0},\tau,\mu,c,\gamma,q_{\rm cut},\sigma)$ in our case.
$\mathbf{\Theta}^{k}$ are the physical parameters of the $k$-th source in the
catalogue, and we use $\\{\mathbf{\Theta}\\}$ to denote the whole catalogue.

The posterior probability of $\mathcal{B}$ given an observed catalogue
$\\{\mathbf{\Theta}\\}_{\rm obs}$:

$p(\mathcal{B}|\\{\mathbf{\Theta}\\}_{\rm obs})\propto
p(\\{\mathbf{\Theta}\\}_{\rm obs}|\mathcal{B})p(\mathcal{B}),$ (3)

where $p(\\{\mathbf{\Theta}\\}_{\rm obs}|\mathcal{B})$ is the likelihood and
$p(\mathcal{B})$ is the prior. The subscript “obs” is to distinguish the
observed parameters from the true parameters of the sources. The likelihood of
detecting a catalogue $\\{\mathbf{\Theta}\\}_{\rm obs}$ given the parameters
$\mathcal{B}$ is a known result as an inhomogeneous Poisson Process (Foreman-
Mackey et al., 2014):

$p(\\{\mathbf{\Theta}\\}_{\rm
obs}|\mathcal{B})=N!\prod_{k=1}^{N}p(\mathbf{\Theta}^{k}_{\rm
obs}|\mathcal{B})p(N|\lambda(\mathcal{B})),$ (4)

where $p(N|\lambda(\mathcal{B}))$ is the probability of detecting $N$ sources
if the expected total number is $\lambda$, which is a Poisson function. The
probability of detecting a source $k$ with parameters $\mathbf{\Theta}^{k}$
is:

$p(\mathbf{\Theta}^{k}_{\rm obs}|\mathcal{B})=\int p(\mathbf{\Theta}^{k}_{\rm
obs}|\mathbf{\Theta}^{k}_{\rm true})p(\mathbf{\Theta}^{k}_{\rm
true}|\mathcal{B})d\mathbf{\Theta}^{k}_{\rm true}.$ (5)

We further assume that observational errors are Gaussian, therefore we have
$p(\mathbf{\Theta}^{k}_{\rm obs}|\mathbf{\Theta}^{k}_{\rm
true})=p(\mathbf{\Theta}^{k}_{\rm true}|\mathbf{\Theta}^{k}_{\rm obs})$.
Taking this symmetric relation into equation 5, we have:

$p(\mathbf{\Theta}^{k}_{\rm obs}|\mathcal{B})=\int p(\mathbf{\Theta}^{k}_{\rm
true}|\mathbf{\Theta}^{k}_{\rm obs})p(\mathbf{\Theta}^{k}_{\rm
true}|\mathcal{B})d\mathbf{\Theta}^{k}_{\rm true}.$ (6)

The above integration is equivalent to the average of
$p(\mathbf{\Theta}^{k}_{\rm true}|\mathcal{B})$ over all possible
$\Theta^{k}$. Therefore:

$p(\mathbf{\Theta}_{\rm
obs}|\mathcal{B})\approx\left<p(\mathbf{\Theta}^{k}_{\rm
true}|\mathcal{B})\right>,$ (7)

where the $\left<\cdots\right>$ denotes an average among a sample of
$\mathbf{\Theta}^{k}$ which is drawn from a multivariate Gaussian
characterized by the observation uncertainties. Taking equation (7) into
equation (4), and noting that $p(\mathbf{\Theta}^{k}|\mathcal{B})\equiv
N_{D}(\mathbf{\Theta}^{k}|\mathcal{B})/\lambda(\mathcal{B})$, and
$p(N|\lambda(\mathcal{B}))$ is the Poisson distribution:

$p(\\{\mathbf{\Theta}\\}_{\rm
obs}|\mathcal{B})=\prod_{k=1}^{N}\left<N_{D}(\mathbf{\Theta}^{k}|\mathcal{B})\right>\exp(-\lambda(\mathcal{B})).$
(8)

For the prior distribution $p(\mathcal{B})$, we use the log normal
distributions for $\tau$, $R_{0}$, $c$, $\gamma$ and $\mu$, which center at
the true values and large enough standard deviations. The posterior
probability distributions of the parameters are obtained with MCMC sampling
from the posterior in equation (3). We keep the nuisance parameters
$q_{\rm{cut}}$ and $\sigma$ fixed to the true value, in order to reduce the
necessary length of chain. The posterior of the interesting parameters would
not be influenced significantly by leaving these free.

Figure 5: The MCMC chains sampled from the posterior distribution of the
population parameters. The blue line is for aLIGO-10 years observation; The
green line is for half-ET’s 1 month observation, and the orange line is for
full ET’s 1 month observation. The black dashed lines mark the true values of
the corresponding parameters, which were used to generate the observed
catalogue.

In Figure 5, we plot the MCMC chains sampled from the posterior distribution
of the population parameters, where we can see the chains are well mixed, and
thus can serve as good representation of the underlying posterior
distribution. In Figure 6, we plot the correlation between the parameters of
the red-shift dependence ($R_{\rm{n}}$, $\tau$) and among the parameters of
the mass function ($\mu$, $c$, $\gamma$). As can be found, there are strong
correlations between the $\mathcal{R}(z)$ parameters $\tau$ and $R_{\rm{n}}$,
and also among the $p(m_{1})$ parameters $\gamma$, $c$ and $\mu$. The
correlation between parameters of the red-shift dependence and those of the
mass function are small and therefore not plotted.

Figure 6: The correlation between the parameters of red-shift dependence
($R_{\rm{n}}$, $\tau$, left panel) and among the parameters of mass function
($\mu$, $c$, $\gamma$, middle and right panel). Blue markers correspond to 10
years observation with aLIGO; Green markers correspond to 1 month observation
with half-ET, and orange markers correspond to 1 month observation with ET.
The red cross in each panel marks the location of the true values, which were
used to generate the observed catalogue.

In figure 7, we plot the inferred differential merger rate as function of
redshift and mass of the primary BH. The bands correspond to the 95%
confidence level of the posterior distribution of the population parameters.
We can see from the panels of figure 7 that, although ten years of aLIGO
observation and 1-month of ET observation result in a similar number of
detections, ET can still obtain much better constraints on both
$\mathcal{R}(z)$ and $p(M_{1})$. That results from 1): the much smaller
uncertainties on parameters with ET observation compared to aLIGO, and 2): The
catalogue observed by ET covers larger redshift and mass ranges than that of
aLIGO. For half-ET and ET, at low red-shift and high mass, both detectors
result in a similar number of detections. Therefore, they give comparable
constraints on $\mathcal{R}(z)$ at low redshift and $p(m_{1})$ at large mass.
While at high red-shift and low mass, the full ET observed more events and
with lower uncertainties, therefore gives better constrains than that of half-
ET.

We see from figure 2 that the maximum red-shift in aLIGO’s catalogue is $z\sim
1.5$, while $\mathcal{R}(z)$ can still be constrained for $z>1.5$, and even be
more stringent towards larger red-shift. This constraint on $\mathcal{R}(z)$
at large redshift is not from the observed catalogue, but is intrinsically
imposed by the assumed model in the Bayesian method: we have fixed our star
formation rate as function of redshift so that $\psi(z_{f})$ peaks at $z\sim
2$, the peak of BBH merger rate can therefore only appear later.

Figure 7: $\mathcal{R}(z)$ (upper panel) and $p(m_{1})$ (lower panel)
reconstructed with the Bayesian method, from the synthetic aLIGO (blue), ET
(orange) and half-ET (green) catalogues. The black dashed curves are the model
used to generate the catalogues.

###  3.2 Nonparametric method

From equation (1), we know that:

$\dot{n}(\mathbf{\Theta})=\frac{N_{\rm{D}}(\mathbf{\Theta})}{\mathcal{D(\mathbf{\Theta})}}\frac{1+z}{T}/\frac{dV_{c}}{dz},$
(9)

and $N_{\rm{D}}(\mathbf{\Theta})$ can be inferred from the observed catalogue.

Marginalising over $m_{1}$ and $m_{2}$ in equation (9) on both sides, we have:

$\mathcal{R}(z)=I(z)|_{w_{i}=1/\mathcal{D}_{i}}\frac{1+z}{T}/\frac{dV_{c}}{dz},$
(10)

where $I(z)$ is the 1-D intensity function (number density) over $z$, of the
observed catalogue, with weights $w_{i}$ equal to the corresponding
$1/\mathcal{D}(\Theta_{i})$. The intensity can be calculated with an empirical
density estimator of the (synthetic) catalogue or with the (non-normalized)
weighted Kernel Density Estimator (KDE, see appendix A). By simulating a
number of different realisations of the Universe population inferred from
$I(z)$, the confidence intervals of the estimated $I(z)$ are obtained.

Here, in order to include the observational uncertainties and give at the same
time an estimate of the Confidence Intervals of the KDE, a specific bootstrap
is employed, as follows:

  1. 1.

The GW-Universe Toolbox returns a catalogue
$\\{\mathbf{\Theta}\\}_{\rm{true}}$, where $\mathbf{\Theta}^{k}_{\rm{true}}$
denotes the true parameters of the $k$-th events in the catalogue, including
its $m_{1}$, $m_{2}$, $z$ and their uncertainties $\delta m_{1}$, $\delta
m_{2}$, $\delta z$. $k$ goes from 1 to $N$, where $N$ is the total number of
events.

  2. 2.

Shift the true values $\mathbf{\Theta}^{k}_{\rm{true}}$ according to a
multivariate Gaussian centred on the true values and with covariance matrix
given by the corresponding uncertainties. Therefore, a new catalogue
$\\{\mathcal{\mathbf{\Theta}}^{k}\\}_{\rm{obs}}$ is obtained.

  3. 3.

Estimate the detection probabilities
$\mathcal{D}(\mathbf{\Theta}^{k}_{\rm{obs}})$.

  4. 4.

Calculate the (unnormalized) weighted KDE of
$\\{\mathcal{\mathbf{\Theta}}^{k}_{\rm{obs}}\\}$, with weights
$1/\mathcal{D}(\mathbf{\Theta}^{k}_{\rm{obs}})$. The intensity is estimated on
the log and then transformed back to the original support; this ensures that
there isn’t any leakage of the intensity on negative support.

  5. 5.

Take $M\equiv\sum_{i=k}^{N}1/\mathcal{D}(\mathbf{\Theta}^{k}_{\rm{obs}})$, as
an estimate of the total number of mergers in the Universe in the parameters
space.

  6. 6.

Generate a realization $\widetilde{M}$ according to a Poisson with expected
value $M$.

  7. 7.

Draw a sample of size $\widetilde{M}$ from the KDE calculated in step 4. This
is the equivalent of generating a new Universe sample, based on the estimate
of the distribution in the Universe from the observations.

  8. 8.

For each of the events in the above sample, calculate its detection
probability. Draw $\widetilde{M}$ uniform random values, $U$, in between 0 and
1, and drop the observations for whose detectability is less than $U$. A
smaller sample is left, which represents a new realization of the detected
catalogue.

  9. 9.

Estimate the covariance for each observation of the generated data set
$\mathbf{\Sigma_{{\Theta}^{k}}}$ using interpolation of the original
catalogue, and shift them according to a multivariate Gaussian centred on each
observation and with covariance matrix $\mathbf{\Sigma_{{\Theta}^{k}}}$.

  10. 10.

Estimate their detection probability.

  11. 11.

Calculate the (unnormalized) weighted KDE of the above mentioned sample.

  12. 12.

Repeat steps from 6 to 11 for a large number of times $B$ (we find $B$=200 is
large enough). Calculate the confidence intervals according to the assigned
quantiles.

  13. 13.

Based on the results above, we can correct for the bias introduced by the
measurement uncertainties. We estimate the bias between the bootstrap
intensities and the intensity estimated on the observed sample, and use it to
rescale the observed sample’s estimated intensity.

Similarly, the merger rate in the Universe as function of the mass of primary
BH is:

$\mathcal{R}(m_{1})=p(m_{1})\times\dot{n}_{\rm{tot}}=I(m_{1})|_{w_{i}=1/\mathcal{D}_{i}}/T,$
(11)

where $p(m_{1})$ is the normalized mass function of the primary BH,
$\dot{n}_{\rm{tot}}$ is the total number of BBH mergers in the Universe in a
unit local time. $I(m_{1})$ is the 1-D intensity function of $m_{1}$,
calculated from the observed catalogue. The procedures of obtaining $I(m_{1})$
is the same as that in the calculation of $I(z)$.

In figures 8 and 10 we plot the $\mathcal{R}(z)$ and $\mathcal{R}(m_{1})$
reconstructed with equations 10 and 11 from the synthetic aLIGO, ET and half-
ET catalogues. The shaded areas are the $95\%$ confidence intervals.

We note in figure 8 that, for aLIGO’s observation, the obtained
$\mathcal{R}(z)$ in $z>0.3$ is an underestimate of the underlying model. The
reason is that aLIGO’s detectable range of BBH masses is increasingly
incomplete towards higher redshift. An intuitive illustration is in Figure 9.
There, we plot the two dimensional merger rate density as function of $m_{1}$
and $z$. The underlying black region marks the parameters space region where
the detectability is essentially zero, i.e., the region beyond aLIGO’s
horizon. When plotting Figure 9, we set the cut at
$\mathcal{D}(z,m_{1},m_{2}=m_{1})<10^{-2}$. We can see the fraction of the
black region in the mass spectrum is increasing towards high redshift. This
corresponds to the fact that BBH with small masses cannot be detected with
aLIGO at higher redshift.

Similarly, in figure 10, we see $\mathcal{R}(m_{1})$ is deviating the theory
towards low $m_{1}$. To compare the estimation with the model, we limit both
the aLIGO catalogue and the theory in the range $z<0.6$. Therefore, the
meaning of $\mathcal{R}(m_{1})$ is no longer the merger rate density of the
complete Universe, but that of a local Universe up to 0.6, which is calculated
with:

$\mathcal{R}(m_{1})|_{z<0.6}=p(m_{1})\int^{0.6}_{0}\mathcal{R}(z)\frac{dV}{dz}/(1+z)dz,$
(12)

where $p(m_{1})$ is the normalized mass function.

It is unlike the results with the Bayesian inference method, where constraints
can still be placed in regimes with no detection. The constraints on
$\mathcal{R}(z)$ and the mass function are also less tight for the
nonparametric method comparing with the parametric one. The difference is
rooted in the degree of prior knowledge we assumed with the population model
in both methods: In the parametric method, we assume the function forms which
are exactly the same as the underlying population model that was used to
generate simulated catalogue. It represents a too optimistic extreme. For in
reality, we would never be able to guess a perfect function representation of
the underlying population. On the other hand, in the nonparametric method, we
do not impose any assumption on the properties of the population model. It
represents a conservative extreme. For in reality, we more or less know the
general shape or trends of the population.

Figure 8: The merger rate as function of the redshift, reconstructed with
equation (11) from the synthetic aLIGO, ET and half-ET catalogues. Figure 9:
Two-dimensional theoretical merger rate as function of redshift and primary
masses. The black region marks the region in the parameters space that is not
reachable by aLIGO (single aLIGO detector at its nominal sensitivity). Figure
10: The merger rate as function of the primary BH mass, reconstructed with
equation (11) from the synthetic aLIGO, ET and half-ET catalogues.

##  4 Conclusion and Discussion

In this paper, we describe the details of how the GW-Universe Toolbox
simulates observations of BBH mergers with different ground-based GW
detectors. As a demonstration of its scientific application, we use the
synthetic catalogues to reconstruct and constrain the cosmic merger rate and
mass function of the BBH population. The comparison amongst the results from
one month (half) ET observation and ten years aLIGO observations shows the
overwhelming advantages of the 3rd generation detectors over the 2nd
generation, especially at a redshift higher than $\sim 2$, around where the
cosmic merger rate is believed to peak.

The reconstruction and constraining on $\mathcal{R}(z)$ and $p(m_{1})$ are
performed with two methods, namely, 1) a Bayesian method, assuming a
parameteric formula of $\mathcal{R}(z)$ and $p(m_{1})$ and 2) a KDE method
with non-parameteric $\mathcal{R}(z)$ and $p(m_{1})$. For the catalogue of ET
observations, the results from both methods are qualitatively the same.
However, for aLIGO, the Bayesian method can put constraints on
$\mathcal{R}(z)$ beyond the detecting limit, where the non-parameteric method
lost its ability completely. The $m_{1}$ dependence is also much more
accurately reconstructed by the Bayesian method than the non-parameteric
method for aLIGO’s catalogue. The difference is due to the extra information
placed by assuming a specific parameterisation of the population model in the
Bayesian method. In our example, the underlying $\mathcal{R}(z)$ and
$p(m_{1})$ used to generate the catalogues are the same with those assumed in
the Bayesian method; while in the KDE method, no assumption of the general
shape of $\mathcal{R}(z)$ and $f(m_{1})$ are placed, not even the assumption
of its smoothness. These represent the two extreme situations in the
differential merger rate reconstruction.

In reality, the underlying shape of $\mathcal{R}(z)$ and $f(m_{1})$ are
unknown and they could be very different with what we used. The parameterised
Bayesian method may give biased inference and underestimated uncertainties.
Furthermore, any unexpected structures in $\mathcal{R}(z)$ and $f(m_{1})$
cannot be recovered with the parameterised Bayesian method. For instance,
there can be additional peaks in the BH spectrum corresponding to pulsational
pair-instability supernovae and primordial BHs. On the other hand, although
the function form of the merger rate is unknown, the general trends should be
more or less known. A better reconstruction method should lie in between these
two extreme methods.

An alternative method is to parameterise $\mathcal{R}(z)$ and $p(m_{1})$ as
piecewise constant functions and use the Bayesian inference to obtain the
posterior of the values in each bin, as in Vitale et al. (2019). It has the
advantage that no assumptions on the shapes of $\mathcal{R}(z)$ and $p(m_{1})$
is needed, like the Kernel density estimator method. The disadvantage of this
method is that MCMC sampling in high dimension needs to be performed. If one
assumes that the $p(m_{1})$ is independent of $z$, there are more than $N+M$
free parameters, where $N$ and $M$ are the number of bins in $z$ and $m_{1}$
respectively. If we want the resolutions to be $0.5$ in $z$ and $1\,M_{\odot}$
in $m_{1}$, then the dimension of the parameter space is $>80$ for ET
catalogue. In the more general case, when one includes the dependence of
$p(m_{1})$ on $z$, the number of parameters becomes $N\times M$, therefore the
dimension of the parameter space goes catastrophically to $>1200$.

Another conclusion we draw about the early phase ET (half-ET) is that, half-ET
can also give a good constraint on $\mathcal{R}(z)$ at high redshift Universe.
However, both $\mathcal{R}(z)$ and primary mass function are less accurately
constrained compared with the full ET, even using a larger number of detected
events from a longer observation duration. The reason is that the full ET
determines smaller uncertainties on the physical parameters of the GW events.

In this paper, we compared the results from a single aLIGO detector with that
of ET. In reality, there will always be multiple detectors working as a
network (e.g. LIGO-Virgo-KAGRA), which is expected to have better results than
that a single aLIGO. However, we do not expect a qualitatively different
conclusion from that ET (as a representative of the 3rd generation GW
detectors) has an overwhelming ability in constraining the BBH population
properties over the 2nd generation detectors, even in a form of a network. We
also anticipate that the same conclusion applies to other 3rd generation GW
detectors, e.g., the Cosmic Explorer.

## References

  * Abbott et al. (2016) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, Phys. Rev. Lett., 116, 061102. doi:10.1103/PhysRevLett.116.061102
  * Abbott et al. (2019a) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019, ApJ, 882, L24. doi:10.3847/2041-8213/ab3800
  * Abbott et al. (2019b) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019, Physical Review X, 9, 031040
  * Abbott et al. (2020) Abbott, R., Abbott, T. D., Abraham, S., et al. 2021, Physical Review X, 11, 021053. doi:10.1103/PhysRevX.11.021053
  * Abbott et al. (2021a) Abbott, R., Abbott, T. D., Abraham, S., et al. 2021, ApJ, 913, L7. doi:10.3847/2041-8213/abe949
  * Abbott et al. (2021b) Abbott, R., Abbott, T. D., Abraham, S., et al. 2021, ApJ, 913, L7. doi:10.3847/2041-8213/abe949
  * Arca Sedda et al. (2021) Arca Sedda, M., Mapelli, M., Benacquista, M., et al. 2021, arXiv:2109.12119
  * Ajith et al. (2011) Ajith, P., Hannam, M., Husa, S., et al. 2011, Phys. Rev. Lett., 106, 241101
  * Banerjee (2021) Banerjee, S. 2022, Phys. Rev. D, 105, 023004. doi:10.1103/PhysRevD.105.023004
  * Belczynski et al. (2002) Belczynski, K., Kalogera, V., & Bulik, T. 2002, ApJ, 572, 407. doi:10.1086/340304
  * Belczynski et al. (2010) Belczynski, K., Dominik, M., Bulik, T., et al. 2010, ApJ, 715, L138. doi:10.1088/2041-8205/715/2/L138
  * Belczynski et al. (2016) Belczynski, K., Holz, D. E., Bulik, T., et al. 2016, Nature, 534, 512. doi:10.1038/nature18322
  * Carilli et al. (2008) Carilli, C. L., Lee, N., Capak, P., et al. 2008, ApJ, 689, 883. doi:10.1086/592319
  * Chruślińska et al. (2021) Chruślińska, M., Nelemans, G., Boco, L., et al. 2021, MNRAS, 508, 4994. doi:10.1093/mnras/stab2690
  * Corral-Santana et al. (2016) Corral-Santana, J. M., Casares, J., Muñoz-Darias, T., et al. 2016, A&A, 587, A61
  * Curtis-Lake et al. (2013) Curtis-Lake, E., McLure, R. J., Dunlop, J. S., et al. 2013, MNRAS, 429, 302. doi:10.1093/mnras/sts338
  * DiNardo et al. (1996) DiNardo, J., Fortin, N., Lemieux, T. 1996, Econometrica, 64, 5, 1001-44
  * Doherty et al. (2006) Doherty, M., Bunker, A., Sharp, R., et al. 2006, MNRAS, 370, 331. doi:10.1111/j.1365-2966.2006.10484.x
  * Farr et al. (2019) Farr, W. M., Fishbach, M., Ye, J., et al. 2019, ApJ, 883, L42
  * Farmer et al. (2020) Farmer, R., Renzo, M., de Mink, S. E., et al. 2020, ApJ, 902, L36. doi:10.3847/2041-8213/abbadd
  * Foreman-Mackey et al. (2014) Foreman-Mackey, D., Hogg, D. W., & Morton, T. D. 2014, ApJ, 795, 64. doi:10.1088/0004-637X/795/1/64
  * Goto et al. (2011) Goto, T., Arnouts, S., Inami, H., et al. 2011, MNRAS, 410, 573. doi:10.1111/j.1365-2966.2010.17466.x
  * Hild et al. (2011) Hild, S., Abernathy, M., Acernese, F., et al. 2011, Classical and Quantum Gravity, 28, 094013. doi:10.1088/0264-9381/28/9/094013
  * Iyer et al. (2011) Iyer, B., et al. 2011, URL https://dcc.ligo.org/LIGO-M1100296/public
  * Kagra Collaboration et al. (2019) Kagra Collaboration, Akutsu, T., Ando, M., et al. 2019, Nature Astronomy, 3, 35. doi:10.1038/s41550-018-0658-y
  * Kennicutt & Evans (2012) Kennicutt, R. C., & Evans, N. J. 2012, ARA&A, 50, 531
  * Kovetz (2017) Kovetz, E. D. 2017, Phys. Rev. Lett., 119, 131301
  * Kovetz et al. (2017) Kovetz, E. D., Cholis, I., Breysse, P. C., et al. 2017, Phys. Rev. D, 95, 103010
  * Kruckow et al. (2018) Kruckow, M. U., Tauris, T. M., Langer, N., et al. 2018, MNRAS, 481, 1908. doi:10.1093/mnras/sty2190
  * Ly et al. (2011) Ly, C., Lee, J. C., Dale, D. A., et al. 2011, ApJ, 726, 109. doi:10.1088/0004-637X/726/2/109
  * Madau & Dickinson (2014) Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415
  * Maiolino & Mannucci (2019) Maiolino, R. & Mannucci, F. 2019, A&A Rev., 27, 3. doi:10.1007/s00159-018-0112-2
  * Magdis et al. (2010) Magdis, G. E., Elbaz, D., Daddi, E., et al. 2010, ApJ, 714, 1740. doi:10.1088/0004-637X/714/2/1740
  * Mapelli et al. (2017) Mapelli, M., Giacobbo, N., Ripamonti, E., et al. 2017, MNRAS, 472, 2422. doi:10.1093/mnras/stx2123
  * Mapelli & Giacobbo (2018) Mapelli, M. & Giacobbo, N. 2018, MNRAS, 479, 4391. doi:10.1093/mnras/sty1613
  * Mapelli et al. (2021) Mapelli, M., Bouffanais, Y., Santoliquido, F., et al. 2022, MNRAS, 511, 5797. doi:10.1093/mnras/stac422
  * Punturo et al. (2010) Punturo, M., Abernathy, M., Acernese, F., et al. 2010, Classical and Quantum Gravity, 27, 194002. doi:10.1088/0264-9381/27/19/194002
  * Regimbau et al. (2012) Regimbau, T., Dent, T., Del Pozzo, W., et al. 2012, Phys. Rev. D, 86, 122001. doi:10.1103/PhysRevD.86.122001
  * Reitze et al. (2019) Reitze, D., Adhikari, R. X., Ballmer, S., et al. 2019, Astro2020: Decadal Survey on Astronomy and Astrophysics, APC white papers, no. 35; Bulletin of the American Astronomical Society, Vol. 51, Issue 7, id. 35 (2019)
  * Rujopakarn et al. (2010) Rujopakarn, W., Eisenstein, D. J., Rieke, G. H., et al. 2010, ApJ, 718, 1171. doi:10.1088/0004-637X/718/2/1171
  * Rodriguez et al. (2013) Rodriguez, C. L., Farr, B., Farr, W. M., et al. 2013, Phys. Rev. D, 88, 084013
  * Singh et al. (2021) Singh, N., Bulik, T., Belczynski, K., et al. 2021, arXiv:2112.04058
  * Scott (1987) Scott, D. W., and Terrell, G. R., 1987, Journal of the American Statistical Association 82, 1131–46.
  * Scott (1992) Scott, D.W. (1992), Multivariate Density Estimation: Theory, Practice and Visualization. John Wiley & Sons, Inc., New York. http://dx.doi.org/10.1002/9780470316849
  * Sedgwick et al. (2011) Sedgwick, C., Serjeant, S., Pearson, C., et al. 2011, MNRAS, 416, 1862. doi:10.1111/j.1365-2966.2011.19156.x
  * Serjeant et al. (2002) Serjeant, S., Gruppioni, C., & Oliver, S. 2002, MNRAS, 330, 621. doi:10.1046/j.1365-8711.2002.05109.x
  * Silverman (1986) Silverman, B.W., 1986, Chapman and Hall/CRC.
  * Thompson et al. (2006) Thompson, R. I., Eisenstein, D., Fan, X., et al. 2006, ApJ, 647, 787. doi:10.1086/505568
  * van Son et al. (2021) van Son, L. A. C., de Mink, S. E., Callister, T., et al. 2021, arXiv:2110.01634
  * Vitale & Evans (2017) Vitale, S., & Evans, M. 2017, Phys. Rev. D, 95, 064052
  * Vitale et al. (2019) Vitale, S., Farr, W. M., Ng, K. K. Y., et al. 2019, ApJ, 886, L1
  * Yen-Chi Chen (2017) Yen-Chi Chen, 2017, Biostatistics & Epidemiology, 1, 161-187
  * You et al. (2020) You, Z.-Q., Zhu, X.-J., Ashton, G., et al. 2021, ApJ, 908, 215. doi:10.3847/1538-4357/abd4d4
  * Yi et al. (2021) Yi, S.-X., Nelemans, G., Brinkerink, C., et al. 2021, arXiv:2106.13662

##  Appendix A Kernel density estimation

Kernel density estimation (KDE) is one of the most famous method for density
estimation. For the univariate case is defined as:

$\hat{f}(x,h)=\frac{1}{nh}\sum_{i=1}^{n}K(\frac{x-X_{i}}{h}),$ (13)

where $K$ is the chosen kernel and $h$ is the so-called bandwidth parameter.
In the case of a Gaussian kernel, $h$ can be thought as the standard deviation
of the normal density centred on each data point. The choice of the kernel
doesn’t really affect the density estimate, if the number of observations is
high enough; what does instead affect it is the choice of the bandwidth
parameter $h$, that has to be estimated with one of the available methods. The
most common methods are rule-of-thumb like Scott’s rule (Scott (1992))
$H=n^{-2/(m+4)}\Sigma$ and Silverman’s rule (Silverman (1986))
$H=\frac{4}{n(m+2)}^{2/(m+4)}\Sigma$, where $m$ is the number of the
variables, n is the sample size and $\Sigma$ is the empirical covariance
matrix. A more accurate but computationally expensive method is Cross-
validation bandwidth selection described in Scott (1987).

As introduced in DiNardo et al. (1996), KDE can be adapted to the case where
each observation is attached to a weight that reflects its contribution to the
sample. The weighted kernel density estimator can be written as:

$\hat{f}(x,h)=\frac{1}{nh}\sum_{i=1}^{n}w_{i}K(\frac{x-X_{i}}{h}),$ (14)

to ensure that the result of the density estimation has still the properties
of a density, it is necessary to normalize the weights so that
$\sum_{i=1}^{n}w_{i}=1$.

KDE is easily extendable to the multivariate case, and used to estimate
densities $f$ in $R^{p}$ based on the same principle: perform an average of
densities centred at the data points. For the multivariate case, kernel
density estimation is defined as:

$\hat{f}(x,H)=\frac{1}{n|H|^{\frac{1}{2}}}\sum_{i=1}^{n}K(H^{-\frac{1}{2}}(x-X_{i})),$
(15)

where $K$ is the multivariate kernel, a p-variate density that depends on the
bandwidth matrix $H$, a $p\times p$ symmetric and positive definite matrix.