Shifted SNR distribution when realizing an event with Gaussian noise multiple times
When running a matched-filter search on an event multiple times with Gaussian noise, we expect to see each interferometer's SNR distributed as a Gaussian about the optimal SNR with unit variance (see fig. 3 in [1]). For a network of N (in this case, four) detectors, the resulting network SNR^2 should be distributed as a noncentral chi-square distribution with N degrees of freedom and the non-centrality parameter equal to the network optimal SNR^2. (see [2])
However, this isn't the distribution that I'm seeing, for 100 matched-filter searches of the same event using bayestar-realize-coincs
:
Here is code to replicate this example. The exact files can also be found on the CIT cluster at /home/geoffrey.mo/bayestar-test/noise_snr
. My ligo.skymap version is 0.1.16.
lalapps_inspinj -o inj.xml --m-distr fixMasses --fixed-mass1 1.3 --fixed-mass2 1.3 --t-distr uniform --time-step 100 --gps-start-time 1000000000 --gps-end-time 1000000099 --d-distr volume --min-distance 150000 --max-distance 160000 --l-distr random --i-distr uniform --f-lower 20 --disable-spin --waveform IMRPhenom_NRTidal
# Make a copy of this, call it inj_many.xml.
cp inj.xml inj_many.xml
# DO THIS MANUALLY #
# Then copy the line in the sim_inspiral inj_many.xml table many times (in this case, 100).
# Make sure to add the delimiter at the end of the line before copying.
bayestar-sample-model-psd -o psd.xml --H1=aLIGOAdVO4T1800545 --L1=aLIGOAdVO4T1800545 --V1=AdVMidLowSensitivityP1200087 --K1=KAGRAEarlySensitivityT1600593
# Get the zero-noise (optimal) SNR:
bayestar-realize-coincs -o coinc.xml inj.xml --reference-psd psd.xml --detector H1 L1 V1 K1 --snr-threshold 0.0 --waveform IMRPhenomD_NRTidal --min-triggers 1 --net-snr-threshold 1 -j 64
# Now localize the many injs:
bayestar-realize-coincs -o coinc_many.xml inj_many.xml --reference-psd psd.xml --detector H1 L1 V1 K1 --snr-threshold 0.0 --waveform IMRPhenomD_NRTidal --min-triggers 1 --net-snr-threshold 1 -j 64 --measurement-error gaussian-noise
# Run this python script to parse the xmls for the SNRs and plot the distribution
python parse_plot.py
And this is parse_plot.py:
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
import xml.etree.ElementTree as ET
# These don't seem to work due to some ligolw incompatabilities with the
# coinc_event_id column validity. I really don't want to debug this right now
# so I'm just going to do some ugly parsing of the xml to get the SNRs.
# coinc_one = EventTable.read('coinc.xml', tablename='coinc_inspiral')
# coinc_many = EventTable.read('coinc_many.xml', tablename='coinc_inspiral')
single_tree = ET.parse('coinc.xml')
single_root = single_tree.getroot()
optimal_snr = float(single_root[1][-1].text.split(',')[-1].split('\n')[0])
many_tree = ET.parse('coinc_many.xml')
many_root = many_tree.getroot()
events = [event.split(',') for event in many_root[1][-1].text.split('\n')[:-1] if
len(event) > 1]
events[-1].append('') # last row needs fixing to look like the other rows
snrs = np.array([float(event[-2]) for event in events])
xarray = np.linspace(5, 15, 200)
plt.figure(figsize=(12,8))
# this is the distribution we expect for the quadrature sum
# of four Gaussians (one for each detector)
plt.hist(np.random.noncentral_chisquare(4, optimal_snr**2, 100000),
bins=200, histtype='step', density=True, label='Noncentral chi squared')
plt.hist(snrs**2, bins=20, histtype='step', density=True, label='Observed');
plt.axvline(optimal_snr**2, ls='--', color='black', label=r'$\rho_\mathrm{opt}$');
plt.legend(loc='upper left');
plt.xlabel('SNR$^2$');
plt.ylabel('Probability density');
plt.savefig('snr2_dist.png', dpi=300)
[1] https://arxiv.org/pdf/1809.02293.pdf
[2] Mathematica code:
In[1] = TransformedDistribution[
s^2 + t^2 + u^2 + v^2, {s \[Distributed] NormalDistribution[w, 1],
t \[Distributed] NormalDistribution[x, 1],
u \[Distributed] NormalDistribution[y, 1],
v \[Distributed] NormalDistribution[z, 1]}]
Out[1] = NoncentralChiSquareDistribution[4, w^2 + x^2 + y^2 + z^2]