Hierarchical Population Likelihoods

Scrollable formula. Use Left and Right arrow keys to reveal hidden content.

Scrollable table. Use Left and Right arrow keys to reveal additional columns.

Scrollable figure. Use Left and Right arrow keys to inspect hidden content.

A population analysis must combine the full likelihood from every event, not a histogram of posterior medians. Sources are modeled as independent random occurrences, distributed unevenly over their physical parameters; the search retains only some of them. This is an inhomogeneous, selection-thinned Poisson process. Its likelihood retains event uncertainty, divides out the interim priors used for event-level sampling, and accounts for the probability that a source drawn from the proposed population would enter the catalog.

Required background. Probability and Statistical Inference supplies conditional probability, marginalization, and Poisson processes. Helpful background. The Gravitational-Wave Likelihood defines the event likelihood; Selection Functions and Injection Campaigns explains how detection efficiency is measured.

A catalog as a thinned Poisson process

Let $\theta$ denote source parameters in a declared measure, $\Lambda$ population-shape hyperparameters, and $\nu$ the expected number of sources before detection over a stated observing interval and finite source domain. Hyperparameters describe a distribution rather than one event: a typical mass and its population spread are examples. The expected source count $\overline N_{\rm src}$ has density

\frac{d\overline N_{\mathrm{src}}}{d\theta} =\nu\,\pi(\theta\mid\Lambda), \qquad \int d\theta\,\pi(\theta\mid\Lambda)=1,

The normalized density $\pi$ describes the relative distribution of sources; the dimensionless number $\nu$ sets their total abundance. This keeps the statistical derivation separate from the conversion of a physical rate density into a count. The latter requires comoving volume and source-frame time dilation and is derived in Merger Rates and Sensitive Spacetime Volume.

If $p_{\rm det}(\theta)$ is the probability that the entire search-and-catalog procedure selects a source, then

\alpha(\Lambda) =\int d\theta\,p_{\rm det}(\theta)\pi(\theta\mid\Lambda)

is the population-averaged detection efficiency. The expected detected count is $\mu=\nu\alpha(\Lambda)$ . For data $d_i$ from $N$ selected events with negligible terrestrial contamination, the point-process likelihood is

\mathcal L(\{d_i\},N\mid \nu,\Lambda) \propto e^{-\nu\alpha(\Lambda)}\nu^N \prod_{i=1}^{N} \int d\theta_i\, p(d_i\mid\theta_i)\pi(\theta_i\mid\Lambda).

To see each factor, define $\mathcal I_i(\Lambda)=\int p(d_i\mid\theta)\pi(\theta\mid\Lambda)d\theta$ . For data lying inside the detection region, the density of one detected datum is $\mathcal I_i/\alpha$ . Multiplying these $N$ densities by the Poisson count probability gives

\frac{e^{-\nu\alpha}(\nu\alpha)^N}{N!} \prod_i\frac{\mathcal I_i}{\alpha} =\frac{e^{-\nu\alpha}\nu^N}{N!}\prod_i\mathcal I_i.

The data-ordering convention changes only the constant $N!$ . There is no additional factor of $p_{\rm det}(\theta_i)$ inside an event integral when $d_i$ contains the data used for selection. Selection has already acted through the allowed data and Poisson normalization. Adding it again double counts detectability. If selection uses extra information omitted from $d_i$ , that information must instead be included or marginalized consistently. This construction follows Mandel, Farr, and Gair 2019, §§2–4.

Conditioning on the observed count directly gives $\prod_i\mathcal I_i/\alpha^N$ . The same dependence follows for $N\ge1$ by marginalizing $\nu$ with the scale prior $p(\nu)\propto1/\nu$ , because $\int_0^\infty e^{-\alpha\nu}\nu^{N-1}d\nu=\Gamma(N)\alpha^{-N}$ :

\mathcal L(\{d_i\}\mid\Lambda,N) \propto \alpha(\Lambda)^{-N} \prod_i \int d\theta_i\, p(d_i\mid\theta_i)\pi(\theta_i\mid\Lambda).

This conditional-on-count form is useful for population shape, but it discards count-based rate information by construction. A different abundance prior gives a different marginal factor. For example, a prior uniform in $\nu$ gives $\alpha^{-(N+1)}$ . The improper scale prior is also unsuitable for a zero-event rate analysis; one must not copy $\alpha^{-N}$ without stating the conditioning and measure.

Reusing event-level posterior samples

An interim prior is the prior originally used to analyze an event, before fitting the population. Suppose event $i$ used $\pi_0(\theta)$ and has evidence $Z_i$ , so that

p(\theta\mid d_i,\pi_0) =\frac{p(d_i\mid\theta)\pi_0(\theta)}{Z_i}.

Let $\theta_{ij}$ be sample $j$ among the $n_i$ posterior samples retained for event $i$ . Multiplying the posterior density by the target-to-interim prior ratio replaces the event prior with the proposed population:

\int d\theta\,p(d_i\mid\theta)\pi(\theta\mid\Lambda) =Z_i\, \mathbb E_{p(\theta\mid d_i,\pi_0)} \!\left[\frac{\pi(\theta\mid\Lambda)}{\pi_0(\theta)}\right] \approx \frac{Z_i}{n_i}\sum_{j=1}^{n_i} \frac{\pi(\theta_{ij}\mid\Lambda)}{\pi_0(\theta_{ij})}.

