From 0a5ea9d5ecb354e9dcf382900bef2fed56d25a67 Mon Sep 17 00:00:00 2001
From: Colm Talbot <colm.talbot@ligo.org>
Date: Mon, 7 Mar 2022 19:27:05 +0000
Subject: [PATCH] Review tutorials
---
.gitlab-ci.yml | 15 +-
.pre-commit-config.yaml | 19 ++
bilby/bilby_mcmc/sampler.py | 7 +-
bilby/core/sampler/dynesty.py | 11 +-
bilby/core/sampler/emcee.py | 16 +-
bilby/core/sampler/ptemcee.py | 45 ++-
docs/conf.py | 2 +-
.../tutorials/basic_ptmcmc_tutorial.ipynb | 283 ---------------
examples/tutorials/compare_samplers.ipynb | 321 ++++++++++--------
.../fitting_with_x_and_y_errors.ipynb | 251 +++++++-------
examples/tutorials/making_priors.ipynb | 147 ++++----
.../tutorials/visualising_the_results.ipynb | 278 ++++++++-------
test/core/sampler/dynesty_test.py | 2 -
test/gw/utils_test.py | 6 +-
test/integration/sampler_run_test.py | 1 +
15 files changed, 596 insertions(+), 808 deletions(-)
delete mode 100644 examples/tutorials/basic_ptmcmc_tutorial.ipynb
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 57415e7d6..3d7920b1f 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -211,17 +211,16 @@ plotting-python-3.9:
docs:
stage: docs
image: containers.ligo.org/lscsoft/bilby/v2-bilby-python39
+ before_script:
+ - conda install -c conda-forge pandoc -y
+ - python -m pip install ipykernel ipython jupyter
+ - python -m ipykernel install
script:
# Make the documentation
- - apt-get -yqq install pandoc
- python -m pip install .
- - cd docs
- - pip install ipykernel ipython jupyter
- - cp ../examples/tutorials/*.ipynb ./
- - rm basic_ptmcmc_tutorial.ipynb
- - rm compare_samplers.ipynb
- - rm visualising_the_results.ipynb
- - jupyter nbconvert --to notebook --execute *.ipynb --inplace
+ - cd examples/tutorials
+ - jupyter nbconvert --to notebook --execute *.ipynb --output-dir ../../docs
+ - cd ../../docs
- make clean
- make html
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 3672d0fc0..137e5813c 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -21,3 +21,22 @@ repos:
- id: isort # sort imports alphabetically and separates import into sections
args: [-w=88, -m=3, -tc, -sp=setup.cfg ]
files: '(^bilby/bilby_mcmc/|^examples/core_examples/examples/gw_examples/data_examples)'
+- repo: https://github.com/datarootsio/databooks
+ rev: 0.1.14
+ hooks:
+ - id: databooks
+ name: databooks
+ description:
+ "Remove notebook metadata using `databooks`."
+ entry: databooks meta
+ language: python
+ minimum_pre_commit_version: 2.9.2
+ types: [ jupyter ]
+ args: [-w]
+- repo: local
+ hooks:
+ - id: jupyter-nb-clear-output
+ name: jupyter-nb-clear-output
+ files: \.ipynb$
+ language: system
+ entry: jupyter nbconvert --ClearOutputPreprocessor.enabled=True --inplace
diff --git a/bilby/bilby_mcmc/sampler.py b/bilby/bilby_mcmc/sampler.py
index a1446e4be..9ea0671e5 100644
--- a/bilby/bilby_mcmc/sampler.py
+++ b/bilby/bilby_mcmc/sampler.py
@@ -108,6 +108,8 @@ class Bilby_MCMC(MCMCSampler):
evidence_method: str, [stepping_stone, thermodynamic]
The evidence calculation method to use. Defaults to stepping_stone, but
the results of all available methods are stored in the ln_z_dict.
+ verbose: bool
+ Whether to print diagnostic output during the run.
"""
@@ -152,6 +154,7 @@ class Bilby_MCMC(MCMCSampler):
diagnostic=False,
resume=True,
exit_code=130,
+ verbose=False,
**kwargs,
):
@@ -189,6 +192,7 @@ class Bilby_MCMC(MCMCSampler):
self.resume_file = "{}/{}_resume.pickle".format(self.outdir, self.label)
self.verify_configuration()
+ self.verbose = verbose
try:
signal.signal(signal.SIGTERM, self.write_current_state_and_exit)
@@ -486,7 +490,8 @@ class Bilby_MCMC(MCMCSampler):
rse = 100 * count / nsamples
msg += f"|rse={rse:0.2f}%"
- print(msg, flush=True)
+ if self.verbose:
+ print(msg, flush=True)
def print_per_proposal(self):
logger.info("Zero-temperature proposals:")
diff --git a/bilby/core/sampler/dynesty.py b/bilby/core/sampler/dynesty.py
index 154866bec..d9f1ae94b 100644
--- a/bilby/core/sampler/dynesty.py
+++ b/bilby/core/sampler/dynesty.py
@@ -107,8 +107,9 @@ class Dynesty(NestedSampler):
`ndim * 10` can be a reasonable rule of thumb for new problems.
dlogz: float, (0.1)
Stopping criteria
- verbose: Bool
- If true, print information information about the convergence during
+ print_progress: Bool
+ If true, print information information about the convergence during.
+ `verbose` has the same effect.
check_point: bool,
If true, use check pointing.
check_point_plot: bool,
@@ -130,7 +131,7 @@ class Dynesty(NestedSampler):
- else: print to `stdout` at every iteration
"""
default_kwargs = dict(bound='multi', sample='rwalk',
- verbose=True, periodic=None, reflective=None,
+ periodic=None, reflective=None,
check_point_delta_t=1800, nlive=1000,
first_update=None, walks=100,
npdim=None, rstate=None, queue_size=1, pool=None,
@@ -230,7 +231,7 @@ class Dynesty(NestedSampler):
if self.kwargs['print_func'] is None:
self.kwargs['print_func'] = self._print_func
print_method = self.kwargs["print_method"]
- if print_method == "tqdm":
+ if print_method == "tqdm" and self.kwargs["print_progress"]:
self.pbar = tqdm(file=sys.stdout)
elif "interval" in print_method:
self._last_print_time = datetime.datetime.now()
@@ -407,7 +408,7 @@ class Dynesty(NestedSampler):
self._close_pool()
# Flushes the output to force a line break
- if self.kwargs["verbose"] and self.kwargs["print_method"] == "tqdm":
+ if self.kwargs["print_progress"] and self.kwargs["print_method"] == "tqdm":
self.pbar.close()
print("")
diff --git a/bilby/core/sampler/emcee.py b/bilby/core/sampler/emcee.py
index 1f9ba786c..d976d391e 100644
--- a/bilby/core/sampler/emcee.py
+++ b/bilby/core/sampler/emcee.py
@@ -39,6 +39,8 @@ class Emcee(MCMCSampler):
The number of autocorrelation times to discard as burn-in
a: float (2)
The proposal scale factor
+ verbose: bool
+ Whether to print diagnostic information during the analysis
"""
@@ -51,7 +53,7 @@ class Emcee(MCMCSampler):
def __init__(self, likelihood, priors, outdir='outdir', label='label',
use_ratio=False, plot=False, skip_import_verification=False,
pos0=None, nburn=None, burn_in_fraction=0.25, resume=True,
- burn_in_act=3, **kwargs):
+ burn_in_act=3, verbose=True, **kwargs):
import emcee
self.emcee = emcee
@@ -69,6 +71,7 @@ class Emcee(MCMCSampler):
self.nburn = nburn
self.burn_in_fraction = burn_in_fraction
self.burn_in_act = burn_in_act
+ self.verbose = verbose
signal.signal(signal.SIGTERM, self.checkpoint_and_exit)
signal.signal(signal.SIGINT, self.checkpoint_and_exit)
@@ -361,10 +364,15 @@ class Emcee(MCMCSampler):
sampler_function_kwargs['p0'] = self.pos0
# main iteration loop
- for sample in tqdm(
- self.sampler.sample(iterations=iterations, **sampler_function_kwargs),
- total=iterations):
+ iterator = self.sampler.sample(
+ iterations=iterations, **sampler_function_kwargs
+ )
+ if self.verbose:
+ iterator = tqdm(iterator, total=iterations)
+ for sample in iterator:
self.write_chains_to_file(sample)
+ if self.verbose:
+ iterator.close()
self.checkpoint()
self.result.sampler_output = np.nan
diff --git a/bilby/core/sampler/ptemcee.py b/bilby/core/sampler/ptemcee.py
index 4e60496fe..bce289e95 100644
--- a/bilby/core/sampler/ptemcee.py
+++ b/bilby/core/sampler/ptemcee.py
@@ -110,7 +110,7 @@ class Ptemcee(MCMCSampler):
Other Parameters
- ------==========
+ ================
nwalkers: int, (200)
The number of walkers
nsteps: int, (100)
@@ -168,6 +168,7 @@ class Ptemcee(MCMCSampler):
pos0="prior",
niterations_per_check=5,
log10beta_min=None,
+ verbose=True,
**kwargs
):
super(Ptemcee, self).__init__(
@@ -249,6 +250,7 @@ class Ptemcee(MCMCSampler):
betas = np.logspace(0, self.log10beta_min, self.ntemps)
logger.warning("Using betas {}".format(betas))
self.kwargs["betas"] = betas
+ self.verbose = verbose
@property
def sampler_function_kwargs(self):
@@ -555,6 +557,7 @@ class Ptemcee(MCMCSampler):
self.tau_list_n,
self.Q_list,
self.mean_log_posterior,
+ verbose=self.verbose,
)
if stop:
@@ -705,6 +708,7 @@ def check_iteration(
tau_list_n,
Q_list,
mean_log_posterior,
+ verbose=True,
):
""" Per-iteration logic to calculate the convergence check
@@ -716,6 +720,8 @@ def check_iteration(
A list of the search parameter keys
time_per_check, tau_list, tau_list_n: list
Lists used for tracking the run
+ verbose: bool
+ Whether to print the output
Returns
=======
@@ -754,10 +760,11 @@ def check_iteration(
Q_list.append(Q)
if np.isnan(tau) or np.isinf(tau):
- print_progress(
- iteration, sampler, time_per_check, np.nan, np.nan,
- np.nan, np.nan, np.nan, False, convergence_inputs, Q,
- )
+ if verbose:
+ print_progress(
+ iteration, sampler, time_per_check, np.nan, np.nan,
+ np.nan, np.nan, np.nan, False, convergence_inputs, Q,
+ )
return False, np.nan, np.nan, np.nan, np.nan
# Convert to an integer
@@ -827,19 +834,20 @@ def check_iteration(
tau_usable = False
# Print an update on the progress
- print_progress(
- iteration,
- sampler,
- time_per_check,
- nsamples_effective,
- samples_per_check,
- tau_int,
- gradient_tau,
- gradient_mean_log_posterior,
- tau_usable,
- convergence_inputs,
- Q
- )
+ if verbose:
+ print_progress(
+ iteration,
+ sampler,
+ time_per_check,
+ nsamples_effective,
+ samples_per_check,
+ tau_int,
+ gradient_tau,
+ gradient_mean_log_posterior,
+ tau_usable,
+ convergence_inputs,
+ Q
+ )
stop = converged and tau_usable
return stop, nburn, thin, tau_int, nsamples_effective
@@ -1145,6 +1153,7 @@ def compute_evidence(sampler, log_likelihood_array, outdir, label, discard, nbur
ax2.set_xlabel(r"$\beta_{min}$")
plt.tight_layout()
fig.savefig("{}/{}_beta_lnl.png".format(outdir, label))
+ plt.close(fig)
return lnZ, lnZerr
diff --git a/docs/conf.py b/docs/conf.py
index 8b130c4af..dcd3e028a 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -48,7 +48,7 @@ templates_path = ['templates']
# The suffix(es) of source filenames.
# You can specify multiple suffix as a list of string:
#
-source_suffix = ['.rst', '.md', '.txt', '.ipynb']
+source_suffix = ['.rst', '.md', '.txt']
# The master toctree document.
master_doc = 'index'
diff --git a/examples/tutorials/basic_ptmcmc_tutorial.ipynb b/examples/tutorials/basic_ptmcmc_tutorial.ipynb
deleted file mode 100644
index c36fe65c2..000000000
--- a/examples/tutorials/basic_ptmcmc_tutorial.ipynb
+++ /dev/null
@@ -1,283 +0,0 @@
-{
- "cells": [
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "from __future__ import division, print_function\n",
- "import matplotlib.pyplot as plt\n",
- "import numpy as np\n",
- "import seaborn as sns\n",
- "import bilby\n",
- "%matplotlib inline"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## This notebook will show how to use the PTMCMCSampler, in particular this will highlight how to add custom jump proposals"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## create 150914 like injection"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# Set the duration and sampling frequency of the data segment that we're\n",
- "# going to inject the signal into\n",
- "duration = 4.\n",
- "sampling_frequency = 2048.\n",
- "\n",
- "# Specify the output directory and the name of the simulation.\n",
- "outdir = 'outdir'\n",
- "label = 'basic_tutorial4'\n",
- "bilby.core.utils.setup_logger(outdir=outdir, label=label)\n",
- "\n",
- "# Set up a random seed for result reproducibility. This is optional!\n",
- "np.random.seed(88170235)\n",
- "\n",
- "# We are going to inject a binary black hole waveform. We first establish a\n",
- "# dictionary of parameters that includes all of the different waveform\n",
- "# parameters, including masses of the two black holes (mass_1, mass_2),\n",
- "# spins of both black holes (a, tilt, phi), etc.\n",
- "injection_parameters = dict(\n",
- " mass_1=36., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.5, tilt_2=1.0,\n",
- " phi_12=1.7, phi_jl=0.3, luminosity_distance=2000., theta_jn=0.4, psi=2.659,\n",
- " phase=1.3, geocent_time=1126259642.413, ra=1.375, dec=-1.2108)\n",
- "\n",
- "# Fixed arguments passed into the source model\n",
- "waveform_arguments = dict(waveform_approximant='IMRPhenomP',\n",
- " reference_frequency=50., minimum_frequency=20.)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Inject into data"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# Create the waveform_generator using a LAL BinaryBlackHole source function\n",
- "waveform_generator = bilby.gw.WaveformGenerator(\n",
- " duration=duration, sampling_frequency=sampling_frequency,\n",
- " frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole,\n",
- " waveform_arguments=waveform_arguments)\n",
- "\n",
- "# Set up interferometers. In this case we'll use two interferometers\n",
- "# (LIGO-Hanford (H1), LIGO-Livingston (L1). These default to their design\n",
- "# sensitivity\n",
- "\n",
- "ifos = bilby.gw.detector.InterferometerList(['H1', 'L1'])\n",
- "ifos.set_strain_data_from_power_spectral_densities(\n",
- " sampling_frequency=sampling_frequency, duration=duration,\n",
- " start_time=injection_parameters['geocent_time'] - 3)\n",
- "ifos.inject_signal(waveform_generator=waveform_generator,\n",
- " parameters=injection_parameters)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### For simplicity, we will fix all parameters here to the injected value and only vary over mass1 and mass2,"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "priors = injection_parameters.copy()\n",
- "priors['mass_1'] = bilby.prior.Uniform(name='mass_1', minimum=10, maximum=80, unit=r'$M_{\\\\odot}$')\n",
- "priors['mass_2'] = bilby.prior.Uniform(name='mass_1', minimum=10, maximum=80, unit=r'$M_{\\\\odot}$')"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Here we create arbitrary jump proposals. This will highlight the necessary features of a jump proposal in ptmcmc. That is it takes the current position, x, then outputs a new position , q, and the jump probability i.e. p(x -> q). These will then be passed to the standard metropolis hastings condition. \n",
- "## The two proposals below are probably not very good ones, ideally we would use proposals based upon our kmowledge of the problem/parameter space. In general for these proposals lqxy will certainly not be 0 "
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "class UniformJump(object):\n",
- " def __init__(self, pmin, pmax):\n",
- " \"\"\"Draw random parameters from pmin, pmax\"\"\"\n",
- " self.pmin = pmin\n",
- " self.pmax = pmax\n",
- " \n",
- " def unjump(self, x, it, beta):\n",
- " \"\"\" \n",
- " Function prototype must read in parameter vector x,\n",
- " sampler iteration number it, and inverse temperature beta\n",
- " \"\"\"\n",
- " # log of forward-backward jump probability\n",
- " lqxy = 0\n",
- " \n",
- " # uniformly drawn parameters\n",
- " q = np.random.uniform(self.pmin, self.pmax, len(x))\n",
- " return q, lqxy"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "class NormJump(object):\n",
- " def __init__(self, step_size):\n",
- " \"\"\"Draw random parameters from pmin, pmax\"\"\"\n",
- " self.step_size = step_size\n",
- " \n",
- " def normjump(self, x, it, beta):\n",
- " \"\"\" \n",
- " Function prototype must read in parameter vector x,\n",
- " sampler iteration number it, and inverse temperature beta\n",
- " \"\"\"\n",
- " # log of forward-backward jump probability. this is only zero for simple examples.\n",
- " lqxy = 0\n",
- " \n",
- " # uniformly drawn parameters\n",
- " q = np.random.multivariate_normal(x , self.step_size * np.eye(len(x)) , 1)\n",
- " return q[0], lqxy"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### Below we create a dictionary containing our jump proposals and the relative weight of that proposal in the proposal cycle, these are then passed to bilby.run_sampler under the keyword argument custom_proposals = "
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "normjump = NormJump(1)\n",
- "normweight = 5\n",
- "ujump = UniformJump(20, 40)\n",
- "uweight = 1 \n",
- "custom = {'uniform': [ujump.unjump , uweight],\n",
- " 'normal': [normjump.normjump , normweight]}"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# Initialise the likelihood by passing in the interferometer data (ifos) and\n",
- "# the waveoform generator\n",
- "likelihood = bilby.gw.GravitationalWaveTransient(\n",
- " interferometers=ifos,waveform_generator=waveform_generator)\n",
- "result = bilby.run_sampler(\n",
- " likelihood=likelihood, priors=priors, sampler= 'PTMCMCsampler',custom_proposals = custom , Niter = 10**4 )"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "result.plot_corner()"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### PTMCMC produces the acceptance rate for each of the proposals (including the ones built in). This is taken as an average at a specified checkpoint. This is one (acceptnace rate is certainly not the only/even the best metric here. Think exploration v exploitation problem ) indicators of whether our jump proposal is a good one"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "sampler_meta = result.meta_data['sampler_meta']\n",
- "jumps = sampler_meta['proposals']"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "plt.figure()\n",
- "plt.xlabel('epoch')\n",
- "plt.ylabel('acceptance rate')\n",
- "for i,proposal in enumerate(jumps): \n",
- " plt.plot(jumps[proposal] , label = proposal)\n",
- "plt.legend(loc='best', frameon=True)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## We can generate the 1d chains for each of the parameters too and the likelihood of those points on the chain"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "m2 = result.posterior.mass_2.values\n",
- "m1 = result.posterior.mass_1.values\n",
- "\n",
- "fig, ax = plt.subplots(nrows = 2 , ncols =1 , sharex = True , figsize = (8,8))\n",
- "ax[0].plot(m1 , 'o', label = 'm1' )\n",
- "ax[0].plot(m2 , 'o', label = 'm2' )\n",
- "ax[0].set_ylabel(r'$M_{\\odot}$')\n",
- "ax[0].legend(loc = 'best' , frameon = True , fontsize = 12)\n",
- "ax[1].plot(result.log_likelihood_evaluations)\n",
- "ax[1].set_ylabel(r'$\\mathcal{L}$')\n",
- "ax[1].set_xlabel('iterations')\n",
- "ax[1].set_xscale('log')"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": []
- }
- ],
- "metadata": {},
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/examples/tutorials/compare_samplers.ipynb b/examples/tutorials/compare_samplers.ipynb
index e65eef499..ad214ebb8 100644
--- a/examples/tutorials/compare_samplers.ipynb
+++ b/examples/tutorials/compare_samplers.ipynb
@@ -6,7 +6,9 @@
"source": [
"# Compare samplers\n",
"\n",
- "In this notebook, we'll compare the different samplers implemented in `bilby`. As of this version, we don't compare the outputs, only how to run them and the timings for their default setup.\n",
+ "In this notebook, we'll compare the different samplers implemented in `Bilby` on a simple linear regression problem.\n",
+ "\n",
+ "This is not an exhaustive set of the implemented samplers, nor of the settings available for each sampler.\n",
"\n",
"## Setup"
]
@@ -14,230 +16,251 @@
{
"cell_type": "code",
"execution_count": null,
- "metadata": {
- "execution": {
- "iopub.execute_input": "2021-02-05T22:05:40.710069Z",
- "iopub.status.busy": "2021-02-05T22:05:40.709587Z",
- "iopub.status.idle": "2021-02-05T22:05:43.017567Z",
- "shell.execute_reply": "2021-02-05T22:05:43.018795Z"
- }
- },
+ "metadata": {},
"outputs": [],
"source": [
- "import numpy as np\n",
- "import pylab as plt\n",
- "\n",
- "%load_ext autoreload\n",
- "%autoreload 2\n",
- "\n",
"import bilby\n",
+ "import numpy as np\n",
+ "import matplotlib.pyplot as plt\n",
"\n",
- "bilby.utils.setup_logger()\n",
- "\n",
- "time_duration = 1. # set the signal duration (seconds)\n",
- "sampling_frequency = 4096. # set the data sampling frequency (Hz)\n",
- "\n",
- "injection_parameters = dict(\n",
- "chirp_mass=36., # detector frame (redshifted) primary mass (solar masses)\n",
- "mass_ratio=0.9, # detector frame (redshifted) secondary mass (solar masses)\n",
- "a_1=0, # primary dimensionless spin magnitude\n",
- "a_2=0, # secondary dimensionless spin magnitude\n",
- "tilt_1=0, # polar angle between primary spin and the orbital angular momentum (radians)\n",
- "tilt_2=0, # polar angle between secondary spin and the orbital angular momentum \n",
- "phi_12=0, # azimuthal angle between primary and secondary spin (radians)\n",
- "phi_jl=0, # azimuthal angle between total angular momentum and orbital angular momentum (radians)\n",
- "luminosity_distance=100., # luminosity distance to source (Mpc)\n",
- "theta_jn=0.4, # angle between the total angular momentum (both spin and orbital) and the line of sight\n",
- "phase=1.3, # phase (radians)\n",
- "ra=1.375, # source right ascension (radians)\n",
- "dec=-1.2108, # source declination (radians)\n",
- "geocent_time=1126259642.413, # reference time at geocentre (time of coalescence or peak amplitude) (GPS seconds)\n",
- "psi=2.659 # gravitational wave polarisation angle\n",
- ")\n",
- "\n",
- "\n",
- "# initialise the waveform generator\n",
- "waveform_generator = bilby.gw.waveform_generator.WaveformGenerator(\n",
- " sampling_frequency=sampling_frequency,\n",
- " duration=time_duration,\n",
- " frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole,\n",
- " parameters=injection_parameters)\n",
- "\n",
- "# generate a frequency-domain waveform\n",
- "hf_signal = waveform_generator.frequency_domain_strain()\n",
- "\n",
- "# initialise a single interferometer representing LIGO Hanford\n",
- "H1 = bilby.gw.detector.get_empty_interferometer('H1')\n",
- "# set the strain data at the interferometer\n",
- "H1.set_strain_data_from_power_spectral_density(sampling_frequency=sampling_frequency, duration=time_duration)\n",
- "# inject the gravitational wave signal into the interferometer model\n",
- "H1.inject_signal(injection_polarizations=hf_signal, parameters=injection_parameters)\n",
- "\n",
- "IFOs = [H1]\n",
- "\n",
- "# compute the likelihood on each of the signal parameters\n",
- "likelihood = bilby.gw.likelihood.GravitationalWaveTransient(IFOs, waveform_generator)"
+ "%matplotlib inline"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "label = \"linear_regression\"\n",
+ "outdir = \"outdir\"\n",
+ "bilby.utils.check_directory_exists_and_if_not_mkdir(outdir)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "## Prior\n",
+ "## Define our model\n",
"\n",
- "For this test, we will simply search of the sky position, setting the other parameters to their simulated values."
+ "Here our model is a simple linear fit to some quantity $y = m x + c$."
]
},
{
"cell_type": "code",
"execution_count": null,
- "metadata": {
- "execution": {
- "iopub.execute_input": "2021-02-05T22:05:43.025484Z",
- "iopub.status.busy": "2021-02-05T22:05:43.024858Z",
- "iopub.status.idle": "2021-02-05T22:05:43.101096Z",
- "shell.execute_reply": "2021-02-05T22:05:43.100588Z"
- }
- },
+ "metadata": {},
"outputs": [],
"source": [
- "# set the priors on each of the injection parameters to be a delta function at their given value\n",
- "priors = bilby.gw.prior.BBHPriorDict()\n",
- "for key in injection_parameters.keys():\n",
- " priors[key] = injection_parameters[key]\n",
- "\n",
- "# now reset the priors on the sky position coordinates in order to conduct a sky position search\n",
- "priors['ra'] = bilby.prior.Uniform(0, 2*np.pi, 'ra')\n",
- "priors['dec'] = bilby.prior.Cosine(name='dec', minimum=-np.pi/2, maximum=np.pi/2)"
+ "def model(x, m, c):\n",
+ " return x * m + c"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "## PyMultinest"
+ "## Simulate data\n",
+ "\n",
+ "We simulate observational data.\n",
+ "We assume some uncertainty in the observations and so perturb the observations from the truth."
]
},
{
"cell_type": "code",
"execution_count": null,
- "metadata": {
- "execution": {
- "iopub.execute_input": "2021-02-05T22:05:43.105117Z",
- "iopub.status.busy": "2021-02-05T22:05:43.104639Z",
- "iopub.status.idle": "2021-02-05T22:05:43.272793Z",
- "shell.execute_reply": "2021-02-05T22:05:43.272156Z"
- }
- },
+ "metadata": {},
"outputs": [],
"source": [
- "%%time \n",
- "result = bilby.core.sampler.run_sampler(\n",
- " likelihood, priors=priors, sampler='pymultinest', label='pymultinest',\n",
- " npoints=2000, verbose=False, resume=False)\n",
- "fig = result.plot_corner(save=False)\n",
- "# show the corner plot\n",
+ "injection_parameters = dict(m=0.5, c=0.2)\n",
+ "\n",
+ "sampling_frequency = 10\n",
+ "time_duration = 10\n",
+ "time = np.arange(0, time_duration, 1 / sampling_frequency)\n",
+ "N = len(time)\n",
+ "sigma = np.random.normal(1, 0.01, N)\n",
+ "data = model(time, **injection_parameters) + np.random.normal(0, sigma, N)\n",
+ "\n",
+ "fig, ax = plt.subplots()\n",
+ "ax.plot(time, data, \"x\", label=\"Data\")\n",
+ "ax.plot(time, model(time, **injection_parameters), \"--r\", label=\"Truth\")\n",
+ "ax.set_xlim(0, 10)\n",
+ "ax.set_ylim(-2, 8)\n",
+ "ax.set_xlabel(\"x\")\n",
+ "ax.set_ylabel(\"y\")\n",
+ "ax.legend()\n",
"plt.show()\n",
- "print(result)"
+ "plt.close()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "## dynesty"
+ "## Define the likelihood and prior\n",
+ "\n",
+ "For any Bayesian calculation we need a likelihood and a prior.\n",
+ "\n",
+ "In this case, we take a `GausianLikelihood` as we assume the uncertainty on the data is normally distributed.\n",
+ "\n",
+ "For both of our parameters we take uniform priors."
]
},
{
"cell_type": "code",
"execution_count": null,
- "metadata": {
- "execution": {
- "iopub.execute_input": "2021-02-05T22:05:43.276181Z",
- "iopub.status.busy": "2021-02-05T22:05:43.275690Z",
- "iopub.status.idle": "2021-02-05T22:06:44.384185Z",
- "shell.execute_reply": "2021-02-05T22:06:44.384572Z"
- }
- },
+ "metadata": {},
"outputs": [],
"source": [
- "%%time \n",
- "result = bilby.core.sampler.run_sampler(\n",
- " likelihood, priors=priors, sampler='dynesty', label='dynesty',\n",
- " bound='multi', sample='rwalk', npoints=200, walks=1, verbose=False,\n",
- " update_interval=100)\n",
- "fig = result.plot_corner(save=False)\n",
- "# show the corner plot\n",
- "plt.show()\n",
- "print(result)"
+ "likelihood = bilby.likelihood.GaussianLikelihood(time, data, model, sigma)\n",
+ "\n",
+ "priors = bilby.core.prior.PriorDict()\n",
+ "priors[\"m\"] = bilby.core.prior.Uniform(0, 5, \"m\")\n",
+ "priors[\"c\"] = bilby.core.prior.Uniform(-2, 2, \"c\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "## Dynamic Nested Sampling (Dynesty)\n",
+ "## Run the samplers and compare the inferred posteriors\n",
+ "\n",
+ "We'll use four of the implemented samplers.\n",
+ "\n",
+ "For each one we specify a set of parameters.\n",
"\n",
- "See [the dynesty docs](http://dynesty.readthedocs.io/en/latest/dynamic.html#). Essentially, this methods improves the posterior estimation over that of standard nested sampling."
+ "`Bilby`/the underlying samplers produce quite a lot of output while the samplers are running so we will suppress as many of these as possible.\n",
+ "\n",
+ "After running the analysis, we print a final summary for each of the samplers."
]
},
{
"cell_type": "code",
"execution_count": null,
- "metadata": {
- "execution": {
- "iopub.execute_input": "2021-02-05T22:06:44.389121Z",
- "iopub.status.busy": "2021-02-05T22:06:44.388707Z",
- "iopub.status.idle": "2021-02-05T22:07:12.688768Z",
- "shell.execute_reply": "2021-02-05T22:07:12.689099Z"
- }
- },
+ "metadata": {},
"outputs": [],
"source": [
- "%%time \n",
- "result = bilby.core.sampler.run_sampler(\n",
- " likelihood, priors=priors, sampler='dynesty', label='dynesty_dynamic',\n",
- " bound='multi', nlive=250, sample='unif', verbose=True,\n",
- " update_interval=100, dynamic=True)\n",
- "fig = result.plot_corner(save=False)\n",
- "# show the corner plot\n",
- "plt.show()\n",
- "print(result)"
+ "samplers = dict(\n",
+ " bilby_mcmc=dict(\n",
+ " nsamples=1000,\n",
+ " L1steps=20,\n",
+ " ntemps=10,\n",
+ " printdt=10,\n",
+ " ),\n",
+ " dynesty=dict(\n",
+ " npoints=500,\n",
+ " walks=5,\n",
+ " nact=3,\n",
+ " ),\n",
+ " pymultinest=dict(nlive=500),\n",
+ " nestle=dict(nlive=500),\n",
+ " emcee=dict(nwalkers=20, iterations=500),\n",
+ " ptemcee=dict(ntemps=10, nwalkers=20, nsamples=1000),\n",
+ ")\n",
+ "\n",
+ "results = dict()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "bilby.core.utils.logger.setLevel(\"ERROR\")\n",
+ "\n",
+ "for sampler in samplers:\n",
+ " result = bilby.core.sampler.run_sampler(\n",
+ " likelihood,\n",
+ " priors=priors,\n",
+ " sampler=sampler,\n",
+ " label=sampler,\n",
+ " resume=False,\n",
+ " clean=True,\n",
+ " verbose=False,\n",
+ " **samplers[sampler]\n",
+ " )\n",
+ " results[sampler] = result"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "print(\"=\" * 40)\n",
+ "for sampler in results:\n",
+ " print(sampler)\n",
+ " print(\"=\" * 40)\n",
+ " print(results[sampler])\n",
+ " print(\"=\" * 40)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "## ptemcee"
+ "## Make comparison plots\n",
+ "\n",
+ "We will make two standard comparison plots.\n",
+ "\n",
+ "In the first we plot the one- and two-dimensional marginal posterior distributions in a \"corner\" plot.\n",
+ "\n",
+ "In the second, we show the inferred model that we are fitting along with the uncertainty by taking random draws from the posterior distribution.\n",
+ "This kind of posterior predicitive plot is useful to identify model misspecification."
]
},
{
"cell_type": "code",
"execution_count": null,
- "metadata": {
- "execution": {
- "iopub.execute_input": "2021-02-05T22:07:12.692452Z",
- "iopub.status.busy": "2021-02-05T22:07:12.692016Z",
- "iopub.status.idle": "2021-02-05T22:12:28.559510Z",
- "shell.execute_reply": "2021-02-05T22:12:28.560201Z"
- }
- },
+ "metadata": {},
"outputs": [],
"source": [
- "%%time \n",
- "result = bilby.core.sampler.run_sampler(\n",
- " likelihood, priors=priors, sampler='ptemcee', label='ptemcee',\n",
- " nwalkers=100, nsteps=200, nburn=100, ntemps=2,\n",
- " tqdm='tqdm_notebook')\n",
- "fig = result.plot_corner(save=False)\n",
- "# show the corner plot\n",
+ "_ = bilby.core.result.plot_multiple(\n",
+ " list(results.values()), labels=list(results.keys()), save=False\n",
+ ")\n",
"plt.show()\n",
- "print(result)"
+ "plt.close()"
]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "fig, ax = plt.subplots(figsize=(12, 8))\n",
+ "ax.plot(time, data, \"x\", label=\"Data\", color=\"r\")\n",
+ "ax.plot(\n",
+ " time, model(time, **injection_parameters), linestyle=\"--\", color=\"k\", label=\"Truth\"\n",
+ ")\n",
+ "\n",
+ "for jj, sampler in enumerate(samplers):\n",
+ " result = results[sampler]\n",
+ " samples = result.posterior[result.search_parameter_keys].sample(500)\n",
+ " for ii in range(len(samples)):\n",
+ " parameters = dict(samples.iloc[ii])\n",
+ " plt.plot(time, model(time, **parameters), color=f\"C{jj}\", alpha=0.01)\n",
+ " plt.axhline(-10, color=f\"C{jj}\", label=sampler.replace(\"_\", \" \"))\n",
+ "ax.set_xlim(0, 10)\n",
+ "ax.set_ylim(-2, 8)\n",
+ "ax.set_xlabel(\"Time\")\n",
+ "ax.set_ylabel(\"y\")\n",
+ "ax.legend(loc=\"upper left\")\n",
+ "plt.show()\n",
+ "plt.close()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
}
],
"metadata": {},
"nbformat": 4,
- "nbformat_minor": 2
+ "nbformat_minor": 4
}
diff --git a/examples/tutorials/fitting_with_x_and_y_errors.ipynb b/examples/tutorials/fitting_with_x_and_y_errors.ipynb
index 7318e4200..6cda25778 100644
--- a/examples/tutorials/fitting_with_x_and_y_errors.ipynb
+++ b/examples/tutorials/fitting_with_x_and_y_errors.ipynb
@@ -4,9 +4,22 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "# Fitting a model to data with both x and y errors with `BILBY`\n",
+ "# Fitting a model to data with both x and y errors with `Bilby`\n",
"\n",
- "Usually when we fit a model to data with a Gaussian Likelihood we assume that we know x values exactly. This is almost never the case. Here we show how to fit a model with errors in both x and y."
+ "Usually when we fit a model to data with a Gaussian Likelihood we assume that we know x values exactly. This is almost never the case. Here we show how to fit a model with errors in both x and y.\n",
+ "\n",
+ "Since we are using a very simple model we will use the `nestle` sampler.\n",
+ "This can be installed using\n",
+ "\n",
+ "```console\n",
+ "$ conda install -c conda-forge nestle\n",
+ "```\n",
+ "\n",
+ "or\n",
+ "\n",
+ "```console\n",
+ "$ pip install nestle\n",
+ "```"
]
},
{
@@ -16,15 +29,19 @@
"outputs": [],
"source": [
"import bilby\n",
- "import inspect\n",
- "%pylab inline"
+ "import matplotlib.pyplot as plt\n",
+ "import numpy as np\n",
+ "\n",
+ "%matplotlib inline"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "### First we create the data and plot it"
+ "### Simulate data\n",
+ "\n",
+ "First we create the data and plot it"
]
},
{
@@ -33,39 +50,51 @@
"metadata": {},
"outputs": [],
"source": [
- "#define our model, a line\n",
+ "# define our model, a line\n",
"def model(x, m, c, **kwargs):\n",
- " y = m*x + c\n",
+ " y = m * x + c\n",
" return y\n",
"\n",
- "#make a function to create and plot our data\n",
- "def make_data(points, m , c, xerr, yerr, seed):\n",
+ "\n",
+ "# make a function to create and plot our data\n",
+ "def make_data(points, m, c, xerr, yerr, seed):\n",
" np.random.seed(int(seed))\n",
- " xtrue = np.linspace(0,100,points)\n",
- " ytrue = model(x = xtrue, m = m, c = c)\n",
+ " xtrue = np.linspace(0, 100, points)\n",
+ " ytrue = model(x=xtrue, m=m, c=c)\n",
"\n",
" xerr = xerr * np.random.randn(points)\n",
" yerr = yerr * np.random.randn(points)\n",
" xobs = xtrue + xerr\n",
" yobs = ytrue + yerr\n",
- " \n",
- " plt.errorbar(xobs, yobs, xerr = xerr, yerr = yerr, fmt = 'x')\n",
- " plt.errorbar(xtrue, ytrue, yerr = yerr, color = 'black', alpha = 0.5)\n",
- " plt.xlim(0,100)\n",
+ "\n",
+ " plt.errorbar(xobs, yobs, xerr=xerr, yerr=yerr, fmt=\"x\")\n",
+ " plt.errorbar(xtrue, ytrue, yerr=yerr, color=\"black\", alpha=0.5)\n",
+ " plt.xlim(0, 100)\n",
" plt.show()\n",
- " \n",
- " data = {'xtrue': xtrue, 'ytrue':ytrue, 'xobs':xobs, 'yobs':yobs, 'xerr':xerr, 'yerr':yerr}\n",
- " \n",
+ " plt.close()\n",
+ "\n",
+ " data = {\n",
+ " \"xtrue\": xtrue,\n",
+ " \"ytrue\": ytrue,\n",
+ " \"xobs\": xobs,\n",
+ " \"yobs\": yobs,\n",
+ " \"xerr\": xerr,\n",
+ " \"yerr\": yerr,\n",
+ " }\n",
+ "\n",
" return data\n",
"\n",
- "data = make_data(points = 30, m = 5, c = 10, xerr = 5, yerr = 5, seed = 123)"
+ "\n",
+ "data = make_data(points=30, m=5, c=10, xerr=5, yerr=5, seed=123)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "### Now lets set up the prior and bilby output directory"
+ "### Define our prior and sampler settings\n",
+ "\n",
+ "Now lets set up the prior and bilby output directory/sampler settings"
]
},
{
@@ -74,20 +103,21 @@
"metadata": {},
"outputs": [],
"source": [
- "#setting up bilby priors\n",
- "priors = dict(m=bilby.core.prior.Uniform(0, 30, 'm'),\n",
- " c=bilby.core.prior.Uniform(0, 30, 'c'))\n",
+ "# setting up bilby priors\n",
+ "priors = dict(\n",
+ " m=bilby.core.prior.Uniform(0, 30, \"m\"), c=bilby.core.prior.Uniform(0, 30, \"c\")\n",
+ ")\n",
"\n",
- "outdir = 'outdir'\n",
- "livepoints = 100\n",
- "walks = 100"
+ "sampler_kwargs = dict(priors=priors, sampler=\"nestle\", nlive=1000, outdir=\"outdir\", verbose=False)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "### Our first step is to recover the straight line using a simple Gaussian Likelihood that only takes into account the y errors. Under the assumption we know x exactly. In this case, we pass in xtrue for x"
+ "### Fit with exactly known x-values\n",
+ "\n",
+ "Our first step is to recover the straight line using a simple Gaussian Likelihood that only takes into account the y errors. Under the assumption we know x exactly. In this case, we pass in xtrue for x"
]
},
{
@@ -96,10 +126,14 @@
"metadata": {},
"outputs": [],
"source": [
- "OneDGaussian_knownx = bilby.core.likelihood.GaussianLikelihood(x = data['xtrue'], y = data['yobs'], func = model, sigma = data['yerr'])\n",
- "result_1D_xtrue = bilby.run_sampler(\n",
- " likelihood=OneDGaussian_knownx, priors=priors, sampler='dynesty', npoints=livepoints,\n",
- " walks=walks, outdir=outdir, label='xtrue_1D_Gaussian')"
+ "known_x = bilby.core.likelihood.GaussianLikelihood(\n",
+ " x=data[\"xtrue\"], y=data[\"yobs\"], func=model, sigma=data[\"yerr\"]\n",
+ ")\n",
+ "result_known_x = bilby.run_sampler(\n",
+ " likelihood=known_x,\n",
+ " label=\"known_x\",\n",
+ " **sampler_kwargs,\n",
+ ")"
]
},
{
@@ -108,16 +142,18 @@
"metadata": {},
"outputs": [],
"source": [
- "result_1D_xtrue.plot_corner(truth=dict(m=5, c = 10), titles = True)\n",
- "result_1D_xtrue.plot_with_data(model = model, x = data['xtrue'], y = data['yobs'], ndraws=1000, npoints=100)\n",
- "plt.show()"
+ "_ = result_known_x.plot_corner(truth=dict(m=5, c=10), titles=True, save=False)\n",
+ "plt.show()\n",
+ "plt.close()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "### As expected this is easy to recover and the sampler does a good job. However this was made too easy - by passing in the 'true' values of x. Lets see what happens when we pass in the observed values of x"
+ "### Fit with unmodeled uncertainty in the x-values\n",
+ "\n",
+ "As expected this is easy to recover and the sampler does a good job. However this was made too easy - by passing in the 'true' values of x. Lets see what happens when we pass in the observed values of x"
]
},
{
@@ -126,11 +162,14 @@
"metadata": {},
"outputs": [],
"source": [
- "OneDGaussian_unknownx = bilby.core.likelihood.GaussianLikelihood(x = data['xobs'], y = data['yobs'], \n",
- " func = model, sigma = data['yerr'])\n",
- "result_1D_xobs = bilby.run_sampler(\n",
- " likelihood=OneDGaussian_unknownx, priors=priors, sampler='dynesty', npoints=livepoints,\n",
- " walks=walks, outdir=outdir, label='xobs_1D_Gaussian')"
+ "incorrect_x = bilby.core.likelihood.GaussianLikelihood(\n",
+ " x=data[\"xobs\"], y=data[\"yobs\"], func=model, sigma=data[\"yerr\"]\n",
+ ")\n",
+ "result_incorrect_x = bilby.run_sampler(\n",
+ " likelihood=incorrect_x,\n",
+ " label=\"incorrect_x\",\n",
+ " **sampler_kwargs,\n",
+ ")"
]
},
{
@@ -139,16 +178,23 @@
"metadata": {},
"outputs": [],
"source": [
- "result_1D_xobs.plot_corner(truth=dict(m=5, c = 10), titles = True)\n",
- "result_1D_xobs.plot_with_data(model = model, x = data['xobs'], y = data['yobs'], ndraws=1000, npoints=100)\n",
- "plt.show()"
+ "_ = result_incorrect_x.plot_corner(truth=dict(m=5, c=10), titles=True, save=False)\n",
+ "plt.show()\n",
+ "plt.close()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "### As expected, this is significantly worse. Let us now define a new likelihood which takes into account x errors but you also pass in xtrue"
+ "### Fit with modeled uncertainty in x-values\n",
+ "\n",
+ "This is not good as there is unmodelled uncertainty in our `x` values.\n",
+ "Getting around this requires marginalisation of the true x values or sampling over them. \n",
+ "See discussion in section 7 of https://arxiv.org/pdf/1008.4686.pdf.\n",
+ "\n",
+ "For this, we will have to define a new likelihood class.\n",
+ "By subclassing the base `bilby.core.likelihood.Likelihood` class we can do this fairly simply."
]
},
{
@@ -157,14 +203,12 @@
"metadata": {},
"outputs": [],
"source": [
- "class TwoDGaussianLikelihood_knownxtrue(bilby.Likelihood):\n",
- " def __init__(self, xtrue, xobs, yobs, xerr, yerr, function):\n",
+ "class GaussianLikelihoodUncertainX(bilby.core.likelihood.Likelihood):\n",
+ " def __init__(self, xobs, yobs, xerr, yerr, function):\n",
" \"\"\"\n",
"\n",
" Parameters\n",
" ----------\n",
- " xtrue: array_like \n",
- " The true injected x values\n",
" xobs, yobs: array_like\n",
" The data to analyse\n",
" xerr, yerr: array_like\n",
@@ -172,21 +216,21 @@
" function:\n",
" The python function to fit to the data\n",
" \"\"\"\n",
+ " super(GaussianLikelihoodUncertainX, self).__init__(dict())\n",
" self.xobs = xobs\n",
- " self.xtrue = xtrue\n",
" self.yobs = yobs\n",
" self.yerr = yerr\n",
" self.xerr = xerr\n",
" self.function = function\n",
- " parameters = inspect.getargspec(function).args\n",
- " parameters.pop(0)\n",
- " self.parameters = dict.fromkeys(parameters)\n",
- " self._marginalized_parameters = list()\n",
"\n",
" def log_likelihood(self):\n",
- " resy = self.yobs - self.function(self.xtrue, **self.parameters)\n",
- " resx = self.xobs - self.xtrue\n",
- " return -0.5 * (np.sum(((resy) / self.yerr) ** 2) + np.sum(((resx) / self.xerr) ** 2))"
+ " variance = (self.xerr * self.parameters[\"m\"]) ** 2 + self.yerr**2\n",
+ " model_y = self.function(self.xobs, **self.parameters)\n",
+ " residual = self.yobs - model_y\n",
+ "\n",
+ " ll = -0.5 * np.sum(residual**2 / variance + np.log(variance))\n",
+ "\n",
+ " return -0.5 * np.sum(residual**2 / variance + np.log(variance))"
]
},
{
@@ -195,12 +239,18 @@
"metadata": {},
"outputs": [],
"source": [
- "TwoDGaussian_knownx = TwoDGaussianLikelihood_knownxtrue(xtrue = data['xtrue'], xobs = data['xobs'], \n",
- " yobs = data['yobs'], xerr=data['xerr'], \n",
- " yerr = data['yerr'], function=model)\n",
- "result_2D_knownx = bilby.run_sampler(\n",
- " likelihood=TwoDGaussian_knownx, priors=priors, sampler='dynesty', npoints=livepoints,\n",
- " walks=walks, outdir=outdir, label='knownx_2D_Gaussian')"
+ "gaussian_unknown_x = GaussianLikelihoodUncertainX(\n",
+ " xobs=data[\"xobs\"],\n",
+ " yobs=data[\"yobs\"],\n",
+ " xerr=data[\"xerr\"],\n",
+ " yerr=data[\"yerr\"],\n",
+ " function=model,\n",
+ ")\n",
+ "result_unknown_x = bilby.run_sampler(\n",
+ " likelihood=gaussian_unknown_x,\n",
+ " label=\"unknown_x\",\n",
+ " **sampler_kwargs,\n",
+ ")"
]
},
{
@@ -209,87 +259,20 @@
"metadata": {},
"outputs": [],
"source": [
- "result_2D_knownx.plot_corner(truth=dict(m=5, c = 10), titles = True)\n",
- "result_2D_knownx.plot_with_data(model = model, x = data['xobs'], y = data['yobs'], ndraws=1000, npoints=100)\n",
- "plt.show()"
+ "_ = result_unknown_x.plot_corner(truth=dict(m=5, c=10), titles=True, save=False)\n",
+ "plt.show()\n",
+ "plt.close()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "### This works well, however it still is not realistic as one still needs to 'know' the true x values. Getting around this requires marginalisation of the true x values or sampling over them. See discussion in section 7 of https://arxiv.org/pdf/1008.4686.pdf"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "class TwoDGaussianLikelihood_unknownx(bilby.Likelihood):\n",
- " def __init__(self, xobs, yobs, xerr, yerr, function):\n",
- " \"\"\"\n",
- "\n",
- " Parameters\n",
- " ----------\n",
- " xobs, yobs: array_like\n",
- " The data to analyse\n",
- " xerr, yerr: array_like\n",
- " The standard deviation of the noise\n",
- " function:\n",
- " The python function to fit to the data\n",
- " \"\"\"\n",
- " self.xobs = xobs\n",
- " self.yobs = yobs\n",
- " self.yerr = yerr\n",
- " self.xerr = xerr\n",
- " self.function = function\n",
- " parameters = inspect.getargspec(function).args\n",
- " parameters.pop(0)\n",
- " self.parameters = dict.fromkeys(parameters)\n",
- " self._marginalized_parameters = list()\n",
- "\n",
- " def log_likelihood(self):\n",
- " m = self.parameters['m']\n",
- " v = np.array([-m, 1.0])\n",
- " \n",
- " Sigma2 = (self.xerr*m)**2 + self.yerr**2\n",
- " model_y = self.function(self.xobs, **self.parameters)\n",
- " Delta = self.yobs - model_y\n",
- "\n",
- " ll = -0.5 * np.sum(Delta**2 / Sigma2 + np.log(Sigma2))\n",
- " \n",
- " return ll"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "TwoDGaussian_unknownx = TwoDGaussianLikelihood_unknownx(xobs = data['xobs'], yobs = data['yobs'], \n",
- " xerr= data['xerr'], yerr = data['yerr'],\n",
- " function=model)\n",
- "result_2D_unknownx = bilby.run_sampler(\n",
- " likelihood=TwoDGaussian_unknownx, priors=priors, sampler='dynesty', npoints=livepoints,\n",
- " walks=walks, outdir=outdir, label='unknownx_2D_Gaussian')"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "result_2D_unknownx.plot_corner(truth=dict(m=5, c = 10), titles = True)\n",
- "result_2D_unknownx.plot_with_data(model = model, x = data['xobs'], y = data['yobs'], ndraws=1000, npoints=100)\n",
- "plt.show()"
+ "Success! The inferred posterior is consistent with the true values."
]
}
],
"metadata": {},
"nbformat": 4,
- "nbformat_minor": 2
+ "nbformat_minor": 4
}
diff --git a/examples/tutorials/making_priors.ipynb b/examples/tutorials/making_priors.ipynb
index 8daef3bc1..a16f941c9 100644
--- a/examples/tutorials/making_priors.ipynb
+++ b/examples/tutorials/making_priors.ipynb
@@ -26,19 +26,14 @@
{
"cell_type": "code",
"execution_count": null,
- "metadata": {
- "execution": {
- "iopub.execute_input": "2021-02-05T21:53:59.048160Z",
- "iopub.status.busy": "2021-02-05T21:53:59.047541Z",
- "iopub.status.idle": "2021-02-05T21:54:00.520314Z",
- "shell.execute_reply": "2021-02-05T21:54:00.520637Z"
- }
- },
+ "metadata": {},
"outputs": [],
"source": [
"import bilby\n",
- "%pylab inline\n",
- "%config InlineBackend.figure_format = 'retina'"
+ "import matplotlib.pyplot as plt\n",
+ "import numpy as np\n",
+ "\n",
+ "%matplotlib inline"
]
},
{
@@ -55,46 +50,50 @@
{
"cell_type": "code",
"execution_count": null,
- "metadata": {
- "execution": {
- "iopub.execute_input": "2021-02-05T21:54:00.528510Z",
- "iopub.status.busy": "2021-02-05T21:54:00.528063Z",
- "iopub.status.idle": "2021-02-05T21:54:02.625157Z",
- "shell.execute_reply": "2021-02-05T21:54:02.624711Z"
- }
- },
+ "metadata": {},
"outputs": [],
"source": [
- "fig = figure(figsize=(12, 5))\n",
+ "fig = plt.figure(figsize=(12, 5))\n",
"\n",
"priors = [\n",
" bilby.core.prior.Uniform(minimum=5, maximum=50),\n",
" bilby.core.prior.LogUniform(minimum=5, maximum=50),\n",
- " bilby.core.prior.PowerLaw(name='name', alpha=2, minimum=100, maximum=1000),\n",
- " bilby.gw.prior.UniformComovingVolume(name='luminosity_distance', minimum=100, maximum=1000, latex_label='label'),\n",
+ " bilby.core.prior.PowerLaw(name=\"name\", alpha=2, minimum=100, maximum=1000),\n",
+ " bilby.gw.prior.UniformComovingVolume(\n",
+ " name=\"luminosity_distance\", minimum=100, maximum=1000, latex_label=\"label\"\n",
+ " ),\n",
" bilby.gw.prior.AlignedSpin(),\n",
- " bilby.core.prior.Gaussian(name='name', mu=0, sigma=1, latex_label='label'),\n",
- " bilby.core.prior.TruncatedGaussian(name='name', mu=1, sigma=0.4, minimum=-1, maximum=1,\n",
- " latex_label='label'),\n",
- " bilby.core.prior.Cosine(name='name', latex_label='label'),\n",
- " bilby.core.prior.Sine(name='name', latex_label='label'),\n",
- " bilby.core.prior.Interped(name='name', xx = np.linspace(0, 10, 1000), yy=np.linspace(0, 10, 1000)**4,\n",
- " minimum=3, maximum=5, latex_label='label')\n",
+ " bilby.core.prior.Gaussian(name=\"name\", mu=0, sigma=1, latex_label=\"label\"),\n",
+ " bilby.core.prior.TruncatedGaussian(\n",
+ " name=\"name\", mu=1, sigma=0.4, minimum=-1, maximum=1, latex_label=\"label\"\n",
+ " ),\n",
+ " bilby.core.prior.Cosine(name=\"name\", latex_label=\"label\"),\n",
+ " bilby.core.prior.Sine(name=\"name\", latex_label=\"label\"),\n",
+ " bilby.core.prior.Interped(\n",
+ " name=\"name\",\n",
+ " xx=np.linspace(0, 10, 1000),\n",
+ " yy=np.linspace(0, 10, 1000) ** 4,\n",
+ " minimum=3,\n",
+ " maximum=5,\n",
+ " latex_label=\"label\",\n",
+ " ),\n",
"]\n",
"\n",
"for ii, prior in enumerate(priors):\n",
" fig.add_subplot(2, 5, 1 + ii)\n",
- " hist(prior.sample(100000), bins=100, histtype='step', density=True)\n",
+ " plt.hist(prior.sample(100000), bins=100, histtype=\"step\", density=True)\n",
" if not isinstance(prior, bilby.core.prior.Gaussian):\n",
- " plot(np.linspace(prior.minimum, prior.maximum, 1000),\n",
- " prior.prob(np.linspace(prior.minimum, prior.maximum, 1000)))\n",
+ " plt.plot(\n",
+ " np.linspace(prior.minimum, prior.maximum, 1000),\n",
+ " prior.prob(np.linspace(prior.minimum, prior.maximum, 1000)),\n",
+ " )\n",
" else:\n",
- " plot(np.linspace(-5, 5, 1000), prior.prob(np.linspace(-5, 5, 1000)))\n",
- " xlabel('{}'.format(prior.latex_label))\n",
- " \n",
- "tight_layout()\n",
- "show()\n",
- "close()"
+ " plt.plot(np.linspace(-5, 5, 1000), prior.prob(np.linspace(-5, 5, 1000)))\n",
+ " plt.xlabel(\"{}\".format(prior.latex_label))\n",
+ "\n",
+ "plt.tight_layout()\n",
+ "plt.show()\n",
+ "plt.close()"
]
},
{
@@ -109,65 +108,59 @@
{
"cell_type": "code",
"execution_count": null,
- "metadata": {
- "execution": {
- "iopub.execute_input": "2021-02-05T21:54:02.634470Z",
- "iopub.status.busy": "2021-02-05T21:54:02.634040Z",
- "iopub.status.idle": "2021-02-05T21:54:02.635882Z",
- "shell.execute_reply": "2021-02-05T21:54:02.635491Z"
- }
- },
+ "metadata": {},
"outputs": [],
"source": [
"class Exponential(bilby.core.prior.Prior):\n",
" \"\"\"Define a new prior class where p(x) ~ exp(alpha * x)\"\"\"\n",
- " \n",
+ "\n",
" def __init__(self, alpha, minimum, maximum, name=None, latex_label=None):\n",
- " super(Exponential, self).__init__(name=name, latex_label=latex_label, minimum=minimum, maximum=maximum)\n",
+ " super(Exponential, self).__init__(\n",
+ " name=name, latex_label=latex_label, minimum=minimum, maximum=maximum\n",
+ " )\n",
" self.alpha = alpha\n",
- " \n",
+ "\n",
" def rescale(self, val):\n",
- " return np.log((np.exp(self.alpha * self.maximum) - np.exp(self.alpha * self.minimum)) * val\n",
- " + np.exp(self.alpha * self.minimum)) / self.alpha\n",
- " \n",
+ " return (\n",
+ " np.log(\n",
+ " (np.exp(self.alpha * self.maximum) - np.exp(self.alpha * self.minimum))\n",
+ " * val\n",
+ " + np.exp(self.alpha * self.minimum)\n",
+ " )\n",
+ " / self.alpha\n",
+ " )\n",
+ "\n",
" def prob(self, val):\n",
" in_prior = (val >= self.minimum) & (val <= self.maximum)\n",
- " return self.alpha * np.exp(self.alpha * val) / (np.exp(self.alpha * self.maximum)\n",
- " - np.exp(self.alpha * self.minimum)) * in_prior"
+ " return (\n",
+ " self.alpha\n",
+ " * np.exp(self.alpha * val)\n",
+ " / (np.exp(self.alpha * self.maximum) - np.exp(self.alpha * self.minimum))\n",
+ " * in_prior\n",
+ " )"
]
},
{
"cell_type": "code",
"execution_count": null,
- "metadata": {
- "execution": {
- "iopub.execute_input": "2021-02-05T21:54:02.640725Z",
- "iopub.status.busy": "2021-02-05T21:54:02.640231Z",
- "iopub.status.idle": "2021-02-05T21:54:03.064456Z",
- "shell.execute_reply": "2021-02-05T21:54:03.063925Z"
- }
- },
+ "metadata": {},
"outputs": [],
"source": [
- "prior = Exponential(name='name', alpha=-1, minimum=0, maximum=10)\n",
- "\n",
- "figure(figsize=(12, 5))\n",
- "hist(prior.sample(100000), bins=100, histtype='step', density=True)\n",
- "plot(np.linspace(prior.minimum, prior.maximum, 1000), prior.prob(np.linspace(prior.minimum, prior.maximum, 1000)))\n",
- "xlabel('{}'.format(prior.latex_label))\n",
- "show()\n",
- "close()"
+ "prior = Exponential(name=\"name\", alpha=-1, minimum=0, maximum=10)\n",
+ "\n",
+ "plt.figure(figsize=(12, 5))\n",
+ "plt.hist(prior.sample(100000), bins=100, histtype=\"step\", density=True)\n",
+ "plt.plot(\n",
+ " np.linspace(prior.minimum, prior.maximum, 1000),\n",
+ " prior.prob(np.linspace(prior.minimum, prior.maximum, 1000)),\n",
+ ")\n",
+ "plt.xlabel(\"{}\".format(prior.latex_label))\n",
+ "plt.show()\n",
+ "plt.close()"
]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": []
}
],
"metadata": {},
"nbformat": 4,
- "nbformat_minor": 2
+ "nbformat_minor": 4
}
diff --git a/examples/tutorials/visualising_the_results.ipynb b/examples/tutorials/visualising_the_results.ipynb
index d839005f3..9f0aec59a 100644
--- a/examples/tutorials/visualising_the_results.ipynb
+++ b/examples/tutorials/visualising_the_results.ipynb
@@ -26,67 +26,124 @@
"source": [
"import bilby\n",
"import matplotlib.pyplot as plt\n",
- "%matplotlib inline\n",
"\n",
- "time_duration = 4. # time duration (seconds)\n",
- "sampling_frequency = 2048. # sampling frequency (Hz)\n",
- "outdir = 'visualising_the_results' # directory in which to store output\n",
- "label = 'example' # identifier to apply to output files\n",
+ "%matplotlib inline"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "duration = 4\n",
+ "sampling_frequency = 2048\n",
+ "outdir = \"visualising_the_results\"\n",
+ "label = \"example\"\n",
"\n",
- "# specify injection parameters\n",
"injection_parameters = dict(\n",
- "mass_1=36., # detector frame (redshifted) primary mass (solar masses)\n",
- "mass_2=29., # detector frame (redshifted) secondary mass (solar masses)\n",
- "a_1=0.4, # primary dimensionless spin magnitude\n",
- "a_2=0.3, # secondary dimensionless spin magnitude\n",
- "tilt_1=0.5, # polar angle between primary spin and the orbital angular momentum (radians)\n",
- "tilt_2=1.0, # polar angle between secondary spin and the orbital angular momentum \n",
- "phi_12=1.7, # azimuthal angle between primary and secondary spin (radians)\n",
- "phi_jl=0.3, # azimuthal angle between total angular momentum and orbital angular momentum (radians)\n",
- "luminosity_distance=200., # luminosity distance to source (Mpc)\n",
- "theta_jn=0.4, # inclination angle between line of sight and orbital angular momentum (radians)\n",
- "phase=1.3, # phase (radians)\n",
- "ra=1.375, # source right ascension (radians)\n",
- "dec=-1.2108, # source declination (radians)\n",
- "geocent_time=1126259642.413, # reference time at geocentre (time of coalescence or peak amplitude) (GPS seconds)\n",
- "psi=2.659 # gravitational wave polarisation angle\n",
- ")\n",
- "\n",
+ " mass_1=36.0,\n",
+ " mass_2=29.0,\n",
+ " a_1=0.4,\n",
+ " a_2=0.3,\n",
+ " tilt_1=0.5,\n",
+ " tilt_2=1.0,\n",
+ " phi_12=1.7,\n",
+ " phi_jl=0.3,\n",
+ " luminosity_distance=1000.0,\n",
+ " theta_jn=0.4,\n",
+ " phase=1.3,\n",
+ " ra=1.375,\n",
+ " dec=-1.2108,\n",
+ " geocent_time=1126259642.413,\n",
+ " psi=2.659,\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
"# specify waveform arguments\n",
"waveform_arguments = dict(\n",
- "waveform_approximant='IMRPhenomPv2', # waveform approximant name\n",
- "reference_frequency=50., # gravitational waveform reference frequency (Hz)\n",
+ " waveform_approximant=\"IMRPhenomXP\", # waveform approximant name\n",
+ " reference_frequency=50.0, # gravitational waveform reference frequency (Hz)\n",
")\n",
"\n",
"# set up the waveform generator\n",
"waveform_generator = bilby.gw.waveform_generator.WaveformGenerator(\n",
- " sampling_frequency=sampling_frequency, duration=time_duration,\n",
+ " sampling_frequency=sampling_frequency,\n",
+ " duration=duration,\n",
" frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole,\n",
- " parameters=injection_parameters, waveform_arguments=waveform_arguments)\n",
- "# create the frequency domain signal\n",
- "hf_signal = waveform_generator.frequency_domain_strain()\n",
- "\n",
- "# initialise an interferometer based on LIGO Hanford, complete with simulated noise and injected signal\n",
- "IFOs = [bilby.gw.detector.get_interferometer_with_fake_noise_and_injection(\n",
- " 'H1', injection_polarizations=hf_signal, injection_parameters=injection_parameters, duration=time_duration,\n",
- " sampling_frequency=sampling_frequency, outdir=outdir)]\n",
- "\n",
+ " parameters=injection_parameters,\n",
+ " waveform_arguments=waveform_arguments,\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "ifos = bilby.gw.detector.InterferometerList([\"H1\", \"L1\"])\n",
+ "ifos.set_strain_data_from_power_spectral_densities(\n",
+ " duration=duration,\n",
+ " sampling_frequency=sampling_frequency,\n",
+ " start_time=injection_parameters[\"geocent_time\"] - 2,\n",
+ ")\n",
+ "_ = ifos.inject_signal(\n",
+ " waveform_generator=waveform_generator, parameters=injection_parameters\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
"# first, set up all priors to be equal to a delta function at their designated value\n",
"priors = bilby.gw.prior.BBHPriorDict(injection_parameters.copy())\n",
"# then, reset the priors on the masses and luminosity distance to conduct a search over these parameters\n",
- "priors['mass_1'] = bilby.core.prior.Uniform(20, 50, 'mass_1')\n",
- "priors['mass_2'] = bilby.core.prior.Uniform(20, 50, 'mass_2')\n",
- "priors['luminosity_distance'] = bilby.core.prior.Uniform(100, 300, 'luminosity_distance')\n",
- "\n",
+ "priors[\"mass_1\"] = bilby.core.prior.Uniform(25, 40, \"mass_1\")\n",
+ "priors[\"mass_2\"] = bilby.core.prior.Uniform(25, 40, \"mass_2\")\n",
+ "priors[\"luminosity_distance\"] = bilby.core.prior.Uniform(\n",
+ " 400, 2000, \"luminosity_distance\"\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
"# compute the likelihoods\n",
- "likelihood = bilby.gw.likelihood.GravitationalWaveTransient(interferometers=IFOs, waveform_generator=waveform_generator)\n",
- "\n",
- "result = bilby.core.sampler.run_sampler(likelihood=likelihood, priors=priors, sampler='dynesty', npoints=100,\n",
- " injection_parameters=injection_parameters, outdir=outdir, label=label,\n",
- " walks=5)\n",
- "\n",
- "# display the corner plot\n",
- "plt.show()"
+ "likelihood = bilby.gw.likelihood.GravitationalWaveTransient(\n",
+ " interferometers=ifos, waveform_generator=waveform_generator\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "result = bilby.core.sampler.run_sampler(\n",
+ " likelihood=likelihood,\n",
+ " priors=priors,\n",
+ " sampler=\"dynesty\",\n",
+ " npoints=100,\n",
+ " injection_parameters=injection_parameters,\n",
+ " outdir=outdir,\n",
+ " label=label,\n",
+ " walks=5,\n",
+ " nact=2,\n",
+ ")"
]
},
{
@@ -121,7 +178,8 @@
"outputs": [],
"source": [
"result.plot_corner()\n",
- "plt.show()"
+ "plt.show()\n",
+ "plt.close()"
]
},
{
@@ -135,12 +193,9 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "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\n",
- "[chain consumer](https://samreay.github.io/ChainConsumer/index.html), and all the keyword arguments passed to `plot_corner` are passed through to [the `plot` function of chain consumer](https://samreay.github.io/ChainConsumer/chain_api.html#chainconsumer.plotter.Plotter.plot).\n",
- "\n",
- "### Adding injection parameters to the plot\n",
- "\n",
- "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."
+ "## Waveform Reconstruction plot\n",
+ "Some plots specific to compact binary coalescence parameter estimation results can\n",
+ "be created by re-loading the result as a `CBCResult`:"
]
},
{
@@ -149,17 +204,24 @@
"metadata": {},
"outputs": [],
"source": [
- "result.plot_corner(truth=dict(luminosity_distance=201))\n",
- "plt.show()"
+ "from bilby.gw.result import CBCResult\n",
+ "\n",
+ "cbc_result = CBCResult.from_json(\"visualising_the_results/example_result.json\")\n",
+ "for ifo in ifos:\n",
+ " cbc_result.plot_interferometer_waveform_posterior(\n",
+ " interferometer=ifo, n_samples=500, save=False\n",
+ " )\n",
+ " plt.show()\n",
+ " plt.close()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "### Plot a subset of the corner plot\n",
+ "### Marginal Distribution plots\n",
"\n",
- "Or, to plot just a subset of parameters, just pass a list of the names you want."
+ "These plots just show the 1D histograms for each parameter"
]
},
{
@@ -168,15 +230,20 @@
"metadata": {},
"outputs": [],
"source": [
- "result.plot_corner(parameters=['mass_1', 'mass_2'], filename='{}/subset.png'.format(outdir))\n",
- "plt.show()"
+ "result.plot_marginals()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "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"
+ "## Customizing corner plots\n",
+ "\n",
+ "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](https://github.com/dfm/corner.py), and all the keyword arguments passed to `plot_corner` are passed through.\n",
+ "\n",
+ "### Adding injection parameters to the plot\n",
+ "\n",
+ "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."
]
},
{
@@ -185,16 +252,18 @@
"metadata": {},
"outputs": [],
"source": [
- "!ls visualising_the_results/"
+ "result.plot_corner(truth=dict(luminosity_distance=201))\n",
+ "plt.show()\n",
+ "plt.close()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "## Alternative\n",
+ "### Plot a subset of the corner plot\n",
"\n",
- "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`](http://corner.readthedocs.io/en/latest/) package"
+ "Or, to plot just a subset of parameters, just pass a list of the names you want."
]
},
{
@@ -203,24 +272,18 @@
"metadata": {},
"outputs": [],
"source": [
- "import corner\n",
- "samples = result.samples\n",
- "labels = result.parameter_labels\n",
- "fig = corner.corner(samples, labels=labels)\n",
- "plt.show()"
+ "result.plot_corner(\n",
+ " parameters=[\"mass_1\", \"mass_2\"], filename=\"{}/subset.png\".format(outdir)\n",
+ ")\n",
+ "plt.show()\n",
+ "plt.close()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "## Other plots\n",
- "\n",
- "We also include some other types of plots which may be useful. Again, these are built on chain consumer so you may find it useful to check the [documentation](https://samreay.github.io/ChainConsumer/chain_api.html#plotter-class) to see how these plots can be extended. Below, we show just one example of these.\n",
- "\n",
- "#### Distribution plots\n",
- "\n",
- "These plots just show the 1D histograms for each parameter"
+ "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"
]
},
{
@@ -229,70 +292,35 @@
"metadata": {},
"outputs": [],
"source": [
- "result.plot_marginals()\n",
- "plt.show()"
+ "!ls visualising_the_results/"
]
},
{
"cell_type": "markdown",
- "metadata": {
- "pycharm": {
- "name": "#%% md\n"
- }
- },
+ "metadata": {},
"source": [
- "#### Best-Fit Time Domain Waveform plot\n",
- "Some plots specific to compact binary coalescence parameter estimation results can\n",
- "be created by re-loading the result as a `CBCResult`:"
+ "## Alternative\n",
+ "\n",
+ "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](http://corner.readthedocs.io/en/latest/) package"
]
},
{
"cell_type": "code",
"execution_count": null,
- "metadata": {
- "pycharm": {
- "name": "#%%\n"
- }
- },
- "outputs": [],
- "source": [
- "from bilby.gw.result import CBCResult\n",
- "\n",
- "cbc_result = CBCResult.from_json(\"visualising_the_results/example_result.json\")\n",
- "cbc_result.plot_waveform_posterior()\n",
- "plt.show()"
- ]
- },
- {
- "cell_type": "markdown",
"metadata": {},
+ "outputs": [],
"source": [
- "Again, notice that the plot is saved as a \"waveform.png\" in the output dir.\n",
- "\n",
+ "import corner\n",
"\n",
- "\n"
+ "samples = result.samples\n",
+ "labels = result.parameter_labels\n",
+ "fig = corner.corner(samples, labels=labels)\n",
+ "plt.show()\n",
+ "plt.close()"
]
}
],
- "metadata": {
- "kernelspec": {
- "display_name": "Python 3",
- "language": "python",
- "name": "python3"
- },
- "language_info": {
- "codemirror_mode": {
- "name": "ipython",
- "version": 3
- },
- "file_extension": ".py",
- "mimetype": "text/x-python",
- "name": "python",
- "nbconvert_exporter": "python",
- "pygments_lexer": "ipython3",
- "version": "3.8.8"
- }
- },
+ "metadata": {},
"nbformat": 4,
- "nbformat_minor": 2
+ "nbformat_minor": 4
}
diff --git a/test/core/sampler/dynesty_test.py b/test/core/sampler/dynesty_test.py
index d1c42b099..8cbab6d31 100644
--- a/test/core/sampler/dynesty_test.py
+++ b/test/core/sampler/dynesty_test.py
@@ -33,7 +33,6 @@ class TestDynesty(unittest.TestCase):
sample="rwalk",
periodic=None,
reflective=None,
- verbose=True,
check_point_delta_t=1800,
nlive=1000,
first_update=None,
@@ -90,7 +89,6 @@ class TestDynesty(unittest.TestCase):
sample="rwalk",
periodic=[],
reflective=[],
- verbose=True,
check_point_delta_t=1800,
nlive=1000,
first_update=None,
diff --git a/test/gw/utils_test.py b/test/gw/utils_test.py
index d8c3286ac..c5f3af46b 100644
--- a/test/gw/utils_test.py
+++ b/test/gw/utils_test.py
@@ -101,12 +101,16 @@ class TestGWUtils(unittest.TestCase):
self.assertEqual(mfsnr, 25.510869054168282)
def test_get_event_time(self):
+ from urllib3.exceptions import NewConnectionError
events = [
"GW150914",
"GW170104",
]
for event in events:
- self.assertTrue(isinstance(gwutils.get_event_time(event), float))
+ try:
+ self.assertTrue(isinstance(gwutils.get_event_time(event), float))
+ except NewConnectionError:
+ return
with self.assertRaises(ValueError):
gwutils.get_event_time("GW010290")
diff --git a/test/integration/sampler_run_test.py b/test/integration/sampler_run_test.py
index 0c4010720..2bf1d355e 100644
--- a/test/integration/sampler_run_test.py
+++ b/test/integration/sampler_run_test.py
@@ -27,6 +27,7 @@ class TestRunningSamplers(unittest.TestCase):
self.kwargs = dict(
save=False,
conversion_function=self.conversion_function,
+ verbose=True,
)
bilby.core.utils.check_directory_exists_and_if_not_mkdir("outdir")
--
GitLab