We will be using the PyCBC library, which is used to study gravitational-wave data, find astrophysical sources due to compact binary mergers, and study their parameters. These are some of the same tools that the LIGO and Virgo collaborations use to find gravitational waves in LIGO/Virgo data
In this tutorial we will walk through estimating the significance of a peak in a signal, given a simplified search (with some assumptions as noted along the way). We will also make use of one of the standard signal consistency tests to help remove some non-Gaussian transient noise from the background.
import sys
!{sys.executable} -m pip install pycbc ligo-common --no-cache-dir
Requirement already satisfied: pycbc in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (1.14.4) Requirement already satisfied: lalsuite in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (6.62) Requirement already satisfied: ligo-common in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (1.0.3) Requirement already satisfied: jinja2 in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from pycbc) (2.10) Requirement already satisfied: astropy>=2.0.3; python_version > "3.0" in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from pycbc) (3.0.4) Requirement already satisfied: beautifulsoup4>=4.6.0 in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from pycbc) (4.6.3) Requirement already satisfied: cython in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from pycbc) (0.29.14) Requirement already satisfied: mpld3>=0.3 in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from pycbc) (0.3) Requirement already satisfied: ligo-segments in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from pycbc) (1.2.0) Requirement already satisfied: matplotlib>=1.5.1 in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from pycbc) (3.0.0) Requirement already satisfied: Mako>=1.0.1 in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from pycbc) (1.1.0) Requirement already satisfied: six>=1.10.0 in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from pycbc) (1.11.0) Requirement already satisfied: requests>=1.2.1 in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from pycbc) (2.22.0) Requirement already satisfied: tqdm in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from pycbc) (4.26.0) Requirement already satisfied: scipy>=0.16.0; python_version >= "3.5" in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from pycbc) (1.1.0) Requirement already satisfied: h5py>=2.5 in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from pycbc) (2.8.0) Requirement already satisfied: emcee==2.2.1 in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from pycbc) (2.2.1) Requirement already satisfied: decorator>=3.4.2 in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from pycbc) (4.3.0) Requirement already satisfied: lscsoft-glue>=1.59.3 in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from pycbc) (2.0.0) Requirement already satisfied: numpy>=1.13.0; python_version > "3.0" in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from pycbc) (1.16.2) Requirement already satisfied: pillow in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from pycbc) (5.3.0) Requirement already satisfied: python-dateutil in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from lalsuite) (2.8.1) Requirement already satisfied: MarkupSafe>=0.23 in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from jinja2->pycbc) (1.1.0) Requirement already satisfied: cycler>=0.10 in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from matplotlib>=1.5.1->pycbc) (0.10.0) Requirement already satisfied: kiwisolver>=1.0.1 in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from matplotlib>=1.5.1->pycbc) (1.0.1) Requirement already satisfied: pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.1 in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from matplotlib>=1.5.1->pycbc) (2.3.0) Requirement already satisfied: idna<2.9,>=2.5 in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from requests>=1.2.1->pycbc) (2.7) Requirement already satisfied: certifi>=2017.4.17 in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from requests>=1.2.1->pycbc) (2018.10.15) Requirement already satisfied: chardet<3.1.0,>=3.0.2 in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from requests>=1.2.1->pycbc) (3.0.4) Requirement already satisfied: urllib3!=1.25.0,!=1.25.1,<1.26,>=1.21.1 in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from requests>=1.2.1->pycbc) (1.23) Requirement already satisfied: pyOpenSSL in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from lscsoft-glue>=1.59.3->pycbc) (18.0.0) Requirement already satisfied: setuptools in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from kiwisolver>=1.0.1->matplotlib>=1.5.1->pycbc) (41.6.0) Requirement already satisfied: cryptography>=2.2.1 in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from pyOpenSSL->lscsoft-glue>=1.59.3->pycbc) (2.3.1) Requirement already satisfied: asn1crypto>=0.21.0 in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from cryptography>=2.2.1->pyOpenSSL->lscsoft-glue>=1.59.3->pycbc) (0.24.0) Requirement already satisfied: cffi!=1.11.3,>=1.7 in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from cryptography>=2.2.1->pyOpenSSL->lscsoft-glue>=1.59.3->pycbc) (1.11.5) Requirement already satisfied: pycparser in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from cffi!=1.11.3,>=1.7->cryptography>=2.2.1->pyOpenSSL->lscsoft-glue>=1.59.3->pycbc) (2.19)
We will estimate the significance of signal-to-noise peak observed in the Virgo instrument near in time to the large peaks observed in the LIGO-Hanford and LIGO-Livingston observatories.
For this purpose we will consider the existence of a gravitational wave signal to have been confirmed based on the signal from the LIGO detectors alone, as was in fact the case for the matched-filtering based analyses of GW170814, as they did not incorporate any information from the Virgo data.
The question we will ask can be phrased as follows. What is the probability that noise can produce a peak as large or larger than the largest peak observed in the Virgo data, within a consistent lightspeed travel time between all three observatories? This is a form of null hypothesis testing, where we create a p-value.
For the purpose of this notebook, we may a few additional simplifying assumptions, and those will be stated as we go along.
In this section, we will read in a short segment of data round GW170814, and do some basic preconditioning, as also demonstrated in previous tutorials. We will also calculate the power spectrum of the data.
Notably, we are making the assumption here that the power spectral estimate of the data is constant over this short stretch of time, and isn't biased by the fact that we chose to center the estimate (very roughly) on the event time. We do not assume that the data is actually stationary, Gaussian, or is free from non-astrophysical transient artefacts.
%matplotlib inline
import pylab
from pycbc.filter import resample_to_delta_t, highpass
from pycbc.catalog import Merger
from pycbc.psd import interpolate, inverse_spectrum_truncation
m = Merger("GW170814")
ifos = ['H1', 'L1', 'V1']
data = {}
psd = {}
pylab.figure(figsize=[10, 5])
for ifo in ifos:
# Read in and precondition the data
ts = m.strain(ifo).highpass_fir(15, 512)
data[ifo] = resample_to_delta_t(ts, 1.0/2048).crop(2, 2)
# Estimate the power spectral density of the data
# This chooses to use 4s samples in the PSD estimate.
# One should note that the tradeoff in segment length is that
# resolving narrow lines becomes more difficult.
p = data[ifo].psd(4)
p = interpolate(p, data[ifo].delta_f)
p = inverse_spectrum_truncation(p, int(4 * data[ifo].sample_rate), low_frequency_cutoff=20.0)
psd[ifo] = p
pylab.plot(psd[ifo].sample_frequencies, psd[ifo], label=ifo)
pylab.yscale('log')
pylab.xscale('log')
pylab.ylim(1e-47, 1e-41)
pylab.xlim(20, 1024)
pylab.ylabel('$Strain^2 / Hz$')
pylab.xlabel('Frequency (Hz)')
pylab.grid()
pylab.legend()
pylab.show()
To calculate the signal-to-noise time series, we need to generate an estimate of the signal. For this purpose we will assume the source black holes are non spinning, have equal mass, and the agree with the total mass estimate for the system as a whole. A better method would be to use the maximum likelihood estimate from an analysis of the LIGO data alone, however, this is sufficient for the purposes of this tutorial.
from pycbc.waveform import get_fd_waveform
from pycbc.filter import matched_filter
from pycbc.conversions import mass1_from_mchirp_q
# Calculate the component mass of each black hole in the detector frame
mchirp= m.median1d("mchirp") # This is in the source frame
zfac = (1 + m.median1d("redshift")) # apply redshift to get to the detector frame
mchirp *= zfac
mass1 = mass2 = mass1_from_mchirp_q(mchirp, 1)
# This is a frequency domain waveform generator. It has a very similar syntax to the time domain
# waveform function used in prior tutorials. This function returns both a plus and a cross
# polarization waveform, but we will just use the plus polarization in building our template
# as these are only different by a phase offset in this specific case.
hp, _ = get_fd_waveform(approximant="IMRPhenomD",
mass1=mass1, mass2=mass2,
f_lower=20.0, delta_f=data[ifo].delta_f)
hp.resize(len(psd[ifo]))
# For each observatory use this template to calculate the SNR time series
snr = {}
for ifo in ifos:
snr[ifo] = matched_filter(hp, data[ifo], psd=psd[ifo], low_frequency_cutoff=25)
snr[ifo] = snr[ifo].crop(5, 4)
# Show a couple sizes
for w, title in [(7, 'Wide View'), (.15, 'Close to GW170814')]:
pylab.figure(figsize=[14, 4])
for ifo in ifos:
pylab.plot(snr[ifo].sample_times, abs(snr[ifo]), label=ifo)
pylab.legend()
pylab.title(title)
pylab.grid()
pylab.xlim(m.time - w, m.time + w)
pylab.ylim(0, 15)
pylab.xlabel('Time (s)')
pylab.ylabel('Signal-to-noise (SNR)')
pylab.show()
We see in the SNR time series plots above, that while there are nice peaks around GW170814 in each detector, there are also some large peaks at other times. LIGO / Virgo data, does contain transient (i.e limited duration) noise artefacts that an analyses must deal with to search LIGO data with high sensitivity. One approach for dealing with this is outlined later in this tutorial.
One of the ways we can test how well the data actual fits the models to use a $\chi^2$-based signal consistency test. We employ a version of the test described in this paper. Schematically, we chop up our template into $p$ number of bins and see how much each contributes to the SNR ($\rho_i$). We can then calculate our statistic as the difference between the SNR in one bin, and the expected fraction of the total SNR ($\rho$).
$ \chi^2 = \sum^p_{i=0} (\rho_i - \rho / p)^2 $
This will have $2p-2$ degrees of freedom as each SNR is complex representing both possible orthogonal phases the signal could have contributions from. There is also a constraint due to the fact that the sum of each bin must each the total SNR by definition. In this notebook we will normalize this statistic by dividing by the number of degrees of freedom, producing $\chi^2_r$.
We expect that this statistic will be high when the template does not match well the data, and near unity when the data either is Gaussian noise, or it contains the expected signal in addition to Gaussian noise.
from pycbc.vetoes import power_chisq
chisq = {}
for ifo in ifos:
# The number of bins to use. In principle, this choice is arbitrary. In practice,
# this is empirically tuned.
nbins = 26
chisq[ifo] = power_chisq(hp, data[ifo], nbins, psd[ifo], low_frequency_cutoff=20.0)
chisq[ifo] = chisq[ifo].crop(5, 4)
dof = nbins * 2 - 2
chisq[ifo] /= dof
pylab.figure(figsize=[14, 4])
for ifo in ifos:
pylab.plot(chisq[ifo].sample_times, chisq[ifo], label=ifo)
pylab.legend()
pylab.grid()
pylab.xlim(m.time -0.15, m.time + 0.15)
pylab.ylim(0, 5)
pylab.ylabel('$chi^2_r$')
pylab.show()
There are some notable features in the $\chi^2_r$ time series. We see that there is a dip in the value at the time of the peak in the SNR in each observatory. We expect this as the template now aligns with the signal in the data. Also, the values climb just around this minima. This occurs because the template is starting to slide against the true signal in the data but is not perfectly aligned with it.
One approach we can take is to down-weight the times where the data does not appear as either Guassian noise or Gaussian noise + our template. We can do this be combining the SNR time series and our $\chi^2_r$ time series as follows. This is a method used to re-weight SNR since initial LIGO, and has been employed in the first two Advanced LIGO observing runs. In this tutorial we will choose to rank our events by this statistic.
$\hat{\rho} = \frac{\rho}{ [1 + (\chi^2_r)^3]^{1/6}}$ where $\chi^2 > 1$, otherwise $\rho$
For reference on how we rank coincident (i.e. occuring in multiple detector) events in Advanced LIGO, there is a description here.
from pycbc.events.ranking import newsnr
# The rho-hat term above is named "newsnr" here
nsnr = {ifo:newsnr(abs(snr[ifo]), chisq[ifo]) for ifo in ifos}
# Show a couple sizes
for w, title in [(7, 'Wide View'), (.15, 'Close to GW170814')]:
pylab.figure(figsize=[14, 4])
for ifo in ifos:
pylab.plot(snr[ifo].sample_times, nsnr[ifo], label=ifo)
pylab.legend()
pylab.title(title)
pylab.grid()
pylab.xlim(m.time - w, m.time + w)
pylab.ylim(0, 15)
pylab.xlabel('Time (s)')
pylab.ylabel('Re-weighted Signal-to-noise')
pylab.show()
We can see above that there are still peaks around GW170814 in all detectors, at roughly the same signal strength, but that at other times, where that had been peaks in the time series, there are no longer large statistic values.
In this section we will determine how significant the peak in the virgo re-weighted SNR time series is.
We will do this first by determining where one might expect a peak relative to the LIGO observed peaks. This is set by the constraint that an astrophysical source can only cause delays between observatories no larger than the light travel time between them. The pycbc.detector.Detector
class provides some convenient methods to ask these sorts of questions.
We will then calculate the peak in the SNR for this window around the LIGO observed peaks. This is our "on-source".
Finally, to determine the significance of the on-source we will compare how likely it is for a peak as large or larger to appear in the background. Our background will be empirically measured by taking portions of the SNR time series from the "off-source" i.e. times that do not overlap the on-source. An important criteria to avoid a biased significance estimate is that the background and experiment be performed in the same manner.
import numpy
from pycbc.detector import Detector
# Calculate the time of flight between the Virgo detectors and each LIGO observatory
d = Detector("V1")
tof = {}
tof['H1'] = d.light_travel_time_to_detector(Detector("H1"))
tof['L1'] = d.light_travel_time_to_detector(Detector("L1"))
# Record the time of the peak in the LIGO observatories
ptime = {}
pylab.figure(figsize=[14, 4])
for ifo in ifos:
# shade the region around each LIGO peak that could have a peak in Virgo if from
# an astrophysical source
if ifo is not 'V1':
ptime[ifo] = snr[ifo].sample_times[nsnr[ifo].argmax()]
pylab.axvspan(ptime[ifo] - tof[ifo], ptime[ifo] + tof[ifo], alpha=0.2, lw=10)
pylab.plot(snr[ifo].sample_times, nsnr[ifo], label=ifo)
# Calculate the span of time that a Virgo peak could in principle happen in from time of flight
# considerations.
start = ptime['H1'] - tof['H1']
end = ptime['L1'] + tof['L1']
# convert the times to indices along with how large the region is in number of samples
window_size = int((end - start) * snr['V1'].sample_rate)
sidx = int((start - snr['V1'].start_time) * snr['V1'].sample_rate)
eidx = sidx + window_size
# Calculate the "on-source" peak re-weighted (newsnr) statistic value.
onsource = nsnr['V1'][sidx:eidx].max()
pylab.legend()
pylab.grid()
pylab.xlim(m.time - .05, m.time + .10)
pylab.ylim(0, 15)
pylab.xlabel('Time (s)')
pylab.ylabel('Re-weighted Signal-to-noise')
pylab.show()
print('Virgo Peak has a statistic value of {}'.format(onsource))
Virgo Peak has a statistic value of 4.121405829558103
In the plot above we see the re-weighted SNR time series. On top of that we've shaded the regions which are consistent with a Virgo signal based on the peaks in the LIGO observatories. Only in the darker region, is it possible to have a peak in the SNR that is consistent with both LIGO observatories.
# Now that we've calculate the onsource peak, we should calculate the background peak values.
# We do this by chopping up the time series into chunks that are the same size as our
# onsource and repeating the same peak finding (max) procedure.
# Walk through the data in chunks and calculate the peak statistic value in each.
peaks = []
i = 0
while i + window_size < len(nsnr['V1']):
p = nsnr['V1'][i:i+window_size].max()
peaks.append(p)
i += window_size
# Skip past the onsource time
if abs(i - sidx) < window_size:
i += window_size * 2
peaks = numpy.array(peaks)
# The p-value is just the number of samples observed in the background with a
# value equal or higher than the onsource divided by the number of samples.
# We can make the mapping between statistic value and p-value using our background
# samples.
pcurve = numpy.arange(1, len(peaks)+1)[::-1] / float(len(peaks))
peaks.sort()
pvalue = (peaks > onsource).sum() / float(len(peaks))
pylab.figure(figsize=[10, 7])
pylab.scatter(peaks, pcurve, label='Off-source (Noise Background)', color='black')
pylab.axvline(onsource, label='On-source', color='red')
pylab.axhline(pvalue, color='red')
pylab.legend()
pylab.yscale('log')
pylab.grid()
pylab.ylim(1e-3, 1e0)
pylab.ylabel('p-value')
pylab.xlabel('Re-weighted Signal-to-noise')
pylab.xlim(2, 5)
pylab.show()
print("The p-value associated with the GW170814 peak is {}".format(pvalue))
The p-value associated with the GW170814 peak is 0.007228915662650603
In this tutorial, we find a peak in Virgo as large as the obseved one has an approximately 1% chance of occuring due to the noise alone. Given the simplifications of this tutorial, we find a result in agreement with the GW170814 discovery paper which reported a p-value of 0.3%.