\section{\label{sec:motivation}Motivation}
The phenomenon of Malmquist selection biases in brightness-limited astronomical surveys has been known for nearly a century~\citep{Malmquist1922, Malmquist1925}, and methods for mitigating this bias have remained an active area of research in astronomy, e.g.,~\cite{Tiwari2018, Farr2019}.
Gravitational-wave searches are affected by Malmquist bias as the parameters of the merging binaries determine the luminosity of the source.
More massive binaries produce larger strain, all else equal, and so they can be observed at greater distances, at least until the mass becomes so large that the signal begins to shift out of the observing band.
Additionally, binaries with black holes that are spinning aligned to the angular momentum are visible to further distances~\citep{Campanelli2006, Ng2018}.

Over the past several years, there has been increasing interest in population studies, which seek to measure the distribution of astrophysical parameters such as the mass, spins, and distance of merging compact objects using events observed with Advanced LIGO/Virgo~\citep{aLIGO, AdVirgo}, e.g.,~\cite{O3aPop} and references therein.
To perform unbiased inference on the distribution of astrophysical parameters, it is necessary to account for selection biases when performing population inference, see, e.g.,~\cite{Loredo2004, Farr2015,Mandel2019,Thrane2019,Vitale2020}.
The standard method employed in gravitational-wave searches requires computing the total observed spacetime volume for the population model.
Evaluating the observed volume for different population parameters involves integrating the observable volume as a function of all binary parameters and time over the whole observing time.

Astrophysical population inference is performed using Markov Chain Monte Carlo (MCMC)~\citep{Metropolis1953, Hastings1970} or nested sampling algorithms~\citep{Skilling04}, which typically require $O(10^6)$ likelihood evaluations for a well converged run, and until now this integral has been evaluated on the fly at every iteration.
As the binary black hole catalog grows this integral must be evaluated with increasing precision~\citep{Farr2019}, and correspondingly increased computational cost.
Additionally, as the catalog grows, so does our resolving power, meaning that sub-dominant effects, e.g., the effect of black hole spin on the observed volume, must be considered.
In this work, we apply a combination of supervised and unsupervised machine learning methods to increase the precision and accuracy of our sensitive volume estimates with manageable computation cost.

The rest of this paper is structured as follows.
In the next section, we outline methods for accounting for Malmquist bias in gravitational-wave searches.
We then briefly summarize some preliminaries for gravitational-wave population inference in Section~\ref{sec:popinf}.
Following this, in Section~\ref{sec:density}, we describe the problem of density estimation and discuss various commonly used methods.
In Sections~\ref{sec:implementation} and~\ref{sec:population-sensitivity}, we estimate the gravitational-wave transient selection function as a function of binary parameters and population hyper-parameters using a Gaussian mixture model and a neural network respectively.
After this, we apply our methods to the binary black hole systems identified in~\cite{GWTC2}.
Some closing thoughts are then provided.

\section{\label{sec:vt}Sensitive volume estimation}
Most gravitational-wave population analyses impose a detection threshold on the analyzed triggers to avoid contamination from terrestrial noise sources.
Applying this threshold leads to a selection bias in the observed sample.
We quantify this by considering the probability that a signal with parameters $\theta$ would surpass our threshold $\rho_{th}$\footnote{The specific choice of threshold is irrelevant so long as it is robustly defined.}
\begin{equation}
    p_{\text{det}}(\theta)
    = \int_{\rho > \rho_{\text{th}}} d{\bf d} \, p({\bf d}|\theta).
\label{eq:pdet}
\end{equation}
The integral is over all observed data surpassing the threshold.
For population analyses, we require the fraction of all sources which are detectable, for a given population model, characterized by parameters, $\Lambda$,
\begin{equation}
        P_{\text{det}}(\Lambda) = \int \, d\theta \, p(\theta | \Lambda) \, p_{\text{det}}(\theta),
\label{eq:Pdet}
\end{equation}
where $p(\theta|\Lambda)$ is a conditional prior for $\theta$ given population (hyper-) parameters $\Lambda$.
For a detailed derivation of these quantities see, e.g.,~\cite{Farr2015, Tiwari2018, Thrane2019, Mandel2019}.

We emphasize that all population analyses which apply a threshold necessarily have a corresponding selection bias that must be accounted for, including analyses that explicitly model contamination of the sample from terrestrial sources~\citep{Gaebel2019, Galaudage2020, Roulet2020}.
However, see~\cite{TBS-PE} for a method that avoids thresholds entirely.

The integral over ${\bf d}$ in Eq.~\ref{eq:pdet} requires that we understand the sensitivity of gravitational-wave searches throughout the observing history.
In practice, there are currently two widely used methods to compute this integral: inject simulated signals into the data and see how many of them are recovered by the search pipelines; or use a semi-analytic approximation based on the power spectral density of the interferometers, e.g.,~\citep{Finn1993}.
However, there have been several recent methods to leverage machine learning methods to estimate $p_\text{det}$~\citep{Gerosa2020, Wong2020b}.

The former method gives the most faithful representation of the search sensitivity.
However, the latter has several computational advantages.
Because of the large parameter space which must be covered, the injection and recovery procedure gives us only the parameter values of the found/missed signals, whereas the semi-analytic approach can efficiently generate a numerical value for $p_{\det}$ marginalized over specific nuisance parameters.
Thus, the semi-analytic approach can also be performed much more computationally cheaply due to the cost of performing and recovering injections.
Previous methods to improve the reliability of semi-analytic estimates include calibration of semi-analytic estimates with the output of injection campaigns~\citep{WysockiVT}.
In this work, we use density estimation on the set of found injections to provide a continuous, generative, model for $p_\text{det}$.

The integral over $\theta$ in Eq.~\ref{eq:Pdet} marginalizes over all the parameters describing the source---15 parameters to completely characterize a quasi-circular binary black hole merger.
In practice, many of the parameters are not modeled in current population analyses, the most complex models considered currently fit for the distribution of seven of these parameters, the two component masses, spin magnitudes, spin-tilt angles, and redshift, requiring the evaluation of a seven-dimensional integral within each likelihood evaluation.
The other parameters are assumed to be well described by the prior distributions used during sampling.
These are mostly geometric parameters describing the position and orientation of the binary, although it is possible that some of these parameters may deviate from isotropy.
For example, we could search for deviations from isotropy over the sky position, e.g.~\cite{Payne2020, Stiskalek2021}, or features in the distribution of the azimuthal spin parameters due be influenced by spin-orbit responances~\citep{Schnittman2004, Gerosa2018b}.
This integral is, therefore, recast as a Monte Carlo integral over the set of found injections~\citep{Tiwari2018, Farr2019}
\begin{equation}
    P_{\text{det}}(\Lambda) = \frac{1}{N_{\text{inj}}} \sum_{i=1}^{N_{\text{found}}} \frac{p(\theta_i | \Lambda)}{p(\theta_i | \Lambda_0)}.
\label{eq:Pdetmc}
\end{equation}
Here, $p(\theta|\Lambda_0)$ is the distribution of the injected signals, $N_{\text{inj}}$ is the total number of injected signals and $N_{\text{found}}$ is the number of found injections.
To ensure sufficient convergence of the Monte Carlo integral we must have an effective sample size of at least four times the number of observed events~\citep{Farr2019}.
This means that to fit tightly peaked distributions we need a large number of samples for the distribution of found injections or a continuous representation of $p_{\text{det}}$.
Performing more injections to increase the number of recovered injections quickly becomes computationally prohibitive.
In this work, we resolve this issue by performing density estimation using the set of found injections.
Using these density estimates, we can directly evaluate $p_{\text{det}}$ and/or generate additional samples from the distribution of found injections.

\section{\label{sec:popinf}Population inference}

\subsection{\label{sec:models}Models}

