Skip to content
Snippets Groups Projects
Commit 5ae22fb5 authored by Colm Talbot's avatar Colm Talbot
Browse files

EXAMPLES: fix broken examples

parent e73045f8
No related branches found
No related tags found
1 merge request!1182EXAMPLES: fix broken examples
Hide whitespace changes
Inline Side-by-side
.gitignore
+ 1
− 1
View file @ 5ae22fb5
......@@ -13,6 +13,6 @@ MANIFEST
*.dat
*.version
*.ipynb_checkpoints
outdir/*
**/outdir
.idea/*
bilby/_version.py
bilby/core/sampler/dynesty.py
+ 8
− 1
View file @ 5ae22fb5
......@@ -725,7 +725,7 @@ class Dynesty(NestedSampler):
def _run_test(self):
import pandas as pd
self.sampler = self.sampler_class(
self.sampler = self.sampler_init(
loglikelihood=self.log_likelihood,
prior_transform=self.prior_transform,
ndim=self.ndim,
......@@ -734,7 +734,14 @@ class Dynesty(NestedSampler):
sampler_kwargs = self.sampler_function_kwargs.copy()
sampler_kwargs["maxiter"] = 2
if self.print_method == "tqdm" and self.kwargs["print_progress"]:
from tqdm.auto import tqdm
self.pbar = tqdm(file=sys.stdout, initial=self.sampler.it)
self.sampler.run_nested(**sampler_kwargs)
if self.pbar is not None:
self.pbar = self.pbar.close()
print("")
N = 100
self.result.samples = pd.DataFrame(self.priors.sample(N))[
self.search_parameter_keys
......
examples/core_examples/conditional_prior.py deleted 100644 → 0
+ 0
− 59
View file @ e73045f8
"""
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 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
"""
import bilby
import numpy as np
class ZeroLikelihood(bilby.core.likelihood.Likelihood):
"""Flat likelihood. This always returns 0.
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, maximum=1, latex_label="$y$"
)
priors["z"] = bilby.core.prior.ConditionalUniform(
condition_func=condition_func_z, minimum=-1, 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",
nlive=5000,
label="conditional_prior",
outdir="outdir",
resume=False,
clean=True,
)
res.plot_corner()
examples/gw_examples/injection_examples/relative_binning.py
+ 27
− 5
View file @ 5ae22fb5
......@@ -7,6 +7,7 @@ This example estimates the masses using a uniform prior in both component masses
and distance using a uniform in comoving volume prior on luminosity distance
between luminosity distances of 100Mpc and 5Gpc, the cosmology is Planck15.
"""
from copy import deepcopy
import bilby
import numpy as np
......@@ -109,6 +110,17 @@ priors["fiducial"] = 0
# Perform a check that the prior does not extend to a parameter space longer than the data
priors.validate_prior(duration, minimum_frequency)
# Set up the fiducial parameters for the relative binning likelihood to be the
# injected parameters. Note that because we sample in chirp mass and mass ratio
# but injected with mass_1 and mass_2, we need to convert the mass parameters
fiducial_parameters = injection_parameters.copy()
m1 = fiducial_parameters.pop("mass_1")
m2 = fiducial_parameters.pop("mass_2")
fiducial_parameters["chirp_mass"] = bilby.gw.conversion.component_masses_to_chirp_mass(
m1, m2
)
fiducial_parameters["mass_ratio"] = m2 / m1
# Initialise the likelihood by passing in the interferometer data (ifos) and
# the waveform generator
likelihood = bilby.gw.likelihood.RelativeBinningGravitationalWaveTransient(
......@@ -116,15 +128,15 @@ likelihood = bilby.gw.likelihood.RelativeBinningGravitationalWaveTransient(
waveform_generator=waveform_generator,
priors=priors,
distance_marginalization=True,
fiducial_parameters=injection_parameters,
fiducial_parameters=fiducial_parameters,
)
# Run sampler. In this case we're going to use the `dynesty` sampler
# Run sampler. In this case, we're going to use the `nestle` sampler
result = bilby.run_sampler(
likelihood=likelihood,
priors=priors,
sampler="nestle",
npoints=100,
npoints=1000,
injection_parameters=injection_parameters,
outdir=outdir,
label=label,
......@@ -158,5 +170,15 @@ print(
)
print(f"Binned vs unbinned log Bayes factor {np.log(np.mean(weights)):.2f}")
# Make a corner plot.
# result.plot_corner()
# Generate result object with the posterior for the regular likelihood using
# rejection sampling
alt_result = deepcopy(result)
keep = weights > np.random.uniform(0, max(weights), len(weights))
alt_result.posterior = result.posterior.iloc[keep]
# Make a comparison corner plot.
bilby.core.result.plot_multiple(
[result, alt_result],
labels=["Binned", "Reweighted"],
filename=f"{outdir}/{label}_corner.png",
)
examples/tutorials/conditional_priors.ipynb 0 → 100644
+ 229
− 0
View file @ 5ae22fb5
%% Cell type:markdown id: tags:
# Conditional prior demonstration
Conditional priors enable inference to be performed with priors that correlate different parameters.
In this notebook, we demonstrate two uses of this: maintaining a two-dimensional distribution while changing the parameterization to be more efficient for sampling, and enforcing an ordering between parameters.
Many cases where `Conditional` priors are useful can also be expressed with `Constraint` priors, however the conditional approach can improve sampling efficiency.
%% Cell type:code id: tags:
```
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from bilby.core.prior import (
Prior, PriorDict, ConditionalPriorDict,
Uniform, ConditionalUniform, Constraint,
)
from corner import corner
from scipy.stats import semicircular
%matplotlib inline
```
%% Cell type:markdown id: tags:
### Sampling from a disc
Our first example is sampling uniformly from a disc
$$p(x, y) = \frac{1}{\pi}; x^2 + y+2 \leq 1.$$
A naive implementation of this would define a uniform prior over a square over `[-1, 1]` and then reject points that don't satisfy the radius constraint.
Naive sampling from this parameterization would have an efficiency of $\pi / 4$.
If we instead consider the marginal distribution $p(x)$ and conditional distribution $p(y | x)$, we can achieve a sampling efficiency of 100%.
$$
p(x) = \int_{-1 + \sqrt{x^2 + y^2}}^{1 - \sqrt{x^2 + y^2}} dy p(x, y) = \frac{2 (1 - \sqrt{x^2 + y^2})}{\pi} \\
p(y | x) = \frac{1}{2 (1 - \sqrt{x^2})}
$$
The marginal distribution for $x$ is the [Wigner semicircle distribution](https://en.wikipedia.org/wiki/Wigner_semicircle_distribution), this distribution is not currently defined in `Bilby`, but we can wrap the `scipy` implementation.
The conditional distribution for $y$ is implemented in `Bilby` as the `ConditionUnifrom`, we just need to define the condition function.
%% Cell type:code id: tags:
```
class SemiCircular(Prior):
def __init__(self, radius=1, center=0, name=None, latex_label=None, unit=None, boundary=None):
super(SemiCircular, self).__init__(
minimum=center - radius,
maximum=center + radius,
name=name,
latex_label=latex_label,
unit=unit,
boundary=boundary,
)
self.radius = radius
self.center = center
self._dist = semicircular(loc=center, scale=radius)
def prob(self, val):
return self._dist.pdf(val)
def ln_prob(self, val):
return self._dist.logpdf(val)
def cdf(self, val):
return self._dist.cdf(val)
def rescale(self, val):
return self._dist.ppf(val)
def conditional_func_y(reference_parameters, x):
condition = np.sqrt(reference_parameters["maximum"]-x**2)
return dict(minimum=-condition, maximum=condition)
```
%% Cell type:markdown id: tags:
#### Sample from the distribution
To demonstrate the equivalence of the two methods, we will draw samples from the distribution using the two methods and verify that they agree.
%% Cell type:code id: tags:
```
N = int(2e4)
CORNER_KWARGS = dict(
plot_contours=False,
plot_density=False,
fill_contours=False,
max_n_ticks=3,
verbose=False,
use_math_text=True,
)
def convert_to_radial(parameters):
p = parameters.copy()
p['r'] = p['x']**2 + p['y']**2
return p
def sample_circle_with_constraint():
d = PriorDict(
dictionary=dict(
x=Uniform(-1, 1),
y=Uniform(-1, 1),
r=Constraint(0, 1),
),
conversion_function=convert_to_radial
)
return pd.DataFrame(d.sample(N))
def sample_circle_with_conditional():
d = ConditionalPriorDict(
dictionary=dict(
x=SemiCircular(),
y=ConditionalUniform(
condition_func=conditional_func_y,
minimum=-1, maximum=1
)
)
)
return pd.DataFrame(d.sample(N))
s1 = sample_circle_with_constraint()
s2 = sample_circle_with_conditional()
fig = corner(s1.values, **CORNER_KWARGS, color="tab:blue", labels=["$x$", "$y$"])
corner(s2.values, **CORNER_KWARGS, color="tab:green", fig=fig)
plt.show()
plt.close()
```
%% Cell type:markdown id: tags:
### Sampling from ordered distributions
As our second example, we demonstrate defining a prior distribution over a set of strictly ordered parameters.
We note that in this case, we do not require that the marginal distributions for each of the parameters are independently and identically disributed, although this can be fairly simply remedied.
%% Cell type:code id: tags:
```
class BoundedUniform(ConditionalUniform):
"""Conditional Uniform prior where prior sample < previous prior sample
This is ensured by fixing the maximum bound to be the previous prior sample value.
"""
def __init__(self, idx: int, minimum, maximum, name=None, latex_label=None,
unit=None, boundary=None):
super(BoundedUniform, self).__init__(
minimum=minimum, maximum=maximum, name=name,
latex_label=latex_label, unit=unit,
boundary=boundary, condition_func=self.bounds_condition
)
self.idx = idx
self.previous_name = f"{name[:-1]}{self.idx - 1}"
self._required_variables = [self.previous_name]
# this is used in prior.sample(... **required_variables)
def bounds_condition(self, reference_params, **required_variables):
previous_sample = required_variables[self.previous_name]
return dict(maximum=previous_sample)
def make_uniform_conditonal_priordict(n_priors=3):
priors = ConditionalPriorDict()
for i in range(n_priors):
if i==0:
priors[f"uni{i}"] = Uniform(minimum=0, maximum=1, name=f"uni{i}")
else:
priors[f"uni{i}"] = BoundedUniform(idx=i, minimum=0, maximum=1, name=f"uni{i}")
return priors
```
%% Cell type:code id: tags:
```
samples = pd.DataFrame(make_uniform_conditonal_priordict(3).sample(10000))
fig = corner(samples.values, **CORNER_KWARGS, color="tab:blue", labels=[f"$A_{ii}$" for ii in range(3)])
plt.show()
plt.close()
```
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment