diff --git a/.gitignore b/.gitignore
index 48e4c01a1387a87c5af4a7381f1f5d157171e3c0..437f5654e8ab5c75ef9e57427ae65d39c467b37c 100644
--- a/.gitignore
+++ b/.gitignore
@@ -13,6 +13,6 @@ MANIFEST
*.dat
*.version
*.ipynb_checkpoints
-outdir/*
+**/outdir
.idea/*
bilby/_version.py
diff --git a/bilby/core/sampler/dynesty.py b/bilby/core/sampler/dynesty.py
index d9c6ab6c01180226ce3bb91e58e0928df9cc00d7..a4e459f59089b723fc9a0747dd6c1aa49edce10b 100644
--- a/bilby/core/sampler/dynesty.py
+++ b/bilby/core/sampler/dynesty.py
@@ -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
diff --git a/examples/core_examples/conditional_prior.py b/examples/core_examples/conditional_prior.py
deleted file mode 100644
index 9dc236a87c5c2986c9bb5320def58aed1c182c39..0000000000000000000000000000000000000000
--- a/examples/core_examples/conditional_prior.py
+++ /dev/null
@@ -1,59 +0,0 @@
-"""
-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()
diff --git a/examples/gw_examples/injection_examples/relative_binning.py b/examples/gw_examples/injection_examples/relative_binning.py
index 5a848f003b0a88262d0380272c85b624ef306111..bc224ca6208a7734b2495c7291830ba4bd640427 100644
--- a/examples/gw_examples/injection_examples/relative_binning.py
+++ b/examples/gw_examples/injection_examples/relative_binning.py
@@ -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",
+)
diff --git a/examples/tutorials/conditional_priors.ipynb b/examples/tutorials/conditional_priors.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..0838e5da9cb205f467fcf58c94acef9aea36e08c
--- /dev/null
+++ b/examples/tutorials/conditional_priors.ipynb
@@ -0,0 +1,229 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Conditional prior demonstration\n",
+ "\n",
+ "Conditional priors enable inference to be performed with priors that correlate different parameters.\n",
+ "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.\n",
+ "\n",
+ "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",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import matplotlib.pyplot as plt\n",
+ "import numpy as np\n",
+ "import pandas as pd\n",
+ "from bilby.core.prior import (\n",
+ " Prior, PriorDict, ConditionalPriorDict,\n",
+ " Uniform, ConditionalUniform, Constraint, \n",
+ ")\n",
+ "from corner import corner\n",
+ "from scipy.stats import semicircular\n",
+ "\n",
+ "\n",
+ "%matplotlib inline"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Sampling from a disc\n",
+ "\n",
+ "Our first example is sampling uniformly from a disc\n",
+ "\n",
+ "$$p(x, y) = \\frac{1}{\\pi}; x^2 + y+2 \\leq 1.$$\n",
+ "\n",
+ "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.\n",
+ "Naive sampling from this parameterization would have an efficiency of $\\pi / 4$.\n",
+ "\n",
+ "If we instead consider the marginal distribution $p(x)$ and conditional distribution $p(y | x)$, we can achieve a sampling efficiency of 100%.\n",
+ "\n",
+ "$$\n",
+ "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} \\\\\n",
+ "p(y | x) = \\frac{1}{2 (1 - \\sqrt{x^2})}\n",
+ "$$\n",
+ "\n",
+ "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.\n",
+ "The conditional distribution for $y$ is implemented in `Bilby` as the `ConditionUnifrom`, we just need to define the condition function."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "class SemiCircular(Prior):\n",
+ "\n",
+ " def __init__(self, radius=1, center=0, name=None, latex_label=None, unit=None, boundary=None):\n",
+ " super(SemiCircular, self).__init__(\n",
+ " minimum=center - radius,\n",
+ " maximum=center + radius,\n",
+ " name=name,\n",
+ " latex_label=latex_label,\n",
+ " unit=unit,\n",
+ " boundary=boundary,\n",
+ " )\n",
+ " self.radius = radius\n",
+ " self.center = center\n",
+ " self._dist = semicircular(loc=center, scale=radius)\n",
+ "\n",
+ " def prob(self, val):\n",
+ " return self._dist.pdf(val)\n",
+ "\n",
+ " def ln_prob(self, val):\n",
+ " return self._dist.logpdf(val)\n",
+ "\n",
+ " def cdf(self, val):\n",
+ " return self._dist.cdf(val)\n",
+ "\n",
+ " def rescale(self, val):\n",
+ " return self._dist.ppf(val)\n",
+ "\n",
+ "\n",
+ "def conditional_func_y(reference_parameters, x):\n",
+ " condition = np.sqrt(reference_parameters[\"maximum\"]-x**2)\n",
+ " return dict(minimum=-condition, maximum=condition)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "#### Sample from the distribution\n",
+ "\n",
+ "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",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "N = int(2e4)\n",
+ "\n",
+ "CORNER_KWARGS = dict(\n",
+ " plot_contours=False,\n",
+ " plot_density=False,\n",
+ " fill_contours=False,\n",
+ " max_n_ticks=3,\n",
+ " verbose=False,\n",
+ " use_math_text=True,\n",
+ ")\n",
+ "\n",
+ "\n",
+ "def convert_to_radial(parameters):\n",
+ " p = parameters.copy()\n",
+ " p['r'] = p['x']**2 + p['y']**2\n",
+ " return p\n",
+ "\n",
+ "def sample_circle_with_constraint():\n",
+ " d = PriorDict(\n",
+ " dictionary=dict(\n",
+ " x=Uniform(-1, 1),\n",
+ " y=Uniform(-1, 1),\n",
+ " r=Constraint(0, 1),\n",
+ " ),\n",
+ " conversion_function=convert_to_radial\n",
+ " )\n",
+ " return pd.DataFrame(d.sample(N))\n",
+ "\n",
+ "\n",
+ "def sample_circle_with_conditional():\n",
+ " d = ConditionalPriorDict(\n",
+ " dictionary=dict(\n",
+ " x=SemiCircular(),\n",
+ " y=ConditionalUniform(\n",
+ " condition_func=conditional_func_y, \n",
+ " minimum=-1, maximum=1\n",
+ " )\n",
+ " )\n",
+ " )\n",
+ " return pd.DataFrame(d.sample(N))\n",
+ "\n",
+ "\n",
+ "s1 = sample_circle_with_constraint()\n",
+ "s2 = sample_circle_with_conditional()\n",
+ "fig = corner(s1.values, **CORNER_KWARGS, color=\"tab:blue\", labels=[\"$x$\", \"$y$\"])\n",
+ "corner(s2.values, **CORNER_KWARGS, color=\"tab:green\", fig=fig)\n",
+ "plt.show()\n",
+ "plt.close()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Sampling from ordered distributions\n",
+ "\n",
+ "As our second example, we demonstrate defining a prior distribution over a set of strictly ordered parameters.\n",
+ "\n",
+ "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",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "class BoundedUniform(ConditionalUniform):\n",
+ " \"\"\"Conditional Uniform prior where prior sample < previous prior sample\n",
+ " \n",
+ " This is ensured by fixing the maximum bound to be the previous prior sample value.\n",
+ " \"\"\"\n",
+ " def __init__(self, idx: int, minimum, maximum, name=None, latex_label=None,\n",
+ " unit=None, boundary=None):\n",
+ " super(BoundedUniform, self).__init__(\n",
+ " minimum=minimum, maximum=maximum, name=name, \n",
+ " latex_label=latex_label, unit=unit,\n",
+ " boundary=boundary, condition_func=self.bounds_condition\n",
+ " )\n",
+ " self.idx = idx\n",
+ " self.previous_name = f\"{name[:-1]}{self.idx - 1}\"\n",
+ " self._required_variables = [self.previous_name] \n",
+ " # this is used in prior.sample(... **required_variables)\n",
+ "\n",
+ "\n",
+ " def bounds_condition(self, reference_params, **required_variables):\n",
+ " previous_sample = required_variables[self.previous_name]\n",
+ " return dict(maximum=previous_sample)\n",
+ "\n",
+ "\n",
+ "def make_uniform_conditonal_priordict(n_priors=3):\n",
+ " priors = ConditionalPriorDict()\n",
+ " for i in range(n_priors):\n",
+ " if i==0:\n",
+ " priors[f\"uni{i}\"] = Uniform(minimum=0, maximum=1, name=f\"uni{i}\")\n",
+ " else:\n",
+ " priors[f\"uni{i}\"] = BoundedUniform(idx=i, minimum=0, maximum=1, name=f\"uni{i}\")\n",
+ " return priors\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "samples = pd.DataFrame(make_uniform_conditonal_priordict(3).sample(10000))\n",
+ "fig = corner(samples.values, **CORNER_KWARGS, color=\"tab:blue\", labels=[f\"$A_{ii}$\" for ii in range(3)])\n",
+ "plt.show()\n",
+ "plt.close()"
+ ]
+ }
+ ],
+ "metadata": {},
+ "nbformat": 4,
+ "nbformat_minor": 2
+}