Binary Vision: The Mass Distribution of Merging Binary Black Holes via Iterative Density Estimation

Our general approach to population inference can be considered as similar to maximum likelihood, with uncertainties quantified via empirical bootstrap methods (Efron 1979). Given a set of observed events, if we neglect measurement uncertainty in each event's parameters, our population estimate is a KDE where the kernel bandwidth for each event is adjusted (Breiman et al. 1977; Abramson 1982; Terrell & Scott 1992; Sain & Scott 1996) using an adaptive scheme to maximize the cross-validated likelihood (Sadiq et al. 2022). The adaptive-bandwidth KDE (Wang & Wang 2011) computes a density estimate $\hat{f}$ from observations X_i , i = 1, ..., n via

$$\begin{eqnarray}&&\hat{f}(x)={n}^{-1}\displaystyle \sum _{i=1}^{n}\displaystyle \frac{1}{h{\lambda }_{i}}K\left(\displaystyle \frac{x-{X}_{i}}{h{\lambda }_{i}}\right),\end{eqnarray} \tag{ 1 }$$

where K( · ) is the standard Gaussian kernel,

$$\begin{eqnarray}&&K(z)=\displaystyle \frac{1}{\sqrt{2\pi }}\exp \left(\displaystyle \frac{-{z}^{2}}{2}\right),\end{eqnarray} \tag{ 2 }$$

n is the total number of samples, and h λ_i takes the role of a local bandwidth, with h being the global bandwidth.

The local bandwidths are determined by first computing a pilot estimate ${\hat{f}}_{0}$ setting λ_i = 1 for all i, a standard fixed-bandwidth KDE; based on this pilot density, we then set

$$\begin{eqnarray}&&{\lambda }_{i}={\left(\displaystyle \frac{{\hat{f}}_{0}({X}_{i})}{g}\right)}^{-\alpha },\end{eqnarray} \tag{ 3 }$$

where α is the sensitivity parameter of the local bandwidth (0 < α ≤ 1) and g is a normalization factor,

$$\begin{eqnarray}&&\mathrm{log}g={n}^{-1}\displaystyle \sum _{i=1}^{n}\mathrm{log}{\hat{f}}_{0}({X}_{i}).\end{eqnarray} \tag{ 4 }$$

Finally the adaptive KDE $\hat{f}(x)$ is obtained by evaluating Equation (1) with the variable (local) bandwidth h λ_i . The equations are written for the case of one-dimensional data X_i , but the method may be applied in more dimensions by linearly transforming the data to have zero mean and unit covariance along each dimension and using an N-dimensional unit Gaussian kernel scaled by h λ_i .

The KDE hyperparameters, global bandwidth h, and sensitivity parameter α are determined by grid search using maximum likelihood as a figure of merit. A naïve "maximum likelihood" KDE is not well defined, as the likelihood increases indefinitely in the limit of small-bandwidth kernels centered on the observations, i.e., delta functions (Silverman 1986, Section 2.8). We prevent this collapse and address the bias–variance tradeoff via leave-one-out cross-validation (Silverman 1986, Section 3.4.4; see also Hastie et al. 2001). The cross-validated (log) likelihood is

$$\begin{eqnarray}&&\mathrm{log}{{ \mathcal L }}_{\mathrm{LOO}}=\displaystyle \sum _{i=1}^{n}\mathrm{log}{\hat{f}}_{\mathrm{LOO},i}({X}_{i}),\end{eqnarray} \tag{ 5 }$$

where ${\hat{f}}_{\mathrm{LOO},i}$ is the KDE constructed from all samples except X_i . Being linear in the logarithm of the estimate at observed values, this choice will penalize relative errors. Since we wish to obtain an accurate estimate of densities over a large dynamic range, the log likelihood is more suitable than considering absolute error or squared absolute error.

We then quantify counting uncertainties for the underlying inhomogeneous Poisson process using generalized bootstrap resampling (Chamandy et al. 2012): we take a number of PE samples from each detected event that is a random variable distributed as Poisson(1).

The new aspect of this work concerns the choice of mass sample values to treat as KDE input data X_i . We noted in Sadiq et al. (2022) that for nontrivial uncertainties in component masses this method, in addition to possible intrinsic biases due to the choice of a Gaussian KDE, will be biased toward an overdispersed estimate of the true distribution. This is expected because, with the uninformed priors used in PE, the posterior over mass for any given event will be randomly displaced relative to the (unknown) true mass due to detector noise. Here we motivate and present our strategy for correcting this overdispersion bias.

