From e636e4f6f3127f6d032a910ce23c8fee30e21b45 Mon Sep 17 00:00:00 2001
From: Moritz <email@moritz-huebner.de>
Date: Thu, 24 Oct 2019 11:32:17 +1100
Subject: [PATCH] Created a generic conditional prior example
---
examples/core_examples/conditional_prior.py | 11 ++-
.../injection_examples/conditional_prior.py | 90 -------------------
2 files changed, 7 insertions(+), 94 deletions(-)
delete mode 100644 examples/gw_examples/injection_examples/conditional_prior.py
diff --git a/examples/core_examples/conditional_prior.py b/examples/core_examples/conditional_prior.py
index 3e4679983..2475e479d 100644
--- a/examples/core_examples/conditional_prior.py
+++ b/examples/core_examples/conditional_prior.py
@@ -3,31 +3,33 @@ import numpy as np
# This tutorial demonstrates how we can sample a prior in the shape of a ball
# Note that this will not end up sampling uniformly in that space, constraint priors are more suitable for that.
-# This implementation will draw a value for the x-coordinate, and given that draw a value for the
-# y-coordinate, and given that draw a value for the z-coordinate.
+# This implementation will draw a value for the x-coordinate from p(x), and given that draw a value for the
+# y-coordinate from p(y|x), and given that draw a value for the z-coordinate from p(z|x,y).
# Only the x-coordinate will end up being uniform for this
class ZeroLikelihood(bilby.core.likelihood.Likelihood):
""" Flat likelihood. This always returns 0.
- This way we can see if we correctly sampled uniformly in the prior space"""
+ This way our posterior distribution is exactly the prior distribution."""
def log_likelihood(self):
return 0
def condition_func_y(reference_params, x):
+ """ Condition function for our p(y|x) prior."""
radius = 0.5 * (reference_params['maximum'] - reference_params['minimum'])
y_max = np.sqrt(radius**2 - x**2)
return dict(minimum=-y_max, maximum=y_max)
def condition_func_z(reference_params, x, y):
- """"""
+ """ Condition function for our p(z|x, y) prior."""
radius = 0.5 * (reference_params['maximum'] - reference_params['minimum'])
z_max = np.sqrt(radius**2 - x**2 - y**2)
return dict(minimum=-z_max, maximum=z_max)
+# Set up the conditional priors and the flat likelihood
priors = bilby.core.prior.ConditionalPriorDict()
priors['x'] = bilby.core.prior.Uniform(minimum=-1, maximum=1, latex_label="$x$")
priors['y'] = bilby.core.prior.ConditionalUniform(condition_func=condition_func_y, minimum=-1,
@@ -36,6 +38,7 @@ priors['z'] = bilby.core.prior.ConditionalUniform(condition_func=condition_func_
maximum=1, latex_label="$z$")
likelihood = ZeroLikelihood(parameters=dict(x=0, y=0, z=0))
+# Sample the prior distribution
res = bilby.run_sampler(likelihood=likelihood, priors=priors, sampler='dynesty', npoints=5000, walks=100,
label='conditional_prior', outdir='outdir', resume=False, clean=True)
res.plot_corner()
diff --git a/examples/gw_examples/injection_examples/conditional_prior.py b/examples/gw_examples/injection_examples/conditional_prior.py
deleted file mode 100644
index 4774c9aca..000000000
--- a/examples/gw_examples/injection_examples/conditional_prior.py
+++ /dev/null
@@ -1,90 +0,0 @@
-import matplotlib.pyplot as plt
-import numpy as np
-
-import bilby.gw.prior
-
-
-def condition_function(reference_params, mass_1):
- return dict(minimum=reference_params['minimum'], maximum=mass_1)
-
-
-mass_1_min = 5
-mass_1_max = 100
-
-mass_1 = bilby.core.prior.Uniform(minimum=mass_1_min, maximum=mass_1_max, name='mass_1',
- latex_label='$m_1$', boundary='reflective')
-mass_2 = bilby.core.prior.ConditionalUniform(minimum=mass_1_min, maximum=mass_1_max, name='mass_2',
- latex_label='$m_2$', condition_func=condition_function,
- boundary='reflective')
-
-conditional_dict = bilby.core.prior.ConditionalPriorDict(dictionary=dict(mass_1=mass_1, mass_2=mass_2))
-
-
-res = conditional_dict.sample(100000)
-
-plt.hist(res['mass_1'], bins='fd', alpha=0.6, density=True, label='Sampled')
-plt.plot(np.linspace(2, 100, 200), conditional_dict['mass_1'].prob(np.linspace(2, 100, 200)), label='Uniform prior')
-plt.xlabel('$m_1$')
-plt.ylabel('$p(m_1)$')
-plt.legend()
-plt.tight_layout()
-plt.show()
-plt.clf()
-
-
-plt.hist(res['mass_2'], bins='fd', alpha=0.6, density=True, label='Sampled')
-plt.xlabel('$m_2$')
-plt.ylabel('$p(m_2 | m_1)$')
-plt.legend()
-plt.tight_layout()
-plt.show()
-plt.clf()
-
-
-duration = 4.
-sampling_frequency = 2048.
-outdir = 'outdir'
-label = 'conditional_prior'
-bilby.core.utils.setup_logger(outdir=outdir, label=label)
-
-np.random.seed(88170235)
-
-injection_parameters = dict(
- mass_1=30., mass_2=30, 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=1500., theta_jn=0.4, psi=2.659,
- phase=1.3, geocent_time=1126259642.413, ra=1.375, dec=-1.2108)
-
-waveform_arguments = dict(waveform_approximant='IMRPhenomPv2',
- reference_frequency=50., minimum_frequency=20.)
-
-waveform_generator = bilby.gw.WaveformGenerator(
- duration=duration, sampling_frequency=sampling_frequency,
- frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole,
- waveform_arguments=waveform_arguments)
-
-ifos = bilby.gw.detector.InterferometerList(['H1', 'L1'])
-ifos.set_strain_data_from_power_spectral_densities(
- sampling_frequency=sampling_frequency, duration=duration,
- start_time=injection_parameters['geocent_time'] - 3)
-ifos.inject_signal(waveform_generator=waveform_generator,
- parameters=injection_parameters)
-
-priors = bilby.core.prior.ConditionalPriorDict()
-for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'phi_12', 'phi_jl', 'psi', 'ra',
- 'dec', 'geocent_time', 'phase', 'theta_jn', 'luminosity_distance']:
- priors[key] = injection_parameters[key]
-priors['mass_1'] = mass_1
-priors['mass_2'] = mass_2
-
-# Initialise the likelihood by passing in the interferometer data (ifos) and
-# the waveform generator
-likelihood = bilby.gw.GravitationalWaveTransient(
- interferometers=ifos, waveform_generator=waveform_generator)
-
-# Run sampler. In this case we're going to use the `dynesty` sampler
-result = bilby.run_sampler(
- likelihood=likelihood, priors=priors, sampler='dynesty', npoints=100,
- injection_parameters=injection_parameters, outdir=outdir, label=label, clean=True, resume=False)
-
-# Make a corner plot.
-result.plot_corner()
--
GitLab