Population Fisher forecasts

Whereas gwforge_fisher forecasts how precisely a detector network can measure the parameters of an individual source, gwforge_population_fisher forecasts how precisely it can constrain the hyper-parameters describing an astrophysical population, such as the mass function, merger-rate history, and spin distribution.

The implementation computes Term I of the population hyper-parameter Fisher matrix derived by Gair, Antonelli & Barbieri (2022), Eq. 21. For a catalogue of \(N_{\rm det}\) detected events with parameters \(\theta_k\), and population hyper-parameters \(\Lambda\), the score for event \(k\) is

\[ s_k^i = \frac{\partial \ln p(\theta_k \mid \Lambda)} {\partial \Lambda^i}, \]

and the Fisher matrix, together with its inverse covariance matrix, is

\[ \Gamma^{ij} = \sum_{k=1}^{N_{\rm det}} \big(s_k^i-\langle s\rangle^i\big) \big(s_k^j-\langle s\rangle^j\big), \qquad \Sigma = \Gamma^{-1}. \]

Here, \(\langle s\rangle^i\) denotes the sample mean of the scores over the detected catalogue. The selection function is not evaluated explicitly: the detected catalogue itself provides the Monte Carlo sample used to estimate the Fisher matrix. Consequently, the SNR threshold is an input to the forecast, not a separate selection-correction step. For a fixed population and detection threshold, the forecast uncertainties scale as \(T_{\rm obs}^{-1/2}\), where \(T_{\rm obs}\) is the observing time.

Configuration

A representative configuration for a joint mass, redshift, and spin forecast is shown below.

[Catalogue]
population-file = bbh_population.h5
snr-file = bbh_snr.h5
detectors = ['CE40', 'CE20', 'ET']
snr-thresholds = [10, 20]

[Model]
blocks = ['mass', 'redshift', 'spin']
maximum-redshift = 10
cosmology = Planck18
mmin = 4.489078237
m-high = 300.0
maximum-mass = 300.0
mass-parameters = {'alpha_1': 1.456442737, 'alpha_2': 5.100400428, 'm_break': 37.9806382, 'lam_0': 0.4206178692, 'lam_1': 0.5228817702, 'mpp_1': 9.989384122, 'sigpp_1': 0.6601839103, 'mpp_2': 33.2656028, 'sigpp_2': 4.582221363, 'delta_m': 3.123416382, 'beta': 0.8049660633}
redshift-parameters = {'gamma': 2.7, 'kappa': 5.6, 'z_peak': 1.9}
spin-parameters = {'mu_chi': 0.0751318233, 'sigma_chi': 0.3667469243, 'mu_t': 0.2788662234, 'sigma_t': 0.9661688072, 'xi_spin': 0.6650168221}

[Fit]
fit-redshift = True
fit-mass = True
fit-spin = True

[Fisher]
free-parameters = ['alpha_1', 'alpha_2', 'm_break', 'lam_0', 'lam_1', 'mpp_1', 'sigpp_1', 'mpp_2', 'sigpp_2', 'beta', 'gamma', 'kappa', 'z_peak', 'mu_chi', 'sigma_chi', 'mu_t', 'sigma_t', 'xi_spin']
score-method = analytic

[Output]
output-directory = population_fisher_output
label = mass_redshift_spin
plot-corner = True
plot-parameters = ['beta', 'gamma', 'kappa', 'z_peak', 'mu_chi']

[Catalogue]: population and detector network

population-file specifies the catalogue generated by gwforge_population, while snr-file points to the corresponding SNRs computed by gwforge_optimal_snr. The detectors entry specifies the detector network used to calculate the network SNR, and snr-thresholds defines the detection thresholds to consider. A separate forecast is produced for each threshold.

[Model]: population model and fiducial parameters

The blocks entry selects the population components to include: mass, redshift, and spin. For a detector-frame spectral-siren forecast, include the cosmology block to constrain \(H_0\), \(\Omega_{\rm m}\), and \(w_0\).

The population densities are those implemented in GWForge.population, from which gwforge_population draws its injections. Thus, the forecast model and the generated catalogue describe the same population by construction. The parameter names in the configuration are the names used by the population-generation code.

Fiducial parameters can be supplied explicitly through mass-parameters, redshift-parameters, and spin-parameters. Any omitted fiducials fall back to the GWTC-5.0 median values in GWForge.population.mass.BGP_PARAMETERS or the defaults in GWForge.population.spin.DEFAULT_BBH_SPIN_PARAMETERS, as appropriate.

The parameters mmin, m-high, and maximum-mass define hard mass cutoffs rather than population hyper-parameters. Their centred scores vanish identically, so these cutoffs are fixed by construction and must match those used to generate the catalogue.

Note: sigma_chi denotes a standard deviation. The Beta spin models instead use a variance, named sigma_squared_chi. Configurations using the old name are rejected rather than silently reinterpreted.

