Extracting distribution parameters from multiple uncertain observations with selection biases

We derive a Bayesian framework for incorporating selection effects into population analyses. We allow for both measurement uncertainty in individual measurements and, crucially, for selection biases on the population of measurements, and show how to extract the parameters of the underlying distribution based on a set of observations sampled from this distribution. We illustrate the performance of this framework with an example from gravitational-wave astrophysics, demonstrating that the mass ratio distribution of merging compact-object binaries can be extracted from Malmquist-biased observations with substantial measurement uncertainty.

💡 Research Summary

The paper presents a rigorous Bayesian hierarchical framework that simultaneously accounts for measurement uncertainties and selection effects when inferring the underlying population distribution of astrophysical sources. The authors begin by defining individual event parameters (θ) and population hyper‑parameters (λ), separating λ into an overall rate (N) and shape parameters (λ₀) that describe the functional form of the distribution (e.g., means and covariances of a multivariate Gaussian or histogram bin weights).

A “bottom‑up” derivation is employed, starting from the simplest case of perfect measurements and no selection bias, and progressively adding complexity. In the idealized scenario, the likelihood for a set of observed events is simply the product of the population density evaluated at each true θ_i, divided by the integral of the density over the full parameter space (which equals one if the density is properly normalized).

Selection bias is introduced via a detection probability p_det(θ), which is expressed as an integral over data space of the probability that the data exceed a detection threshold (the indicator function I(d)). This formulation captures both intrinsic selection effects (e.g., larger chirp masses yielding higher signal‑to‑noise ratios) and stochastic fluctuations due to instrumental noise. The detection probability appears only in the normalization factor α(λ₀)=∫p_det(θ)p_pop(θ|λ₀)dθ, representing the fraction of the total astrophysical population that would be observed for a given set of hyper‑parameters.

When measurement uncertainties are present, the observed data d_i are linked to the true parameters through a likelihood p(d_i|θ_i). By invoking Bayes’ theorem for each single‑event analysis, p(d_i|θ_i)=p(θ_i|d_i)p(d_i)/π(θ), the authors rewrite the full hierarchical likelihood in terms of posterior samples {θ_i^j} obtained from standard parameter‑estimation pipelines. The integral over θ in the likelihood is replaced by a discrete sum over these samples, yielding a computationally tractable expression:

p({d_i}|λ₀) = ∏{i=1}^{N_obs} (1/S_i) Σ{j=1}^{S_i} p_pop(θ_i^j|λ₀) / α(λ₀)

where S_i is the number of posterior samples for event i. This expression automatically incorporates both measurement error (through the spread of posterior samples) and selection bias (through α).

The framework is extended to include the overall event rate N. The expected number of detections is N_det = N·α(λ₀), and the probability of observing N_obs events follows a Poisson distribution e^{-N_det} N_det^{N_obs}/N_obs!. By assigning a non‑informative prior π(N)∝1/N, the rate can be analytically marginalized, leaving the same posterior for λ₀ as in the rate‑free case, up to a constant that does not affect inference.

To demonstrate practicality, two examples are provided. The first revisits a classic flux‑limited survey, showing how the luminosity function can be recovered when only the brightest sources are detected. The second, more astrophysically relevant, concerns the mass‑ratio distribution of merging binary neutron stars observed via gravitational waves. Here, the detection probability scales as the 5/6 power of the chirp mass (p_det∝M_c^{5/6}), reflecting the volume increase for higher‑mass systems. Simulated third‑generation detector data with a signal‑to‑noise threshold of 20 indicate that roughly one thousand detections are required to constrain the underlying mass‑ratio distribution with useful precision.

In summary, the authors deliver a complete mathematical derivation, clarify common misconceptions (such as incorrectly inserting the detection probability directly into the per‑event likelihood without proper normalization), and supply a ready‑to‑use algorithmic recipe. The methodology is broadly applicable to any astronomical field where selection effects and measurement errors coexist, including galaxy surveys, X‑ray cluster catalogs, and future gravitational‑wave population studies.