Population likelihood

Gravitational-wave transient surveys provide a biased sample of the astrophysical population. The likelihood function used for population inference is given be

\[ \begin{align}\begin{aligned}{\cal L}(\{d_i\} | \Lambda) &= \prod_i {\cal L}(d_i | \Lambda, {\rm det})\\&= \prod_i \frac{{\cal L}(d_i | \Lambda)}{P_{\rm det}(\Lambda)}\\&= \frac{1}{P_{\rm det}^{N}(\Lambda)} \prod_i \int d\theta_i p(d_i | \theta_i) \pi(\theta_i | \Lambda).\end{aligned}\end{align} \]

The quantity \(P_{\rm det}(\Lambda)\) is the detection probability for a single source (see selection.html).

The integrals over the per-event parameters \(\theta_i\) are typically performed using Monte Carlo integration

\[\hat{{\cal L}}(d_i | \Lambda) = \frac{1}{K} \sum_{k=1}^K \frac{\pi(\theta_k | \Lambda)}{\pi(\theta_k | \varnothing)}.\]

The full approximate log-likelihood is then given by

\[\ln \hat{\cal L}(\{d_i\} | \Lambda) = \sum_i \ln \hat{{\cal L}}(d_i | \Lambda) - N \ln \hat{P}_{\rm det}(\Lambda).\]

This approximation is implemented in gwpopulation.hyperpe.HyperparameterLikelihood.

There is another related expression for the likelihood as the result of an inhomoegeneous Poisson process. In this case the likelihood is given by

\[\ln {\cal L}(\{d_i\} | \Lambda) = N \ln R - N_{\rm exp}(\Lambda) + \sum_i \ln \hat{{\cal L}}(d_i | \Lambda)\]

Here \(R\) is the total merger rate and \(T\) is the total observation time and \(N_{\rm exp}(\Lambda) = RT\hat{P}_{\rm det}(\Lambda)\). This is implemented in gwpopulation.hyperpe.RateLikelihood.

Each of these Monte Carlo integrals have associated uncertainties which are propagated through the likelihood calculation and can be calculated using gwpopulation.hyperpe.HyperparameterLikelihood.ln_likelihood_and_variance().