[Fit]: fitting the fiducial population

The fit-redshift, fit-mass, and fit-spin flags control the population-fit checks.

GWForge.population.redshift.Redshift models the merger-rate history by convolving the Madau–Dickinson star-formation rate with a time-delay distribution. Accordingly, fit-redshift = True fits the direct parameters of this model to the population, rather than fitting the distribution of detected events. This distinction is essential: fitting the detected distribution directly would absorb selection effects into the inferred population hyper-parameters.

The mass and spin fits serve as consistency checks rather than corrections to the fiducial model. They should recover the parameter values used to generate the catalogue; both the generation and fitted values are reported in the summary.

[Fisher]: free parameters and score calculation

free-parameters specifies the hyper-parameters included in the Fisher matrix. Any parameter omitted from this list is held fixed at its configured fiducial value. Fixing parameters can be particularly useful when the likelihood or population model contains sharp features or discontinuities that complicate local derivative estimates.

In particular, all shipped configurations hold delta_m fixed. This parameter controls the low-mass taper, a sharp feature at a fixed source-frame mass. Since it also couples strongly to \(H_0\), its influence should be explored by rerunning the forecast with different fixed values rather than including it in the standard Fisher matrix.

score-method = analytic selects the production implementation. The finite-difference option is also available and can be helpful where I or Claudio haven’t implemented the analytical forms.

[Output]: results and plots

output-directory sets the destination for the results, and label supplies the common filename prefix. If plot-corner = True, a corner plot is generated for the parameters listed in plot-parameters. This allows the full Fisher matrix to include many hyper-parameters while restricting the corner plot to a smaller, more interpretable subset.

Running a forecast

First generate the population catalogue and calculate its network SNRs. Then run the population Fisher forecast:

gwforge_population \
    --config-file GWForge/population/population_configuration_files/bgp-gwtc5.ini \
    --output-file bbh_population.h5 \
    --source-type bbh \
    --reference-frequency 5 \
    --seed 250114 \
    --save-config

gwforge_optimal_snr \
    --injection-file bbh_population.h5 \
    --output-file bbh_snr.h5 \
    --ifos CE40 CE20 ET \
    --waveform-approximant IMRPhenomXPHM \
    --minimum-frequency 5 \
    --sampling-frequency 2048 \
    --cores 10

gwforge_population_fisher \
    --config-file mass_redshift_spin.ini

The first command generates the population using the specified configuration and saves the generation configuration alongside the catalogue. The second computes the optimal network SNRs. The final command evaluates the population Fisher forecast using mass_redshift_spin.ini.

Ready-made configurations are available in GWForge/population_fisher/configuration_files/:

  • redshift.ini: redshift-only forecast.

  • mass.ini: mass-only forecast.

  • mass_redshift.ini: joint mass and redshift forecast.

  • mass_redshift_spin.ini: joint mass, redshift, and spin forecast.

  • spectral_siren.ini: detector-frame spectral-siren forecast.

  • mass_lvk.ini: mass forecast for a second-generation detector network.

  • spectral_siren_aplus.ini: spectral-siren forecast for an A+ network.

Outputs and comparison between detector networks

For each SNR threshold, the forecast is saved as a pickle, {label}_snr{threshold}.pkl, in the configured output directory. It holds a dictionary with the Fisher matrix, covariance matrix, parameter names, fiducial values, and detected-event counts.

The run also produces:

  • A corner plot for the parameters specified by plot-parameters.

  • Reconstructions of the primary-mass distribution \(\pi(m_1)\) and redshift distribution \(p(z)\), with 90% uncertainty bands evaluated over the injections.

  • A {label}_summary.txt file containing the fit and forecast tables. The summary is also printed to standard output.

Previously saved forecasts can be compared without recomputing the Fisher matrices, provided the forecasts share the same free parameters. For example:

gwforge_population_fisher \
    --compare "CE40+CE20+ET=population_fisher_output/mass_snr10.pkl" \
              "H1+L1+V1=population_fisher_output/mass_lvk_snr10.pkl" \
    --parameters "['alpha_1', 'alpha_2', 'm_break']" \
    --output-file mass_xg_vs_lvk.pdf

This compares the selected mass hyper-parameters for the XG and LVK networks at an SNR threshold of 10 and saves the comparison to mass_xg_vs_lvk.pdf.

Limitations and interpretation: Term I corresponds to the small-measurement-error limit, in which the uncertainty on each event’s parameters is much narrower than the underlying population distribution. It can therefore be optimistic near sharp features of the mass function. Terms II–V of the full population Fisher matrix are not yet implemented. Moreover, an 18-parameter forecast can contain strongly correlated directions, so marginal uncertainties should be interpreted alongside the correlation structure. Use PopulationFisherResult.correlation_matrix() to inspect these correlations before quoting marginal \(\sigma\) values.