Our motivation is linked to Bayesian hierarchical population inference (Mandel 2010), where measurement errors are treated by considering the true event properties θ _i for detected events labeled i = 1, ..., N as nuisance parameters. Here, the likelihood of a set of gravitational-wave data segments d_i corresponding to the events is

$$\begin{eqnarray}&&{P}_{N}(\{d\}| {\boldsymbol{\lambda }})=\displaystyle \prod _{i=1}^{N}\int P({d}_{i}| {{\boldsymbol{\theta }}}_{i}){p}_{\mathrm{pop}}({{\boldsymbol{\theta }}}_{i}| {\boldsymbol{\lambda }})\,d{{\boldsymbol{\theta }}}_{i},\end{eqnarray} \tag{ 6 }$$

where P(d_i ∣...) is the likelihood for a single data segment and p_pop( θ _i ∣ λ ) is the population distribution over θ _i for population model hyperparameters λ (here, for simplicity, we omit selection effects). Inference is implemented using PE samples that were generated using a standard or fiducial prior π_PE( θ ), often chosen as uniform over parameters of interest (see, e.g., Veitch et al. 2015; Thrane & Talbot 2019). Samples (labeled by k) are distributed as the posterior density p( θ _i ∣...) using this prior, hence

$$\begin{eqnarray}&&{{\boldsymbol{\theta }}}_{i}^{k}\sim p({{\boldsymbol{\theta }}}_{i}| {d}_{i},{\pi }_{\mathrm{PE}}({\boldsymbol{\theta }}))\propto P({d}_{i}| {{\boldsymbol{\theta }}}_{i}){p}_{\mathrm{pop}}({{\boldsymbol{\theta }}}_{i}| {\boldsymbol{\lambda }})\displaystyle \frac{{\pi }_{\mathrm{PE}}({\boldsymbol{\theta }})}{{p}_{\mathrm{pop}}({{\boldsymbol{\theta }}}_{i}| {\boldsymbol{\lambda }})}.\end{eqnarray} \tag{ 7 }$$

Hence, integrals over θ _i , as in Equation (6) may be performed (up to a constant factor) by summing over samples ${{\boldsymbol{\theta }}}_{i}^{k}$ reweighted by the ratio of the population distribution to the PE prior, p_pop( θ _i ∣ λ )/π_PE( θ ).

Here, while not making use of this hierarchical likelihood, we remark that such reweighted PE sample sets give an unbiased estimate of event properties if p_pop( θ _i ) is equal to the true population distribution; conversely, the unweighted PE samples give a biased estimate if π_PE( θ ) differs from the true population distribution. This observation is the basic motivation for our improved method.

A KDE trained on points drawn randomly from PE samples will be biased because these samples are themselves biased by the "uninformative" prior. If we have access to an estimated population distribution ${\hat{p}}_{\mathrm{pop}}({\boldsymbol{\theta }})$ that is closer to the true distribution than the PE prior is, we will obtain more accurate estimates of event properties by drawing samples weighted proportional to ${\hat{p}}_{\mathrm{pop}}({\boldsymbol{\theta }})/{\pi }_{\mathrm{PE}}({\boldsymbol{\theta }})$, as detailed below. The better an estimate of the true distribution we are able to obtain, the smaller will be the bias in event parameters using reweighted PE samples, and ultimately the smaller will be the bias of the KDE.

The above discussion suggests an iterative procedure where, beginning with both biased PE samples and a biased population KDE, one may be improved in turn using the other, finally reaching a stationary state where—ideally—both the sample draws and the corresponding population estimates are unbiased. This iterative strategy is similar to the expectation-maximization algorithm (Dempster et al. 1977), a popular method to estimate parameters for statistical models when there are missing or incomplete data.

Our basic algorithm follows these steps:

• 1.  
  for each GW event, draw Poisson-distributed (with mean 1) PE samples weighted by the current estimate of population density ${\hat{p}}_{\mathrm{pop}}$;
• 2.  
  create an awKDE from this sample set, optimizing the global bandwidth (and sensitivity parameter α, if not fixed);
• 3.  
  update the current population estimate using one or more KDEs and the selection function, and go to step 1.

In more detail, in step 1 we draw PE samples with probability proportional to the ratio of ${\hat{p}}_{\mathrm{pop}}({{\boldsymbol{\theta }}}_{i}^{k})$ to the PE prior distribution. Step 2 reproduces our previous awKDE method. Step 3 relates the KDE of detected events to an estimate of the true population distribution, hence in general it requires us to compensate for the selection function over the event parameter space: i.e., we estimate the true distribution as the KDE of detected events divided by the probability of detection, as detailed in Sadiq et al. (2022, Section 3.1).