For demonstration purposes, we consider two simple population models from within the gravitational-wave literature.
Following~\cite{O3aPop, Talbot2018a}, we model the binary black hole mass distribution as a power law in the larger mass, $m_1$, between the minimum and maximum mass along with a normally distributed component and a power law in the mass ratio, $q = m_2 / m_1$,
\begin{align}
    p(m_1 &| \alpha, m_{\min}, m_{\max}) = (1 - \lambda) \frac{(1-\alpha) m_1^{-\alpha}}{m_{\max}^{1-\alpha} - m_{\min}^{1-\alpha}} \nonumber\\
    &\quad + \frac{\lambda}{\sqrt{2\pi\sigma^2_m}} \exp\left( - \frac{(m_1 - \mu_m)^2}{2\sigma^2_m} \right) \\
    p(q &| m_1, \beta, m_{\min}) =  \frac{(1+\beta) q^{\beta}}{1 - \left(\frac{m_{\min}}{m_1}\right)^{1+\beta}}\label{eq:mass_ratio}.
\end{align}
This is the \textsc{Power-Law + Peak} model in~\cite{O3aPop} without the low-mass smoothing.
We assume that both component spins are drawn from the same distribution.
We model the distribution of spin magnitudes as following a Beta distribution~\citep{Wysocki2019}
\begin{equation}
    p(a_i | \alpha_{\chi}, \beta_{\chi}) = \frac{a_i^{\alpha_{\chi} - 1} (1 - a_i)^{\alpha_{\chi} - 1}}{\mathrm{B}(\alpha_{\chi}, \beta_{\chi})}.
\end{equation}
We model the distribution of spin orientations as a combination of a truncated half-normal and a uniform distribution~\citep{Talbot2017}
\begin{equation}
    p(\cos\theta_i | \sigma_{s}) = \frac{(1 - \xi)}{2} + \xi N(\sigma_s) \exp\left( - \frac{(\cos\theta_i - 1)^2}{2 \sigma^{2}_{s}} \right)
\end{equation}
The factor $N$ ensures that the distribution is properly normalized.
This is the \textsc{Default} model in \cite{O3aPop}.

The reference distribution is the product of these distributions
with population hyper-parameters $\alpha=1.35$, $m_{\min}=2$, $m_{\max}=100$, $\lambda=0$, $\beta=2$, $\alpha_\chi=\beta_\chi=1$, and $\sigma_{s}=\infty$.
We note however that the injections performed in~\cite{O3aPop} did not include the in-plane components of the spins responsible for spin-induced precession of the orbital plane and only include the spin aligned with the orbital angular momentum, $\chi = a \cos\theta$.
The distribution of the aligned-spin component is~\citep{Lange2018}
\begin{equation}
    p(\chi_i | \Lambda_0) = -\frac{1}{2} \ln \left| \chi_i \right|.
\end{equation}

\subsection{\label{sec:poplike}Likelihood}

The standard likelihood used in population inference for gravitational-wave sources in the presence of selection biases is, e.g.,~\citep{Thrane2019, Mandel2019, Vitale2020},
\begin{equation}
    \mathcal{L}(\{d_i\}|\Lambda) = \frac{1}{P_{\text{det}}(\Lambda)^N} \prod^{N}_{i} \int d\theta_i \mathcal{L}(d_i | \theta_i) p(\theta_i | \Lambda).
\label{eq:likelihood}
\end{equation}
Where the product over $i$ runs over the $N$ events with data $d_i$.
The integral over $\theta_i$ is typically performed by importance sampling from the single-event posterior distribution for $p(\theta_i | d_i)$ as is done to calculate $P_{\text{det}}$.
We take the publicly available samples from the single-event posterior distributions from~\cite{GWTC1, GWTC2}.

\section{\label{sec:density}Density estimation}

Reconstructing a function or probability density from a finite set of samples from the distribution is a widespread problem in data analysis.
For example, injection campaigns to determine the sensitivity of gravitational-wave detectors do not give us a continuous description of the sensitivity, but rather a discrete set of samples from the distribution of found injections.

Density estimation methods can be loosely divided into parametric and non-parametric methods.\footnote{We note that the word ``non-parametric'' is something of a misnomer, as often these models involve large numbers of unphysical parameters to perform the fit. An alternative delineation is between models where the parameters are either physically motivated (parametric) or not physically motivated (non-parametric).}
Parametric density estimation involves fitting a parameterized phenomenological model to the data.
An example of this is the method used to reconstruct the population distribution of binary black holes in this work.
To estimate the gravitational-wave selection function we will rely on non-parametric density estimation.

