[1]:
! rm visualising_the_results/*
rm: cannot remove 'visualising_the_results/*': No such file or directory

Visualising the results

In this tutorial, we demonstrate the plotting tools built-in to bilby and how to extend them. First, we run a simple injection study and return the result object.

[2]:
import bilby
import matplotlib.pyplot as plt

%matplotlib inline
[3]:
duration = 4
sampling_frequency = 2048
outdir = "visualising_the_results"
label = "example"

injection_parameters = dict(
    mass_1=36.0,
    mass_2=29.0,
    a_1=0.4,
    a_2=0.3,
    tilt_1=0.5,
    tilt_2=1.0,
    phi_12=1.7,
    phi_jl=0.3,
    luminosity_distance=1000.0,
    theta_jn=0.4,
    phase=1.3,
    ra=1.375,
    dec=-1.2108,
    geocent_time=1126259642.413,
    psi=2.659,
)
[4]:
# specify waveform arguments
waveform_arguments = dict(
    waveform_approximant="IMRPhenomXP",  # waveform approximant name
    reference_frequency=50.0,  # gravitational waveform reference frequency (Hz)
)

# set up the waveform generator
waveform_generator = bilby.gw.waveform_generator.WaveformGenerator(
    sampling_frequency=sampling_frequency,
    duration=duration,
    frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole,
    parameters=injection_parameters,
    waveform_arguments=waveform_arguments,
)
04:48 bilby INFO    : Waveform generator initiated with
  frequency_domain_source_model: bilby.gw.source.lal_binary_black_hole
  time_domain_source_model: None
  parameter_conversion: bilby.gw.conversion.convert_to_lal_binary_black_hole_parameters
[5]:
ifos = bilby.gw.detector.InterferometerList(["H1", "L1"])
ifos.set_strain_data_from_power_spectral_densities(
    duration=duration,
    sampling_frequency=sampling_frequency,
    start_time=injection_parameters["geocent_time"] - 2,
)
_ = ifos.inject_signal(
    waveform_generator=waveform_generator, parameters=injection_parameters
)
04:48 bilby INFO    : Injected signal in H1:
04:48 bilby INFO    :   optimal SNR = 23.78
04:48 bilby INFO    :   matched filter SNR = 23.50+0.79j
04:48 bilby INFO    :   mass_1 = 36.0
04:48 bilby INFO    :   mass_2 = 29.0
04:48 bilby INFO    :   a_1 = 0.4
04:48 bilby INFO    :   a_2 = 0.3
04:48 bilby INFO    :   tilt_1 = 0.5
04:48 bilby INFO    :   tilt_2 = 1.0
04:48 bilby INFO    :   phi_12 = 1.7
04:48 bilby INFO    :   phi_jl = 0.3
04:48 bilby INFO    :   luminosity_distance = 1000.0
04:48 bilby INFO    :   theta_jn = 0.4
04:48 bilby INFO    :   phase = 1.3
04:48 bilby INFO    :   ra = 1.375
04:48 bilby INFO    :   dec = -1.2108
04:48 bilby INFO    :   geocent_time = 1126259642.413
04:48 bilby INFO    :   psi = 2.659
04:48 bilby INFO    : Injected signal in L1:
04:48 bilby INFO    :   optimal SNR = 19.25
04:48 bilby INFO    :   matched filter SNR = 21.64+0.70j
04:48 bilby INFO    :   mass_1 = 36.0
04:48 bilby INFO    :   mass_2 = 29.0
04:48 bilby INFO    :   a_1 = 0.4
04:48 bilby INFO    :   a_2 = 0.3
04:48 bilby INFO    :   tilt_1 = 0.5
04:48 bilby INFO    :   tilt_2 = 1.0
04:48 bilby INFO    :   phi_12 = 1.7
04:48 bilby INFO    :   phi_jl = 0.3
04:48 bilby INFO    :   luminosity_distance = 1000.0
04:48 bilby INFO    :   theta_jn = 0.4
04:48 bilby INFO    :   phase = 1.3
04:48 bilby INFO    :   ra = 1.375
04:48 bilby INFO    :   dec = -1.2108
04:48 bilby INFO    :   geocent_time = 1126259642.413
04:48 bilby INFO    :   psi = 2.659
[6]:
# first, set up all priors to be equal to a delta function at their designated value
priors = bilby.gw.prior.BBHPriorDict(injection_parameters.copy())
# then, reset the priors on the masses and luminosity distance to conduct a search over these parameters
priors["mass_1"] = bilby.core.prior.Uniform(25, 40, "mass_1")
priors["mass_2"] = bilby.core.prior.Uniform(25, 40, "mass_2")
priors["luminosity_distance"] = bilby.core.prior.Uniform(
    400, 2000, "luminosity_distance"
)
[7]:
# compute the likelihoods
likelihood = bilby.gw.likelihood.GravitationalWaveTransient(
    interferometers=ifos, waveform_generator=waveform_generator
)
[8]:
result = bilby.core.sampler.run_sampler(
    likelihood=likelihood,
    priors=priors,
    sampler="dynesty",
    npoints=100,
    injection_parameters=injection_parameters,
    outdir=outdir,
    label=label,
    walks=5,
    nact=2,
)
04:48 bilby INFO    : Running for label 'example', output will be saved to 'visualising_the_results'
04:48 bilby INFO    : Using lal version 7.2.2
04:48 bilby INFO    : Using lal git version Branch: None;Tag: lalsuite-v7.8;Id: f7003a0d1b9811e05caba042bccc6baaa4cdd020;;Builder: Unknown User <>;Repository status: CLEAN: All modifications committed
04:48 bilby INFO    : Using lalsimulation version 4.0.2
04:48 bilby INFO    : Using lalsimulation git version Branch: None;Tag: lalsuite-v7.8;Id: f7003a0d1b9811e05caba042bccc6baaa4cdd020;;Builder: Unknown User <>;Repository status: CLEAN: All modifications committed
04:48 bilby INFO    : Analysis priors:
04:48 bilby INFO    : mass_1=Uniform(minimum=25, maximum=40, name='mass_1', latex_label='$m_1$', unit=None, boundary=None)
04:48 bilby INFO    : mass_2=Uniform(minimum=25, maximum=40, name='mass_2', latex_label='$m_2$', unit=None, boundary=None)
04:48 bilby INFO    : luminosity_distance=Uniform(minimum=400, maximum=2000, name='luminosity_distance', latex_label='$d_L$', unit=None, boundary=None)
04:48 bilby INFO    : a_1=0.4
04:48 bilby INFO    : a_2=0.3
04:48 bilby INFO    : tilt_1=0.5
04:48 bilby INFO    : tilt_2=1.0
04:48 bilby INFO    : phi_12=1.7
04:48 bilby INFO    : phi_jl=0.3
04:48 bilby INFO    : theta_jn=0.4
04:48 bilby INFO    : phase=1.3
04:48 bilby INFO    : ra=1.375
04:48 bilby INFO    : dec=-1.2108
04:48 bilby INFO    : geocent_time=1126259642.413
04:48 bilby INFO    : psi=2.659
04:48 bilby INFO    : Analysis likelihood class: <class 'bilby.gw.likelihood.base.GravitationalWaveTransient'>
04:48 bilby INFO    : Analysis likelihood noise evidence: -8456.481498128032
04:48 bilby INFO    : Single likelihood evaluation took 1.232e-02 s
04:48 bilby INFO    : Using sampler Dynesty with kwargs {'bound': 'multi', 'sample': 'rwalk', 'print_progress': True, 'periodic': None, 'reflective': None, 'check_point_delta_t': 1800, 'nlive': 100, 'first_update': None, 'walks': 5, 'npdim': None, 'rstate': None, 'queue_size': 1, 'pool': None, 'use_pool': None, 'live_points': None, 'logl_args': None, 'logl_kwargs': None, 'ptform_args': None, 'ptform_kwargs': None, 'enlarge': 1.5, 'bootstrap': None, 'vol_dec': 0.5, 'vol_check': 8.0, 'facc': 0.2, 'slices': 5, 'update_interval': 60, 'print_func': <bound method Dynesty._print_func of <bilby.core.sampler.dynesty.Dynesty object at 0x7f5f739ae8e0>>, 'dlogz': 0.1, 'maxiter': None, 'maxcall': None, 'logl_max': inf, 'add_live': True, 'save_bounds': False, 'n_effective': None, 'maxmcmc': 5000, 'nact': 2, 'print_method': 'tqdm'}
04:48 bilby INFO    : Checkpoint every check_point_delta_t = 600s
04:48 bilby INFO    : Using dynesty version 1.0.1
04:48 bilby INFO    : Using the bilby-implemented rwalk sample method with ACT estimated walks
04:48 bilby INFO    : Resume file visualising_the_results/example_resume.pickle does not exist.
04:48 bilby INFO    : Generating initial points from the prior
04:51 bilby INFO    : Written checkpoint file visualising_the_results/example_resume.pickle
04:51 bilby INFO    : Sampling time: 0:03:00.707365

04:51 bilby INFO    : Summary of results:
nsamples: 1088
ln_noise_evidence: -8456.481
ln_evidence: -7955.195 +/-  0.354
ln_bayes_factor: 501.287 +/-  0.354

In running this code, we already made the first plot! In the function bilby.detector.get_interferometer_with_fake_noise_and_injection, the ASD, detector data, and signal are plotted together. This figure is saved under visualsing_the_results/H1_frequency_domain_data.png. Note that visualising_the_result is our outdir where all the output of the run is stored. Let’s take a quick look at that directory now:

[9]:
!ls visualising_the_results/
example_checkpoint_run.png         example_dynesty.pickle
example_checkpoint_stats.png       example_result.json
example_checkpoint_trace.png       example_resume.pickle
example_checkpoint_trace_unit.png

Corner plots

Now lets make some corner plots. You can easily generate a corner plot using result.plot_corner() like this:

[10]:
result.plot_corner()
plt.show()
plt.close()

In a notebook, this figure will display. But by default the file is also saved to visualising_the_result/example_corner.png. If you change the label to something more descriptive then the example here will of course be replaced.

Waveform Reconstruction plot

Some plots specific to compact binary coalescence parameter estimation results can be created by re-loading the result as a CBCResult:

[11]:
from bilby.gw.result import CBCResult

cbc_result = CBCResult.from_json("visualising_the_results/example_result.json")
for ifo in ifos:
    cbc_result.plot_interferometer_waveform_posterior(
        interferometer=ifo, n_samples=500, save=False
    )
    plt.show()
    plt.close()
04:51 bilby INFO    : Generating waveform figure for H1
04:51 bilby INFO    : Waveform generator initiated with
  frequency_domain_source_model: bilby.gw.source.lal_binary_black_hole
  time_domain_source_model: None
  parameter_conversion: bilby.gw.conversion.convert_to_lal_binary_black_hole_parameters
_images/visualising_the_results_15_1.png
04:51 bilby INFO    : Generating waveform figure for L1
04:51 bilby INFO    : Waveform generator initiated with
  frequency_domain_source_model: bilby.gw.source.lal_binary_black_hole
  time_domain_source_model: None
  parameter_conversion: bilby.gw.conversion.convert_to_lal_binary_black_hole_parameters
_images/visualising_the_results_15_3.png

Marginal Distribution plots

These plots just show the 1D histograms for each parameter

[12]:
result.plot_marginals()
04:52 bilby INFO    : Plotting mass_1 marginal distribution
04:52 bilby INFO    : Plotting mass_1 marginal distribution
04:52 bilby INFO    : Plotting mass_2 marginal distribution
04:52 bilby INFO    : Plotting mass_2 marginal distribution
04:52 bilby INFO    : Plotting luminosity_distance marginal distribution
04:52 bilby INFO    : Plotting luminosity_distance marginal distribution
04:52 bilby INFO    : Plotting a_1 marginal distribution
04:52 bilby INFO    : Plotting a_1 marginal distribution
04:52 bilby INFO    : Plotting a_2 marginal distribution
04:52 bilby INFO    : Plotting a_2 marginal distribution
04:52 bilby INFO    : Plotting tilt_1 marginal distribution
04:52 bilby INFO    : Plotting tilt_1 marginal distribution
04:52 bilby INFO    : Plotting tilt_2 marginal distribution
04:52 bilby INFO    : Plotting tilt_2 marginal distribution
04:52 bilby INFO    : Plotting phi_12 marginal distribution
04:52 bilby INFO    : Plotting phi_12 marginal distribution
04:52 bilby INFO    : Plotting phi_jl marginal distribution
04:52 bilby INFO    : Plotting phi_jl marginal distribution
04:52 bilby INFO    : Plotting theta_jn marginal distribution
04:52 bilby INFO    : Plotting theta_jn marginal distribution
04:52 bilby INFO    : Plotting phase marginal distribution
04:52 bilby INFO    : Plotting phase marginal distribution
04:52 bilby INFO    : Plotting ra marginal distribution
04:52 bilby INFO    : Plotting ra marginal distribution
04:52 bilby INFO    : Plotting dec marginal distribution
04:52 bilby INFO    : Plotting dec marginal distribution
04:52 bilby INFO    : Plotting geocent_time marginal distribution
04:52 bilby INFO    : Plotting geocent_time marginal distribution
04:52 bilby INFO    : Plotting psi marginal distribution
04:52 bilby INFO    : Plotting psi marginal distribution
04:52 bilby INFO    : Plotting log_likelihood marginal distribution
04:52 bilby INFO    : Plotting log_likelihood marginal distribution
04:52 bilby INFO    : Plotting log_prior marginal distribution
04:52 bilby INFO    : Plotting log_prior marginal distribution

Customizing corner plots

You may also want to plot a subset of the parameters, or perhaps add the injection_paramters as lines to check if you recovered them correctly. All this can be done through plot_corner. Under the hood, plot_corner uses corner, and all the keyword arguments passed to plot_corner are passed through.

Adding injection parameters to the plot

In the previous plot, you’ll notice bilby added the injection parameters to the plot by default. You can switch this off by setting truth=None when you call plot_corner. Or to add different injection parameters to the plot, just pass this as a keyword argument for truth. In this example, we just add a line for the luminosity distance by passing a dictionary of the value we want to display.

[13]:
result.plot_corner(truth=dict(luminosity_distance=201))
plt.show()
plt.close()

Plot a subset of the corner plot

Or, to plot just a subset of parameters, just pass a list of the names you want.

[14]:
result.plot_corner(
    parameters=["mass_1", "mass_2"], filename="{}/subset.png".format(outdir)
)
plt.show()
plt.close()

Notice here, we also passed in a keyword argument filename=, this overwrites the default filename and instead saves the file as visualising_the_results/subset.png. Useful if you want to create lots of different plots. Let’s check what the outdir looks like now

[15]:
!ls visualising_the_results/
example_1d                         example_corner.png
example_checkpoint_run.png         example_dynesty.pickle
example_checkpoint_stats.png       example_result.json
example_checkpoint_trace.png       example_resume.pickle
example_checkpoint_trace_unit.png  subset.png

Alternative

If you would prefer to do the plotting yourself, you can get hold of the samples and the ordering as follows and then plot with a different module. Here is an example using the corner package

[16]:
import corner

samples = result.samples
labels = result.parameter_labels
fig = corner.corner(samples, labels=labels)
plt.show()
plt.close()
_images/visualising_the_results_25_0.png