...
 
Commits (5)
# These are the default priors for analysing GW150914.
mass_1 = Uniform(name='mass_1', minimum=30, maximum=50, unit='$M_{\odot}$', boundary=None)
mass_2 = Uniform(name='mass_2', minimum=20, maximum=40, unit='$M_{\odot}$', boundary=None)
mass_ratio = Constraint(name='mass_ratio', minimum=0.125, maximum=1)
a_1 = Uniform(name='a_1', minimum=0, maximum=0.8, boundary='reflective')
a_2 = Uniform(name='a_2', minimum=0, maximum=0.8, boundary='reflective')
tilt_1 = Sine(name='tilt_1', boundary='reflective')
tilt_2 = Sine(name='tilt_2', boundary='reflective')
phi_12 = Uniform(name='phi_12', minimum=0, maximum=2 * np.pi, boundary='periodic')
phi_jl = Uniform(name='phi_jl', minimum=0, maximum=2 * np.pi, boundary='periodic')
luminosity_distance = bilby.gw.prior.UniformSourceFrame(name='luminosity_distance', minimum=1e2, maximum=1e3, unit='Mpc', boundary=None)
dec = Cosine(name='dec', boundary='reflective')
ra = Uniform(name='ra', minimum=0, maximum=2 * np.pi, boundary='periodic')
theta_jn = Sine(name='theta_jn', boundary='reflective')
psi = Uniform(name='psi', minimum=0, maximum=np.pi, boundary='periodic')
phase = Uniform(name='phase', minimum=0, maximum=2 * np.pi, boundary='periodic')
geocent_time = Uniform(1126259462.322, 1126259462.522, name='geocent_time', unit='$s$', boundary=None)
# These are the calibration parameters as described in
# https://journals.aps.org/prx/abstract/10.1103/PhysRevX.6.041015
# recalib_H1_frequency_0 = 20
# recalib_H1_frequency_1 = 54
# recalib_H1_frequency_2 = 143
# recalib_H1_frequency_3 = 383
# recalib_H1_frequency_4 = 1024
# recalib_H1_amplitude_0 = Gaussian(mu=0, sigma=0.048, name='recalib_H1_amplitude_0), '$\\delta A_{H0}$'
# recalib_H1_amplitude_1 = Gaussian(mu=0, sigma=0.048, name='recalib_H1_amplitude_1), '$\\delta A_{H1}$'
# recalib_H1_amplitude_2 = Gaussian(mu=0, sigma=0.048, name='recalib_H1_amplitude_2), '$\\delta A_{H2}$'
# recalib_H1_amplitude_3 = Gaussian(mu=0, sigma=0.048, name='recalib_H1_amplitude_3), '$\\delta A_{H3}$'
# recalib_H1_amplitude_4 = Gaussian(mu=0, sigma=0.048, name='recalib_H1_amplitude_4), '$\\delta A_{H4}$'
# recalib_H1_phase_0 = Gaussian(mu=0, sigma=0.056, name='recalib_H1_phase_0', '$\\delta \\phi_{H0}$')
# recalib_H1_phase_1 = Gaussian(mu=0, sigma=0.056, name='recalib_H1_phase_1', '$\\delta \\phi_{H1}$')
# recalib_H1_phase_2 = Gaussian(mu=0, sigma=0.056, name='recalib_H1_phase_2', '$\\delta \\phi_{H2}$')
# recalib_H1_phase_3 = Gaussian(mu=0, sigma=0.056, name='recalib_H1_phase_3', '$\\delta \\phi_{H3}$')
# recalib_H1_phase_4 = Gaussian(mu=0, sigma=0.056, name='recalib_H1_phase_4', '$\\delta \\phi_{H4}$')
# recalib_L1_frequency_0 = 20
# recalib_L1_frequency_1 = 54
# recalib_L1_frequency_2 = 143
# recalib_L1_frequency_3 = 383
# recalib_L1_frequency_4 = 1024
# recalib_L1_amplitude_0 = Gaussian(mu=0, sigma=0.082, name='recalib_L1_amplitude_0), '$\\delta A_{L0}$'
# recalib_L1_amplitude_1 = Gaussian(mu=0, sigma=0.082, name='recalib_L1_amplitude_1), '$\\delta A_{L1}$'
# recalib_L1_amplitude_2 = Gaussian(mu=0, sigma=0.082, name='recalib_L1_amplitude_2), '$\\delta A_{L2}$'
# recalib_L1_amplitude_3 = Gaussian(mu=0, sigma=0.082, name='recalib_L1_amplitude_3), '$\\delta A_{L3}$'
# recalib_L1_amplitude_4 = Gaussian(mu=0, sigma=0.082, name='recalib_L1_amplitude_4), '$\\delta A_{L4}$'
# recalib_L1_phase_0 = Gaussian(mu=0, sigma=0.073, name='recalib_L1_phase_0', '$\\delta \\phi_{L0}$')
# recalib_L1_phase_1 = Gaussian(mu=0, sigma=0.073, name='recalib_L1_phase_1', '$\\delta \\phi_{L1}$')
# recalib_L1_phase_2 = Gaussian(mu=0, sigma=0.073, name='recalib_L1_phase_2', '$\\delta \\phi_{L2}$')
# recalib_L1_phase_3 = Gaussian(mu=0, sigma=0.073, name='recalib_L1_phase_3', '$\\delta \\phi_{L3}$')
# recalib_L1_phase_4 = Gaussian(mu=0, sigma=0.073, name='recalib_L1_phase_4', '$\\delta \\phi_{L4}$')
......@@ -2,32 +2,32 @@
Examples
========
1. `General inference examples <https://git.ligo.org/lscsoft/bilby/tree/master/examples/other_examples>`_
1. `General inference examples <https://git.ligo.org/lscsoft/bilby/tree/master/examples/core_examples>`_
- `A simple Gaussian likelihood <https://git.ligo.org/lscsoft/bilby/blob/master/examples/core_examples/gaussian_example.py>`_: a good example to see how to write your own likelihood.
- `Linear regression for unknown noise <https://git.ligo.org/lscsoft/bilby/blob/master/examples/core_examples/linear_regression_unknown_noise.py>`_: fitting to general time-domain data.
- `A simple Gaussian likelihood <https://git.ligo.org/lscsoft/bilby/blob/master/examples/other_examples/gaussian_example.py>`_: a good example to see how to write your own likelihood.
- `Linear regression for unknown noise <https://git.ligo.org/lscsoft/bilby/blob/master/examples/other_examples/linear_regression_unknown_noise.py>`_: fitting to general time-domain data.
2. `Examples of injecting and recovering
data <https://git.ligo.org/lscsoft/bilby/tree/master/examples/injection_examples>`__
data <https://git.ligo.org/lscsoft/bilby/tree/master/examples/gw_examples/injection_examples>`__
- `4-parameter CBC
tutorial <https://git.ligo.org/lscsoft/bilby/blob/master/examples/injection_examples/fast_tutorial.py>`__
tutorial <https://git.ligo.org/lscsoft/bilby/blob/master/examples/gw_examples/injection_examples/fast_tutorial.py>`__
- `15-parameter CBC tutorial
<https://git.ligo.org/lscsoft/bilby/blob/master/examples/injection_examples/fast_tutorial.py>`__
<https://git.ligo.org/lscsoft/bilby/blob/master/examples/gw_examples/injection_examples/standard_15d_cbc_tutorial.py>`__
- `Create your own source
model <https://git.ligo.org/lscsoft/bilby/blob/master/examples/injection_examples/create_your_own_source_model.py>`__
model <https://git.ligo.org/lscsoft/bilby/blob/master/examples/gw_examples/injection_examples/create_your_own_source_model.py>`__
- `Create your own time-domain source
model <https://git.ligo.org/lscsoft/bilby/blob/master/examples/injection_examples/create_your_own_time_domain_source_model.py>`__
model <https://git.ligo.org/lscsoft/bilby/blob/master/examples/gw_examples/injection_examples/create_your_own_time_domain_source_model.py>`__
- `How to specify the
prior <https://git.ligo.org/lscsoft/bilby/blob/master/examples/injection_examples/how_to_specify_the_prior.py>`__
prior <https://git.ligo.org/lscsoft/bilby/blob/master/examples/gw_examples/injection_examples/how_to_specify_the_prior.py>`__
- `Using a partially marginalized
likelihood <https://git.ligo.org/lscsoft/bilby/blob/master/examples/injection_examples/marginalized_likelihood.py>`__
likelihood <https://git.ligo.org/lscsoft/bilby/blob/master/examples/gw_examples/injection_examples/marginalized_likelihood.py>`__
3. `Examples using open
data <https://git.ligo.org/lscsoft/bilby/tree/master/examples/open_data_examples>`__
3. `Examples using open
data <https://git.ligo.org/lscsoft/bilby/tree/master/examples/gw_examples/data_examples>`__
- `Analysing the first Binary Black hole detection,
GW150914 <https://git.ligo.org/lscsoft/bilby/blob/master/examples/open_data_examples/GW150914.py>`__
GW150914 <https://git.ligo.org/lscsoft/bilby/blob/master/examples/gw_examples/data_examples/GW150914.py>`__
4. `Notebook-style
tutorials <https://git.ligo.org/lscsoft/bilby/tree/master/examples/tutorials>`__
......
========
Examples
========
Here we host a large number of example scripts separated into core_examples,
which demonstrate the core features of :code:`bilby`, and gw_examples, which
demonstrate the gravitational-wave-specific features of :code:`bilby`. Finally,
there is also a directory containing various tutorials in the form of
notebooks.
The :code:`basic_tutorial.py` example discussed in arXiv:1811.02042 has been
renamed to `gw_examples/injection_examples/fast_tutorial.py
<https://git.ligo.org/lscsoft/bilby/tree/master/examples/gw_examples/injection_examples/fast_tutorial.py>`_.
mass_ratio = Uniform(name='mass_ratio', minimum=0.125, maximum=1, boundary='reflective')
chirp_mass = Uniform(name='chirp_mass', minimum=25, maximum=31, unit='$M_{\odot}$', boundary='reflective')
mass_1 = Constraint(name='mass_1', minimum=10, maximum=80)
mass_2 = Constraint(name='mass_2', minimum=10, maximum=80)
a_1 = Uniform(name='a_1', minimum=0, maximum=0.99, boundary='reflective')
a_2 = Uniform(name='a_2', minimum=0, maximum=0.99, boundary='reflective')
tilt_1 = Sine(name='tilt_1', boundary='reflective')
tilt_2 = Sine(name='tilt_2', boundary='reflective')
phi_12 = Uniform(name='phi_12', minimum=0, maximum=2 * np.pi, boundary='periodic')
phi_jl = Uniform(name='phi_jl', minimum=0, maximum=2 * np.pi, boundary='periodic')
luminosity_distance = PowerLaw(alpha=2, name='luminosity_distance', minimum=50, maximum=2000, unit='Mpc', latex_label='$d_L$')
dec = Cosine(name='dec', boundary='reflective')
ra = Uniform(name='ra', minimum=0, maximum=2 * np.pi, boundary='periodic')
theta_jn = Sine(name='theta_jn', boundary='reflective')
psi = Uniform(name='psi', minimum=0, maximum=np.pi, boundary='periodic')
phase = Uniform(name='phase', minimum=0, maximum=2 * np.pi, boundary='periodic')
#!/usr/bin/env python
"""
Tutorial to demonstrate running parameter estimation on GW150914 using open
data.
Tutorial to demonstrate running parameter estimation on GW150914
This example estimates all 15 parameters of the binary black hole system using
commonly used prior distributions. This will take a few hours to run.
commonly used prior distributions. This will take several hours to run. The
data is obtained using gwpy, see [1] for information on how to access data on
the LIGO Data Grid instead.
[1] https://gwpy.github.io/docs/stable/timeseries/remote-access.html
"""
from __future__ import division, print_function
import bilby
from gwpy.timeseries import TimeSeries
logger = bilby.core.utils.logger
outdir = 'outdir'
label = 'GW150914'
time_of_event = bilby.gw.utils.get_event_time(label)
# This sets up logging output to understand what bilby is doing
bilby.core.utils.setup_logger(outdir=outdir, label=label)
# Here we import the detector data. This step downloads data from the
# LIGO/Virgo open data archives. The data is saved to an `Interferometer`
# object (here called `H1` and `L1`). A Power Spectral Density (PSD) estimate
# is also generated and saved to the same object. To check the data and PSD
# makes sense, for each detector a plot is created in the `outdir` called
# H1_frequency_domain_data.png and LI_frequency_domain_data.png. The two
# objects are then placed into a list.
# For GW170817, 170608 and 170814 add the following line to select clean data.
# kwargs = {"tag": 'CLN'}
interferometers = bilby.gw.detector.get_event_data(label)
# Data set up
trigger_time = 1126259462
roll_off = 0.4 # Roll off duration of tukey window in seconds, default is 0.4s
duration = 4 # Analysis segment duration
post_trigger_duration = 2 # Time between trigger time and end of segment
end_time = trigger_time + post_trigger_duration
start_time = end_time - duration
psd_duration = 32 * duration
psd_start_time = start_time - psd_duration
psd_end_time = start_time
# We now use gwpy to obtain analysis and psd data and create the ifo_list
ifo_list = bilby.gw.detector.InterferometerList([])
for det in ["H1", "L1"]:
logger.info("Downloading analysis data for ifo {}".format(det))
ifo = bilby.gw.detector.get_empty_interferometer(det)
data = TimeSeries.fetch_open_data(det, start_time, end_time)
ifo.strain_data.set_from_gwpy_timeseries(data)
logger.info("Downloading psd data for ifo {}".format(det))
psd_data = TimeSeries.fetch_open_data(det, psd_start_time, psd_end_time)
psd_alpha = 2 * roll_off / duration
psd = psd_data.psd(
fftlength=duration,
overlap=0,
window=("tukey", psd_alpha),
method="median"
)
ifo.power_spectral_density = bilby.gw.detector.PowerSpectralDensity(
frequency_array=psd.frequencies.value, psd_array=psd.value)
ifo_list.append(ifo)
logger.info("Saving data plots to {}".format(outdir))
ifo_list.plot_data(outdir=outdir, label=label)
# We now define the prior.
# We have defined our prior distribution in a file packaged with BILBY.
# We have defined our prior distribution in a local file, GW150914.prior
# The prior is printed to the terminal at run-time.
# You can overwrite this using the syntax below in the file,
# or choose a fixed value by just providing a float value as the prior.
prior = bilby.gw.prior.BBHPriorDict(filename='GW150914.prior')
priors = bilby.gw.prior.BBHPriorDict(filename='GW150914.prior')
# In this step we define a `waveform_generator`. This is out object which
# In this step we define a `waveform_generator`. This is the object which
# creates the frequency-domain strain. In this instance, we are using the
# `lal_binary_black_hole model` source model. We also pass other parameters:
# the waveform approximant and reference frequency.
# the waveform approximant and reference frequency and a parameter conversion
# which allows us to sample in chirp mass and ratio rather than component mass
waveform_generator = bilby.gw.WaveformGenerator(
frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole,
parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_black_hole_parameters,
waveform_arguments={'waveform_approximant': 'IMRPhenomPv2',
'reference_frequency': 50})
# In this step, we define the likelihood. Here we use the standard likelihood
# function, passing it the data and the waveform generator.
likelihood = bilby.gw.likelihood.GravitationalWaveTransient(
interferometers, waveform_generator)
ifo_list, waveform_generator, priors=priors, time_marginalization=True,
phase_marginalization=True, distance_marginalization=True)
# Finally, we run the sampler. This function takes the likelihood and prior
# along with some options for how to do the sampling and how to save the data
result = bilby.run_sampler(likelihood, prior, sampler='dynesty',
outdir=outdir, label=label)
result = bilby.run_sampler(
likelihood, priors, sampler='dynesty', outdir=outdir, label=label,
nlive=1000, walks=100, n_check_point=10000, check_point_plot=True,
conversion_function=bilby.gw.conversion.generate_all_bbh_parameters)
result.plot_corner()