In step 3 we may choose to derive the updated population density ${\hat{p}}_{\mathrm{pop}}({\boldsymbol{\theta }})$ from only the most recently calculated KDE; then the iterative process is a Markov chain, ⁶ and we may characterize it via the autocorrelation of various scalar quantities computed at each iteration. We use the optimized global bandwidth h (and adaptive sensitivity parameter α, if not fixed to unity) for this purpose.

After discarding a small number of initial iterations and then accumulating a number significantly greater than the autocorrelation time, we expect the collection of iterations to provide unbiased (though not necessarily independent) estimates of the population distribution. For subsequent iterations we then use the median of ${\hat{p}}_{\mathrm{pop}}({{\boldsymbol{\theta }}}_{i}^{k})$ over a buffer of previous iterations (usually the previous 100) to determine the sample draw probabilities for the next iteration. This population estimate should be more precise than one using only a single previous KDE; in addition using the buffer estimate, the samples for each successive iteration are essentially independent.

In reality, we do not expect the resulting population density estimate to be entirely unbiased, due to more fundamental limitations in both the PE input (e.g., inaccuracies in the merger waveform model used, see Abbott et al. 2023b) and in the KDE procedure itself. Specifically, the expected bias of the KDE is proportional to the second derivative of the true distribution function (Silverman 1986, Section 3.3), thus, depending on the eventual choice of bandwidth, sharp peak or gap features may not be well represented. This bias might be reduced by replacing the KDE with a more general estimate, for instance Gaussian mixture models that allow the positions and weights of kernels to vary (e.g., Rinaldi & Del Pozzo 2022), although this implies higher complexity and computational cost.

We first test this iterative reweighting method on a simple mock data set. We generate true event parameters for a mixture, drawing 30 events each from a truncated power law p_pop(x) ∼ x^−0.5 and from a Gaussian with mean (standard deviation) of μ = 35 (σ_p = 3). We then simulate "measured" event values with a Gaussian random scatter relative to the true parameters of standard deviation σ_m = 5 (broader than the true Gaussian peak); for each measured value we generate 100 mock parameter samples with mean equal to the measured value and with the same uncertainty σ_m = 5. This procedure mimics PE using an uninformative (flat) prior over x.

First, applying awKDE as in Sadiq et al. (2022) to random draws from these mock parameter samples, we find as expected an overdispersed estimate around the peak (Figure 1, top). Second, applying awKDE with our iterative reweighting algorithm we obtain the bottom plot of Figure 1: here the Gaussian height and standard deviation appear accurately reconstructed and the true distribution is well within the 90% percentiles of iteration samples, except near the step-function truncations of the power law, which cannot be accurately represented by a Gaussian KDE.

Zoom In Zoom Out Reset image size

We verify that the initial Markov process has accumulated several independent samples by plotting the autocorrelation of the optimized global bandwidth h versus lag (separation along the chain) in Figure 2. ⁷ The autocorrelation drops near zero by a lag of ≲10 iterations, thus a buffer of 100 iterations contains several independent population estimates. (The estimate is more noisy for larger lags as fewer iterations are available.)

Zoom In Zoom Out Reset image size

To further investigate possible statistical biases, we performed the same two awKDE analyses (unweighted and iterative reweighted) for 50 realizations of mock data with the same distribution. At any given x value, we can find the percentile (p-value) of the true p_pop relative to the awKDE Monte Carlo samples, and thus produce a probability-probability (P–P) plot. The first panel in Figure 3 for unweighted awKDE shows consistent estimates for values well away from the population peak, but catastrophic under- or overestimation close to the peak, as expected from strong overdispersion. For reweighted awKDE, the second panel shows much smaller biases close to the peak. (We do not show x values close to 0 or 100, as these are subject to large edge effects for both methods.) While not entirely free of bias, the reweighted method offers much improved reconstruction of population features.

Zoom In Zoom Out Reset image size

In the next subsection we verify the improved behavior of the reweighted awKDE for two-dimensional mock data in the presence of correlated parameter uncertainties.

To investigate the performance of iterative reweighted KDE over 2D data, we start with 60 mock events, 50% from a two-dimensional uniform distribution and 50% from a bivariate normal (Gaussian) distribution. The mean of the normal distribution is μ = (50, 50) in arbitrary units and we take the two dimensions to be uncorrelated, each with a variance of 9. We then simulate measurement uncertainty in our mock data values with an anticorrelation between the two dimensions using covariance ${{\rm{\Sigma }}}_{p}=\left(\begin{array}{cc}32.5 & -31.5\\ -31.5 & 32.5\end{array}\right)$ , corresponding to an error ellipse with 8:1 axes at 45°. This is a simple choice to mimic the anticorrelation between m₁ and m₂ along a contour of constant chirp mass for lower-mass BBHs. As for the 1D mock data, 100 parameter samples with the same covariance are generated around each measured value.