Many methods for nonparametric density estimation are commonly used, however, most traditional methods such as binning or kernel density estimation scale poorly as the dimensionality of the problem increases.
More sophisticated density estimation techniques involving the optimization of many parameters, such as Gaussian mixture models or flow-based inference, have proved successful at approximating complex functions in large dimensional spaces, see, e.g.,~\cite{Powell2019, Gabbard2019, Green2020a, Green2021, Wong2020a, Wong2020b, Wong2021a} for applications in gravitational-wave inference.
These models also provide natural ways to generate additional samples from the underlying densities and are therefore sometimes referred to as generative models.

In this work, we approximate $p_{\text{det}}(\theta)$ using a Gaussian mixture model.
A Gaussian mixture model approximates the distribution as a set of multivariate Gaussian distributions each with a unique mean and covariance.
The number of components is a hyper-parameter that can be manually tuned.
Expectation maximization is performed by varying the means and covariances of each component.

\section{\label{sec:implementation}Sensitivity to binaries}

In this section and the next, we develop methods to evaluate Equations~\ref{eq:Pdet} and~\ref{eq:pdet} and generate new samples from the distribution of found injections.
We begin by training a function to estimate $p_{\text{det}}(\theta) p(\theta | \Lambda_0)$ using a set of $\sim 5\times 10^4$ found injections.
We take the found injections from the first half of Advanced LIGO/Virgo's third observing run~\citep{O3aPop, GWTC2} and samples generated assuming a PSD representative of the first two observing runs.

\subsection{\label{sec:preprocessing}Pre-processing}

Gravitational-wave parameters are typically only defined over finite domains and many have significant support at the edges, e.g., spin magnitudes are contained in the unit interval and the majority of observed black holes are consistent with being non-spinning.
However, the algorithms we use for density estimation work best over an infinite domain without sharp boundaries.
Therefore, we begin by performing the following mapping to the found injections:
\begin{enumerate}
    \item Transform the injections from the original distribution to the unit interval.
    \item Map the samples from the unit interval to a unit normal distribution using the probit function $\Phi^{-1}$~\citep{Bliss1934}.
\end{enumerate}

We considered various choices for the initial mapping onto the unit interval.
The simplest method is a simple shift and scale of the original domain onto the unit interval.
The next method we employ is mapping the samples to the cumulative distribution function of $p(\theta|\Lambda_0)$ denoted $F(\theta)$.
The final method we consider is mapping the samples to the cumulative distribution of an analytic approximation to the distribution of found injections $\tilde{p}_{\text{det}}(\theta)p(\theta | \Lambda_0)$.

We use the second option to maximize the simplicitly in evaluating $p_{\text{det}}$.
Mathematically the full transformation is
\begin{equation}
    \theta^{'} = \Phi^{-1}\left(F(\theta)\right) .
\label{eq:regularisation}
\end{equation}
The Jacobian for this transformation is given by
\begin{align}
    \mathcal{J}(\theta) = \frac{d \theta'}{d\theta} = \frac{p(\theta | \Lambda_0)}{\mathcal{N}(\theta';\mu=0,\sigma=1)},
\end{align}
where $\mathcal{N}$ is the normal distribution.

