Compact binary coalescence parameter estimation
In this example, we demonstrate how to generate simulated data for a binary black hole coalescence observed by the two LIGO interferometers at Hanford and Livingston for all parameters in the IMRPhenomPv2 waveform model.
The code will take around 15 hours to run.
For testing, you may prefer to run the 4-parameter CBC tutorial.
1#!/usr/bin/env python
2"""
3Tutorial to demonstrate running parameter estimation on a full 15 parameter
4space for an injected cbc signal. This is the standard injection analysis script
5one can modify for the study of injected CBC events.
6
7This will take many hours to run.
8"""
9import bilby
10import numpy as np
11from bilby.core.utils.random import seed
12
13# Sets seed of bilby's generator "rng" to "123" to ensure reproducibility
14seed(123)
15
16# Set the duration and sampling frequency of the data segment that we're
17# going to inject the signal into
18duration = 4.0
19sampling_frequency = 1024.0
20
21# Specify the output directory and the name of the simulation.
22outdir = "outdir"
23label = "full_15_parameters"
24bilby.core.utils.setup_logger(outdir=outdir, label=label)
25
26
27# We are going to inject a binary black hole waveform. We first establish a
28# dictionary of parameters that includes all of the different waveform
29# parameters, including masses of the two black holes (mass_1, mass_2),
30# spins of both black holes (a, tilt, phi), etc.
31injection_parameters = dict(
32 mass_1=36.0,
33 mass_2=29.0,
34 a_1=0.4,
35 a_2=0.3,
36 tilt_1=0.5,
37 tilt_2=1.0,
38 phi_12=1.7,
39 phi_jl=0.3,
40 luminosity_distance=2000.0,
41 theta_jn=0.4,
42 psi=2.659,
43 phase=1.3,
44 geocent_time=1126259642.413,
45 ra=1.375,
46 dec=-1.2108,
47)
48
49# Fixed arguments passed into the source model
50waveform_arguments = dict(
51 waveform_approximant="IMRPhenomXPHM",
52 reference_frequency=50.0,
53 minimum_frequency=20.0,
54)
55
56# Create the waveform_generator using a LAL BinaryBlackHole source function
57# the generator will convert all the parameters
58waveform_generator = bilby.gw.WaveformGenerator(
59 duration=duration,
60 sampling_frequency=sampling_frequency,
61 frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole,
62 parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_black_hole_parameters,
63 waveform_arguments=waveform_arguments,
64)
65
66# Set up interferometers. In this case we'll use two interferometers
67# (LIGO-Hanford (H1), LIGO-Livingston (L1). These default to their design
68# sensitivity
69ifos = bilby.gw.detector.InterferometerList(["H1", "L1"])
70ifos.set_strain_data_from_power_spectral_densities(
71 sampling_frequency=sampling_frequency,
72 duration=duration,
73 start_time=injection_parameters["geocent_time"] - 2,
74)
75
76ifos.inject_signal(
77 waveform_generator=waveform_generator, parameters=injection_parameters
78)
79
80# For this analysis, we implement the standard BBH priors defined, except for
81# the definition of the time prior, which is defined as uniform about the
82# injected value.
83# We change the mass boundaries to be more targeted for the source we
84# injected.
85# We define priors in the time at the Hanford interferometer and two
86# parameters (zenith, azimuth) defining the sky position wrt the two
87# interferometers.
88priors = bilby.gw.prior.BBHPriorDict()
89
90time_delay = ifos[0].time_delay_from_geocenter(
91 injection_parameters["ra"],
92 injection_parameters["dec"],
93 injection_parameters["geocent_time"],
94)
95priors["H1_time"] = bilby.core.prior.Uniform(
96 minimum=injection_parameters["geocent_time"] + time_delay - 0.1,
97 maximum=injection_parameters["geocent_time"] + time_delay + 0.1,
98 name="H1_time",
99 latex_label="$t_H$",
100 unit="$s$",
101)
102del priors["ra"], priors["dec"]
103priors["zenith"] = bilby.core.prior.Sine(latex_label="$\\kappa$")
104priors["azimuth"] = bilby.core.prior.Uniform(
105 minimum=0, maximum=2 * np.pi, latex_label="$\\epsilon$", boundary="periodic"
106)
107
108# Initialise the likelihood by passing in the interferometer data (ifos) and
109# the waveoform generator, as well the priors.
110# The explicit distance marginalization is turned on to improve
111# convergence, and the posterior is recovered by the conversion function.
112likelihood = bilby.gw.GravitationalWaveTransient(
113 interferometers=ifos,
114 waveform_generator=waveform_generator,
115 priors=priors,
116 distance_marginalization=True,
117 phase_marginalization=False,
118 time_marginalization=False,
119 reference_frame="H1L1",
120 time_reference="H1",
121)
122
123# Run sampler. In this case we're going to use the `dynesty` sampler
124# Note that the `walks`, `nact`, and `maxmcmc` parameter are specified
125# to ensure sufficient convergence of the analysis.
126# We set `npool=4` to parallelize the analysis over four cores.
127# The conversion function will determine the distance posterior in post processing
128result = bilby.run_sampler(
129 likelihood=likelihood,
130 priors=priors,
131 sampler="dynesty",
132 nlive=1000,
133 walks=20,
134 nact=50,
135 maxmcmc=2000,
136 npool=4,
137 injection_parameters=injection_parameters,
138 outdir=outdir,
139 label=label,
140 conversion_function=bilby.gw.conversion.generate_all_bbh_parameters,
141 result_class=bilby.gw.result.CBCResult,
142)
143
144# Plot the inferred waveform superposed on the actual data.
145result.plot_waveform_posterior(n_samples=1000)
146
147# Make a corner plot.
148result.plot_corner()
Running this script will generate data then perform parameter estimation. In doing all of this, it prints information about the matched-filter SNRs in each detector (saved to the log-file). Moreover, it generates a plot for each detector showing the data, amplitude spectral density (ASD) and the signal; here is an example for the Hanford detector:
Finally, after running the parameter estimation. It generates a corner plot:
The solid lines indicate the injection parameters.