Sampling the full hierarchical population posterior distribution in gravitational-wave astronomy

TL;DR

The paper addresses the challenge of GW population inference by performing full hierarchical sampling of the population hyperparameters and all event parameters, yielding a high-dimensional posterior with about $500$ dimensions for the GWTC-3 catalog. It employs probabilistic programming (e.g., PyMC) and Hamiltonian Monte Carlo, using a continuous Gaussian Mixture interpolation of transformed single-event posteriors to enable joint sampling of $oldsymbol{ ho}$ and $ heta_i$ while accounting for selection via the function $\xi(oldsymbol{ ho})$. The authors demonstrate that this approach recovers standard population constraints, provides population-informed event posteriors, and reveals direct event-by-event contributions to population features through correlation metrics $ ho_{ijk}$ and $oldsymbol{oldsymbol{}}_{ik}$, applied to 69 BH events. This method improves accuracy by removing certain Monte Carlo integrations, offers richer physical insights, and is scalable to larger catalogs and future cosmological extensions, with open-source tooling and potential for emulator-based selection functions and dark-siren cosmology.

Abstract

We present a full sampling of the hierarchical population posterior distribution of merging black holes using current gravitational-wave data. We directly tackle the most relevant intrinsic parameter space made of the binary parameters (masses, spin magnitudes, spin directions, redshift) of all the events entering the GWTC-3 LIGO/Virgo/KAGRA catalog, as well as the hyperparameters of the underlying population of sources. This results in a parameter space of about 500 dimensions, in contrast with current investigations where the targeted dimensionality is drastically reduced by marginalizing over all single-event parameters. In particular, we have direct access to (i) population parameters, (ii) population-informed single-event parameters, and (iii) correlations between these two sets of parameters. We quantify the fractional contribution of each event to the constraints on the population hyperparameters. Our implementation relies on modern probabilistic programming languages and Hamiltonian Monte Carlo, with a continuous interpolation of single-event posterior probabilities. Sampling the full hierarchical problem is feasible, as demonstrated here, and advantageous as it removes some (but not all) of the Monte Carlo integrations that enter the likelihood together with the related variances.

Figures (5)

• Figure 1: Posterior distribution of the population hyperparameters, assuming a fixed cosmology. The dark blue distribution refers to our full hierarchical sampling of Eq. (\ref{['h_post']}). The light blue distribution reports results from Ref. KAGRA:2021duu, which are restricted to the marginalized posterior of Eq. (\ref{['h_post_marginal']}). Contours correspond to $68\%$ and $90\%$ credible intervals.
• Figure 2: Single-event marginal posterior distributions for some representative events and some representative parameters. We show results from our full inference run (dark blue, filled), parameter-estimation results obtained with uninformative priors (light blue, filled), and population-informed reweighted results from Ref. KAGRA:2021duu (black, empty)
• Figure 3: Fractional contribution from each event in the GWTC-3 catalog to the population hyperparameters constraints. This is quantified using the averaged correlation coefficient of Eq. (\ref{['omegaik']}). Rows indicate the different hyperameters, either fixing (left) or varying (right) the cosmology. Vertical black lines indicate the fractional contribution of a given event in chronological order. Colors refer to the LVK data-taking periods. For each hyperparameter, the event contributing the most is indicated explicitly.
• Figure 4: Examples of correlations between single-event parameters and population/cosmological parameters. The top panel shows the primary mass of GW190521, the upper mass cutoff, and the Hubble constant. The middle panel shows the primary mass of $\rm GW200224 \_222234$, the location of the Gaussian peak in the mass spectrum, and the Hubble constant. The bottom panel shows the primary spin of $\rm GW190727\_06033$, and the moments of the population spin distribution. Dark blue (light blue) distributions refer to runs where we fix (vary) the cosmological parameters. Contours correspond to $68\%$ and $95\%$ credible intervals.
• Figure 5: Comparison of inference results with different cuts on the log--likelihood variance for a subset of the population parameters. The blue distribution corresponds to the choice made in the rest of the paper, where we threshold on the likelihood variance; the purple and light blue distributions explore variations around that threshold value; the teal distribution instead considers a threshold on the effective number of events. Contours correspond to $68\%$ and $90\%$ credible intervals.