We then investigate whether the reweighted awKDE can reconstruct the true (uncorrelated) population peak by reducing the correlation of the artificial measurement uncertainty. The top panel of Figure 4 for the unweighted awKDE (median over 900 iterations) clearly shows a biased reconstruction with anticorrelation around the peak at (50, 50). ⁸ The bottom panel shows the results using the iterative reweighting scheme: the anticorrelation appears to be removed and the reconstructed peak is as expected for the true distribution.

Zoom In Zoom Out Reset image size

We also compute 1D results from these 2D reweighted awKDE iteration samples by integrating over the second parameter (labeled m₂ in Figure 4) and compare with the corresponding 1D true distribution. Figure 5 shows this comparison: the true Gaussian peak height and width are well recovered by the reweighted awKDE, although as in the 1D case there are evident edge artifacts near the step-function truncations of the uniform distribution.

Zoom In Zoom Out Reset image size

We start by evaluating the effect of the iterative reweighting method on the 1D primary mass distribution, taking 100 random PE samples for each of the 69 BBH events; here, we assume a power-law distribution for secondary mass p(m₂) ∼ q^1.5. We reproduce the awKDE results from Sadiq et al. (2022) and use this estimate to seed the reweighted iteration algorithm. After ∼1000 reweighting iterations we compute the median and symmetric 90% interval from the last 900 rate estimates (the first 100 are used to set up a buffer for population weighting, as above): results are presented in Figure 6.

Zoom In Zoom Out Reset image size

Our estimate is generally consistent with other nonparametric or semiparametric approaches, represented by the Flexible mixtures and Power Law + Spline models in Abbott et al. (2023b), and does not show the overdispersion apparent in Figure 8 of Sadiq et al. (2022); specifically, we find a slightly higher and narrower peak around 35 M_⊙, but no identifiable feature around 20 M_⊙ (compare Tiwari 2022; Toubiana et al. 2023).

Next, we apply our reweighting scheme on PE samples for both component masses and compute the two-dimensional merger rate, using the estimated sensitive volume × time as a function of the two masses. We will first discuss various technical aspects of extending the 1D calculation without assuming any power-law dependence for m₂.

Binary exchange symmetry. Typically, the convention m₁ > m₂ is applied when presenting binary parameter estimates. However, all aspects of binary formation physics and event detection and parameter estimation will be invariant under swapping the component labels, i.e., exchanging m₁ ↔ m₂ (at the same time exchanging the spins). Thus, considering the differential merger rate ${ \mathcal R }({m}_{1},{m}_{2})$ as a function over the whole plane, it must have a reflection symmetry about the line m₁ = m₂. To respect this symmetry and remove biases resulting from the apparent lack of support at m₂ > m₁, if using samples with the m₁ > m₂ convention, we train and evaluate KDEs on reflected sample sets, i.e., after adding copies of the samples with swapped components. (Note also that a power-law m₂ distribution implies the rate is a nondifferentiable function at the equal-mass line, whereas a KDE by construction is smooth and differentiable everywhere.)

Choice of KDE parameters. In Sadiq et al. (2022), we mainly considered a KDE constructed over linear mass (or distance) parameters; however, here we choose the logarithms of component masses. While this choice is not expected to have a large impact on the results, since the kernel bandwidth is free to locally adapt in either case, it is technically preferable for a few reasons. We avoid any possible KDE support at negative masses; there is a generally higher density of events toward lower masses considering the entire 3–100 M_⊙ range; the density of observed events also shows less overall variation over log coordinates; and when evaluating the KDE on a grid with equal spacing, fewer points are required to maintain precision for the low-mass region.

For a 2D KDE we also have a choice of kernel parameters, i.e., the Gaussian covariance matrix: given the similar or identical physical interpretation and range of values between lnm₁ and lnm₂, we choose a covariance proportional to the unit matrix, with an overall factor determined by the local adaptive bandwidth for each event.

PE prior. The PE samples released by LVK use a prior uniform in component masses (LIGO Scientific Collaboration et al. 2021a) up to a factor dependent on cosmological redshift; we currently do not consider reweighting relative to the default cosmological model. As the prior is a density, it transforms with a Jacobian factor when changing variables to $\mathrm{ln}{m}_{1},\mathrm{ln}{m}_{2}$; thus, we must divide the estimated rate ${ \mathcal R }(\mathrm{ln}{m}_{1},\mathrm{ln}{m}_{2})$ by a prior ∝m₁ m₂ when obtaining reweighted draw probabilities.

