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

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.