Analyzing simulated posteriors

A standard test of population inference is to analyze a simulated population to verify that we can recover the true population. The optimal way to do this is to run full parameter estimation on signals injected into real interferometer noise. However, to reduce the computational cost, gwpopulation_pipe implements a simple method for generating Gaussian posteriors based on a specified population model.

We can perform such an analysis with the following bash script.

#!/usr/bin/env bash

RUNDIR=simulation
LOGDIR=$RUNDIR/logs
LABEL=simluation

RESULT_FILE=$RUNDIR/result/${LABEL}_result.json
SAMPLE_FILE=$RUNDIR/result/${LABEL}_samples.hdf5

mkdir $RUNDIR

gwpopulation_pipe_collection empty.txt\
  --run-dir $RUNDIR --log-dir $LOGDIR --label $LABEL --data-label $LABEL\
  --parameters mass_1 --parameters mass_ratio --parameters a_1 --parameters a_2 --parameters redshift --parameters cos_tilt_1 --parameters cos_tilt_2\
  --mass-models b --magnitude-models iid --tilt-models iid\
  --n-simulations 10 --samples-per-posterior 1000 --injection-file samples.json --injection-index 0


gwpopulation_pipe_analysis empty.txt\
  --run-dir $RUNDIR --log-dir $LOGDIR --label $LABEL --data-label $LABEL\
  --models gwpopulation.models.mass.power_law_primary_mass_ratio --models iid_spin\
  --prior-file test.prior --sampler pymultinest --sampler-kwargs "{nlive: 100}"\
  --enforce-minimum-neffective-per-event True\
   --injection-file samples.json --injection-index 0\
  --vt-function ""

gwpopulation_pipe_plot empty.txt --run-dir $RUNDIR --result-file $RESULT_FILE --samples $SAMPLE_FILE --redshift False

gwpopulation_pipe_to_common_format -r $RESULT_FILE -s $SAMPLE_FILE --models gwpopulation.models.mass.power_law_primary_mass_ratio --models iid_spin

Note

In this example, we use the pymultinest sampler rather than the default dynesty. This does not install automatically with gwpopulation_pipe, however it can be simply installed with

$ conda install -c conda-forge pymultinest

pymultinest can’t be simply installed using pypi as it relies on the fortran package MultiNest.

The file samples.json contains the true values of the hyperparameters to simulate the universe for. This should be in a format that can be read in with pandas.read_json.

{
  "alpha":{"0": 3},
  "beta": {"0": 1},
  "mmin": {"0": 6},
  "mmax": {"0": 50},
  "alpha_chi": {"0": 2},
  "beta_chi": {"0": 2},
  "amax": {"0": 1},
  "xi_spin": {"0": 0.5},
  "sigma_spin": {"0": 1}
}

This will create a directory simulation with two subdirectories data and result. The data directory contains the simulated posterior and violin plots of the input parameters. The result directory contains the output of the sampler and various plots produced in post-processing.