With these technical choices, we perform 1500 reweighting iterations in total, the first 600 using the Markov chain (i.e., the immediately preceding rate estimate) for sample draw weights, and the remaining 900 using the buffer median estimate. Figure 7 shows the rate estimate computed with iterative reweighting for BBH events in GWTC-3 (median over the last 900 iterations). The autocorrelations of optimal global bandwidth and sensitivity parameter α for the first 600 iterations are shown in Figure 8: the correlation drops close to 0 at a lag of ∼30 iterations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The mass distribution shows several interesting features in addition to the expected peaks (overdensities) around primary masses of ∼10 M_⊙ and ∼35 M_⊙, with corresponding peaks over secondary mass. For primary masses of ∼35 M_⊙ up to 80 M_⊙, the most likely secondary mass is ∼30–35 M_⊙. Thus, over this range the two component masses appear almost independently chosen. Around the m₁ ∼ 10 M_⊙ peak, some anticorrelation of the two components appears, i.e., higher m₁ favors lower m₂. We investigate the significance of this feature by calculating the median of the m₂ distribution as a function of primary mass, i.e., ${m}_{2}^{50 \% }({m}_{1});$ this is a decreasing function over a range of m₁ from ∼10 through ∼14 M_⊙ for 90%–95% of our sample estimates, depending on the exact range of m₁ chosen. The anticorrelation is thus moderately but not highly significant.

Between the two peaks the distribution of mass ratios appears broader than at either one (as hinted at in Tiwari 2023), although the apparent trend is based on a small number of events. We also see a narrow lower-density region just above the ∼10 M_⊙ peak (cf. the local minimum at chirp mass ∼11 M_⊙ in Tiwari 2023).

We also integrated the 2D KDE rate estimate numerically over both m₁ and m₂ to obtain merger rates over component mass. As shown in Figure 9, we recover features consistent with the 2D estimates and with other methods. Each component mass distribution appears well modeled by a combination of two Gaussian peaks and (broken) power laws.

Zoom In Zoom Out Reset image size

To elucidate features in the 2D distribution, we choose various representative values of primary mass to plot the distribution of m₂ in Figure 10. The similarity between secondary distributions for m₁ ≳ 35 M_⊙ is evident.

Zoom In Zoom Out Reset image size

We may also derive the distribution of mass ratio q from our 2D rate estimate. We plot this for various representative values of primary mass in Figure 11, and compare to a typical power law ∝q^1.5. First, we see that the q distribution varies over primary mass; hence, models where it is forced to the same form over the whole mass range are likely to have nontrivial bias. For some primary masses, p(q) is consistent with a monotonically increasing function such as a positive power, but for others it clearly decreases over some of the range. Roughly, if m₁ is close to a peak then p(q) is consistent with an increasing power law, but for other values the mass ratio rather shows a maximum at intermediate values, down to q ∼ 0.5 for m₁ = 15 or m₁ = 70.

Zoom In Zoom Out Reset image size

This behavior may suggest that the primary and secondary masses are independently drawn from similar distributions, modulo a q-dependent "pairing factor" (Fishbach & Holz 2020) that influences the relative probability of binary merger (recent work by Farah et al. 2023b further explores the case of similar or identical underlying primary and secondary mass distributions). In any case, the preference toward q ∼ 1 seen in previous work is not confirmed here. Callister & Farr (2023) reached a similar conclusion, although assuming a "universal" distribution p(q) over all primary masses.

Results including GW190814. We also estimate the 2D and 1D integrated merger rates using our iterative reweighted KDE method when the outlier event GW190814 (Abbott et al. 2020), which has a mass ratio ${ \mathcal O }(10)$ and a secondary mass barely above the likely neutron star maximum mass, is included in the BBH population. Detailed results are presented in the Appendix: roughly summarizing the trends seen there, the bulk of the estimated distribution remains little changed by the addition of the extra event, although the peak in secondary mass below 10 M_⊙ is shifted toward lower values and becomes both higher and broader; this is likely due to a general increase in KDE bandwidth when optimized with cross-validation. (In parameterized models the estimated mass distribution is also highly sensitive to inclusion of GW190814 ; Abbott et al. 2021b, 2023b.) It is not clear whether other methods for choosing the bandwidth would yield more accurate estimates; higher event statistics in the regime of low m₂ are clearly desirable.