In Figure~\ref{fig:data}, we show the set of found injections in the original ($\theta$) space (top) and the transformed ($\theta'$) space.
The transformation has removed the railing against the boundaries in all the parameters.

\begin{figure*}
    \centering
    \includegraphics[width=\linewidth]{scaling}
    \caption{
    Two-dimensional density plots of the distribution of found injections before (top) and after (bottom) applying the transformation described in Section~\ref{sec:preprocessing}.
    The pre-processing simplifies the density estimation procedure by removing the sharp boundaries.
    }
    \label{fig:data}
\end{figure*}

\subsection{\label{sec:training}Density estimation}

We then take these regularized samples $\{\theta^{'}\}$ and split them into training (80\%) and test (20\%) sets for our estimator.
Training is performed by maximizing the mean natural log-likelihood of the test samples.
For our Gaussian mixture model, we consider models with up to 20 components.
We find that the quality of the fit plateaus when using 10 components and so use a 10 component model throughout.

\subsection{Evaluation}

After training, we can trivially generate new samples from the target distribution by drawing samples $\theta'$ from a unit-normal distribution and apply the inverse of Eq.~\ref{eq:regularisation}.
Alternatively, since these are density estimates, we can also directly evaluate the estimated density
\begin{equation}
    \mathcal{D}(\theta') = \frac{p_{\text{det}}(\theta) p(\theta | \Lambda_0)}{ \mathcal{J}(\theta)}.
\end{equation}
In practice, we want to evaluate the selection function
\begin{equation}
    p_{\text{det}}(\theta) = \frac{\mathcal{J}(\theta) \mathcal{D}(\theta')}{p(\theta | \Lambda_{0})}.
\label{eq:pdet_density}
\end{equation}
We use Eq.~\ref{eq:pdet_density} as an alternate means of computing Eq.~\ref{eq:Pdet} with an equivalent Monte Carlo integral over samples from the population distribution
\begin{equation}
    P_{\text{det}}(\Lambda) = \left< p_{\text{det}}(\theta_i) \right>_{\theta_i \sim p(\theta | \Lambda)}.
    \label{eq:Pdet_density}
\end{equation}

If we perfomed the initial map using the injected cumulative distribution function we find
\begin{equation}
    p_{\text{det}}(\theta) = \frac{\mathcal{D}(\theta')}{\mathcal{N}(\theta';\mu=0,\sigma=1)}.
\label{eq:pdet_density_sensible}
\end{equation}

We note that this method requires an efficient method of generating samples from the population distribution.
This can be trivially performed using inverse-transform sampling if the population model has an analytically invertible cumulative distribution function, or is a sum of such distributions.
For other population models such as those using the low-mass smoothing introduced in~\cite{Talbot2018a} or the redshift distribution introduced in~\cite{Fishbach2018} and used in many population analyses, e.g.,~\citep{O2Pop, O3aPop}, another method is required, see, e.g.~\cite{Wong2020a}.

In practice, we find that this method requires far fewer samples in the Monte Carlo integral than when resampling the found injections, 5000 samples from the population model versus 50000 found injections with the same number of effective samples for each method.
The number of effective samples is defined slightly differently for the two Monte Carlo methods considered here.
For Equation~\ref{eq:Pdet_density} we adopt the usual definition~\citep{Elvira2018}
\begin{equation}
    N_{\text{eff}} = \frac{\left(\sum_{i=1}^{N} p_{\text{det}}(\theta_i)\right)^2}{\sum_{i=1}^{N} p_{\text{det}}(\theta_i)^2}.
\end{equation}
However, for Equation~\ref{eq:Pdetmc} a correction is required to account for the initial injections with $p_{\text{det}}=0$~\citep{Farr2019}.
Following~\cite{Farr2019}, for both of these methods we only allow for samples with $N_{\text{eff}} > 4 N_{\text{events}}$ and marginalize over the statistical uncertainty in $P_{\text{det}}$ in the likelihood.

Another advantage of evaluating Equation~\ref{eq:Pdet_density} is that we can compute the aligned-spin-dependent $P_{\text{det}}$ when employing a parameterization of the spin magnitude and orientation.
This is not possible when using Equation~\ref{eq:Pdetmc} as there is no unique mapping from the set of found injections to the full component spins, as injections were not performed with in-plane spin components.

\section{\label{sec:population-sensitivity}Sensitivity to populations}

In the previous sections, we described two methods to evaluate $P_{\text{det}}(\Lambda)$ using Monte Carlo integrals, one by performing a sum over found injections and one by performing a sum over samples from the population model.
However, both of these methods can require large numbers of samples to converge and these samples can be computationally costly to generate.
In this section, we use a neural network regressor as a high-dimensional interpolant of $P_{\text{det}}(\Lambda)$.

Our training data are $5 \times 10^5$ samples from the prior distribution for $\Lambda$ with 60\% used for training, 20\% for testing, and 20\% for comparison in the following section.
For each sample in the training data, the label is calculated by evaluating $P_{\text{det}}(\Lambda)$ using Equation~\ref{eq:Pdet_density}.
To perform the integral, we generate 5000 samples from the population model at each iteration and evaluate $p_{\text{det}}$ using our Gaussian mixture model estimator.

We then train a fully connected neural network with three hidden layers, each with 64 nodes, to estimate $P_{\text{det}}$.
Training is performed by minimizing the mean absolute error in $\log(P_{\text{det}})$ using the Adam optimization algorithm.

\section{\label{sec:results}Results}

\subsection{\label{sec:validation}Computing $P_{\text{det}}$}

While the mean log-likelihood is a suitable metric for training our estimators, we can perform a stronger test on the samples by considering how well the population-averaged $P_{\text{det}}(\Lambda)$ compares across a range of $\Lambda$ values using our Gaussian mixture model.
We compare our estimators using 100000 samples from the prior distribution for $\{\Lambda^{\pi}_i\}$.
We evaluate $P_{\text{det}}$ four times for each this set of samples.
\begin{enumerate}
    \item Using Equation~\ref{eq:Pdetmc} with the original found injection set.
    \item Using Equation~\ref{eq:Pdetmc} with samples generated using our Gaussian mixture model.
    \item Using Equation~\ref{eq:Pdet_density} with samples from the population model.
    \item Using our neural network.
\end{enumerate}

We find that the dependence of the selection function on the parameters describing the mass distribution is indistinguishable between the different methods.
When using Equation~\ref{eq:Pdetmc} the $P_{\text{det}}$ does not depend on the spin distribution by construction.
We, therefore, see a difference in the dependence of the selection function on the spin parameters between these two cases.

\subsection{\label{sec:hyperpe}Population inference}

\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{spectrum}
    \caption{
    The inferred distribution of primary mass (top) and cosine spin tilt (bottom) when analyzing the binary black hole mergers identified in GWTC-2 using four different methods for estimating the observed volume, $P_{\text{det}}$:
    (blue) using the original found injections, (orange) new samples generated using a Gaussian mixture model fit to the found injections, (green) evaluating $p_{\text{det}}$ using samples from the population models, and (red) using the neural network interpolant.
    The results in blue and green do not include the dependence of $P_{\text{det}}$ on the spin distribution.
    We note that this leads to a slight change in the inferred distribution of spin orientations, when including spin dependence in the selection function there is slightly more preference for isotropic spins, although this is well below statistical uncertainties.
    The solid curves show the posterior predictive distribution and the shaded regions show the symmetric $95\%$ credible region.
    }
    \label{fig:hyperpe}
\end{figure}

We now consider the impact of our different $P_{\text{det}}$ evaluation methods on population inference.
Following~\cite{O3aPop}, we analyze the 44 confident binary black hole mergers with false alarm rate $<1~{\text{yr}^{-1}}$ from GWTC-2.
We perform population inference four times, once with each of the $P_{\text{det}}$ estimation methods described in Section~\ref{sec:validation}.

In Figure~\ref{fig:hyperpe}, we show the inferred astrophysical population of binary black holes when using our two methods to evaluate the selection function $P_{\text{det}}$.
In blue, we evaluate Equation~\ref{eq:Pdetmc} using the original simulated injections.
In orange, we evaluate Equation~\ref{eq:Pdetmc} using samples from our Gaussian mixture model.
In green, we evaluate Equation~\ref{eq:Pdet_density} using our Gaussian mixture model.
In red, we use our neural network interpolant.

We find no difference between the inferred mass distributions across the various methods.
Similarly, when neglecting the impact of spin on the selection function we see no difference between the methods.
When including the effect of the aligned component of the spins on the selection function, we find a modest change to the inferred distribution of spin orientations preferring an isotropic distribution of spin orientations.
We stress that the changes are within the current statistical uncertainties.

\section{\label{sec:discussion}Discussion}

Malmquist biases are ubiquitous in astronomical surveys, and methods of understanding and mitigating these biases are vital to performing astrophysical inference.
Calculating and mitigating this bias is typically done by performing large simulations where synthetic signals are injected into the data and counting the number of recovered signals.
While this method gives an optimal estimate of the performance of searches for signals in real data, they can be difficult to work with and extend to generic population models.
In this work, we apply unsupervised and supervised machine learning techniques to enable flexible, accurate, and efficient estimation of the selection function for gravitational-wave transient surveys.

First, we estimated the gravitational-wave selection function for individual events by employing unsupervised density estimation with a Gaussian mixture model.
To efficiently estimate this density we introduce a pre-processing step that improves the convergence of the density estimate.
Using this density estimate, we tested two methods to compute the population level selection function.
We found that our new density estimation method matches the previous injection resampling method for equivalent population models.

This method is more computationally expensive, especially for complex population models.
To mitigate this increased computational cost, we used the output of the first step and as a training set for supervised neural network regression to produce an interpolant for the population-level selection function.
Using this interpolant, we can quickly and accurately evaluate the population-level selection function for a wide range of population models.

Using our method, it is trivial to compute the fraction of sources which are observed $P_{\text{det}}$ by marginalizing over parameters other than those parameterized in the population model, e.g., evaluate $P_{\text{det}}$ using the parameters which most directly affect the sensitivity (chirp mass, mass ratio, effective aligned/precessing spin) and model the population in terms of parameters with the most intuitive physical meaning (component masses, spin magnitudes, and orientations).
Specifically, we model the distribution of black hole spins with a model for the magnitudes and orientations, while the injection campaign used to evaluate $P_{\text{det}}$ only considered the aligned-spin component of the black hole spins.
This allows us to extend the analysis presented in~\cite{O3aPop} to simultaneously measure the distribution of black hole masses and spins while considering the effect of aligned spin components on the selection function.

We leave a detailed analysis of the best combination of parameters to use for the density estimation to future work.
Our results are consistent with the results presented in~\cite{O3aPop}, however, the uncertainty on the measured selection function is less in most of the space when using our new method.
As the catalog grows it will be vital to understand the systematic error in our estimation of the selection function.

Machine learning methods for density estimation are rapidly gaining popularity in the gravitational-wave data analysis community, e.g.,~\cite{Powell2019, Gabbard2019, Green2020a, Green2021, Wong2021a, Wong2020a, Wong2020b, Cuoco2020}.
Most of these methods require the use of complex neural network-based density estimators which require tuning many more free parameters and thus extremely large training data sets.
The pre-processing method introduced here removes sharp spectral features, e.g., at prior boundaries, and thus enables high-precision estimation of the target distribution using Gaussian mixture models, rather than having to employ deep learning density estimators.
Combining this pre-processing in other density estimation problems may have a similarly simplifying effect.

One limitation of the current method is that the Gaussian mixture model employed in this work provides only a best-fit model and does not provide an indication of uncertainty in the fit over the parameter space.
We leave the exploration of density estimation techniques that model this uncertainty, e.g, Bayesian Gaussian mixture models to a future study.

\acknowledgments{
We thank Maya Fishbach for producing fake injections used in an early version of this work.
We thank Sylvia Biscoveanu, Tom Dent, Reed Essick, Cody Messick, Richard O'Shaughnessy, Alan Weinstein, Daniel Wysocki, and Salvatore Vitale for useful comments and discussions.
This work is supported through the Australian Research Council (ARC) Centre of Excellence CE170100004 and ARC Future Fellowship FT150100281.
This is document LIGO-P2000505.
This work was made possible using the following software: scikit-learn~\citep{scikit-learn}, tensorflow~\citep{tensorflow}, GWPopulation~\citep{Talbot2019}.
This research has made use of data, software, and/or web tools obtained from the Gravitational Wave Open Science Center (https://www.gw-openscience.org), a service of LIGO Laboratory, the LIGO Scientific Collaboration and the Virgo Collaboration.
Computing was performed computing clusters at the California Institute of Technology (LIGO Laboratory) supported by National Science Foundation Grants PHY-0757058 and PHY-0823459 and Swinburne University of Technology (OzSTAR).
This work used publicly available samples from~\cite{O2SampleRelease, O3aSampleRelease, O3aPopulationRelease, O3aSensitivityRelease}.
This work made use of Google Colaboratory.
}