Here $\mathbb E$ denotes an expectation over the event posterior, approximated by the sample average on the right. The $Z_i$ factors cancel in hyperparameter ratios only if the same event likelihood and interim analysis are used throughout. The ratio also requires support: wherever $\pi(\theta\mid\Lambda)\gt0$ , the interim prior must be nonzero. A hypermodel that puts weight outside the event prior cannot be repaired by more posterior samples.

For a concrete two-class example, take an interim prior $(0.5,0.5)$ and an event posterior $(0.8,0.2)$ for classes A and B. A population with fractions $(f,1-f)$ gives

\frac{\mathcal I_i(f)}{Z_i} =0.8\frac{f}{0.5}+0.2\frac{1-f}{0.5} =0.4+1.2f.

This smoothly averages both possible classes; assigning the event entirely to A would discard its uncertainty. The population analysis must still apply the selection and count treatment derived above.

A useful Monte Carlo warning is the effective sample size

n_{{\rm eff},i} =\frac{\left(\sum_j w_{ij}\right)^2}{\sum_j w_{ij}^2}, \qquad w_{ij}=\frac{\pi(\theta_{ij}\mid\Lambda)}{\pi_0(\theta_{ij})}.

Large $n_i$ is not enough when a few weights dominate. Hyperposterior features that appear where event or injection $n_{\rm eff}$ collapses are numerical warnings, not discoveries. Thrane and Talbot give posterior recycling in §5 and posterior prediction in Appendix D; their Appendix F erratum corrects the original discussion of single-event selection Thrane and Talbot 2019, corrected PDF v9. An ordinary source-parameter posterior does not acquire an extra inverse-detection-efficiency weight.

Selection by injections

An injection campaign draws $\theta_k\sim p_{\rm inj}(\theta)$ , embeds signals in representative data, and records $I_k=1$ when the declared pipeline and threshold recover them. Importance sampling estimates

\alpha(\Lambda) \approx \frac{1}{N_{\rm inj}} \sum_{k=1}^{N_{\rm inj}} I_k\frac{\pi(\theta_k\mid\Lambda)}{p_{\rm inj}(\theta_k)}.

The injection distribution must cover the hypermodel, including mass, spin, sky position, inclination, distance or redshift, and observing-state dependence. Search pipelines, thresholds, detector segments, waveform families, calibration, and cosmology are part of this estimator. An analytic signal-to-noise cutoff may teach the geometry, but it does not replace a catalog-specific injection campaign.

The current worked interface is the official GWTC-5.0 event/posterior products together with its corresponding search-sensitivity and injection products, frozen at the evidence cutoff 2026-08-18. Any numerical population conclusion therefore inherits that release, event-inclusion rule, $p_{\rm astro}$ treatment, calibration and waveform assumptions, and injection support LVK 2026, GWTC-5.0 population analysis, §§2–3 and appendices.

Contamination and uncertain membership

Near threshold, treating every trigger as certainly astrophysical biases both rate and shape. A mixture point process instead gives each candidate an intensity sum,

\mathcal L\propto e^{-(\mu_{\rm s}+\mu_{\rm n})} \prod_i\left[\mu_{\rm s}\,\ell_{{\rm s},i}(\Lambda) +\mu_{\rm n}\,\ell_{{\rm n},i}(\eta)\right],

Here $\mu_{\rm s}$ and $\mu_{\rm n}$ are expected selected counts, and $\ell_{{\rm s},i}$ and $\ell_{{\rm n},i}$ are data densities normalized over the selected-data domain; in the signal-only model, $\ell_{{\rm s},i}=\mathcal I_i/\alpha$ . The parameter $\eta$ describes terrestrial-background statistics. Without these normalizations the intensity expression would be ambiguous. A published $p_{\rm astro}$ is not generally a fixed, population-independent weight: it was computed under a particular signal population, noise model, and pipeline. Either model the mixture coherently or state the approximation used to import those probabilities.

Checks that should precede interpretation

Selected mock recovery. Simulate a Poisson population, pass it through the actual selection rule, and verify recovery of abundance $\nu$ and shape $\Lambda$ . Omitting $\alpha$ should visibly reproduce the Malmquist bias.
Interim-prior invariance. Reanalyze or reweight the same controlled events under two sufficiently broad interim priors. The hyperposterior should agree within Monte Carlo error.
Injection support. Map injection $n_{\rm eff}$ over the hyperposterior. Truncate claims where the selection integral is extrapolated.
Misspecification stress test. Generate from a distribution outside the fitted family. Posterior predictive discrepancies should reveal, rather than conceal, the failure.
Nuisance propagation. Vary waveform, calibration, cosmology, catalog threshold, and outlier treatment when they are capable of moving the claimed feature.

These tests distinguish a precisely sampled hyperposterior from a trustworthy population statement. A clean sampler cannot cure an incorrect catalog likelihood.

Exercises

1. Why posterior medians fail. Two events have broad, overlapping likelihoods in a scalar parameter $x$ . Explain why fitting a population width to their posterior medians is biased.

Solution

Each median collapses measurement uncertainty to a point and retains the interim prior’s influence. The dispersion of medians is therefore neither the intrinsic population dispersion nor the distribution of future measurements. The hierarchical integral convolves the proposed population with each event likelihood and divides out the interim prior; it can attribute observed spread to measurement error rather than forcing all spread into the population.

2. A constant-efficiency limit. If $p_{\rm det}(\theta)=p_0$ for all parameters and shape inference is conditioned on $N$ , show that selection cannot change the inferred shape.