From 075986319422df2a58bf8ad8c992eef089f5f415 Mon Sep 17 00:00:00 2001
From: Matthew Pitkin <matthew.pitkin@ligo.org>
Date: Wed, 13 Mar 2019 12:33:11 +0000
Subject: [PATCH] Update the multivariate Gaussian likelihood prior example
---
.../multivariate_gaussian_prior.py | 51 ++++++++++++-------
sampler_requirements.txt | 2 +-
2 files changed, 34 insertions(+), 19 deletions(-)
diff --git a/examples/other_examples/multivariate_gaussian_prior.py b/examples/other_examples/multivariate_gaussian_prior.py
index e3f258c1..1ccbd994 100644
--- a/examples/other_examples/multivariate_gaussian_prior.py
+++ b/examples/other_examples/multivariate_gaussian_prior.py
@@ -7,6 +7,8 @@ Gaussian prior distribution.
from __future__ import division
import bilby
import numpy as np
+from scipy import linalg, stats
+import matplotlib as mpl
from bilby.core.likelihood import GaussianLikelihood
@@ -34,16 +36,15 @@ N = len(time)
data = model(time, 0., 0.) # noiseless data
-# Now lets instantiate a version of our GaussianLikelihood, giving it
-# the time, data, signal model and standard deviation
+# instantiate the GaussianLikelihood
likelihood = GaussianLikelihood(time, data, model, sigma=sigma)
# Create a Multivariate Gaussian prior distribution with two modes
names = ['m', 'c']
-mus = [[-5., -5.], [5., 5.]] # means of two modes
-corrcoefs = [[[1., -0.7], [-0.7, 1.]], [[1., -0.7], [-0.7, 1.]]] # correlation coefficients of the modes
-sigmas = [[1.5, 1.5], [1.5, 1.5]] # standard deviations of the modes
-weights = [0.5, 0.5] # weights of each mode
+mus = [[-5., -5.], [5., 5.]] # means of the two modes
+corrcoefs = [[[1., -0.7], [-0.7, 1.]], [[1., 0.7], [0.7, 1.]]] # correlation coefficients of the two modes
+sigmas = [[1.5, 1.5], [2.1, 2.1]] # standard deviations of the two modes
+weights = [1., 3.] # relative weights of each mode
nmodes = 2
mvg = bilby.core.prior.MultivariateGaussianDist(names, nmodes=2, mus=mus,
corrcoefs=corrcoefs,
@@ -52,20 +53,34 @@ priors = dict()
priors['m'] = bilby.core.prior.MultivariateGaussian(mvg, 'm')
priors['c'] = bilby.core.prior.MultivariateGaussian(mvg, 'c')
-# And run sampler
-# result = bilby.run_sampler(
-# likelihood=likelihood, priors=priors, sampler='pymc3',
-# outdir=outdir, draws=2000, label=label)
-
result = bilby.run_sampler(
likelihood=likelihood, priors=priors, sampler='dynesty', nlive=4000,
outdir=outdir, label=label)
-# result = bilby.run_sampler(
-# likelihood=likelihood, priors=priors, sampler='nestle', nlive=4000,
-# outdir=outdir, label=label)
+fig = result.plot_corner(save=False)
+
+# plot the priors (to show that they look correct)
+axs = fig.get_axes()
+
+# plot the 1d marginal distributions
+x = np.linspace(-12, 12, 5000)
+aidx = [0, 3]
+for j in range(2): # loop over parameters
+ gp = np.zeros(len(x))
+ for i in range(nmodes): # loop over modes
+ gp += weights[i]*stats.norm.pdf(x, loc=mus[i][j], scale=mvg.sigmas[i][j])
+ gp = gp/np.trapz(gp, x) # renormalise
+
+ axs[aidx[j]].plot(x, gp, 'k--', lw=2)
+
+# plot the 2d distribution
+for i in range(nmodes):
+ v, w = linalg.eigh(mvg.covs[i])
+ v = 2. * np.sqrt(2.) * np.sqrt(v)
+ u = w[0] / linalg.norm(w[0])
+ angle = np.arctan(u[1] / u[0])
+ angle = 180. * angle / np.pi # convert to degrees
+ ell = mpl.patches.Ellipse(mus[i], v[0], v[1], 180. + angle, edgecolor='black', facecolor='none', lw=2, ls='--')
+ axs[2].add_artist(ell)
-# result = bilby.run_sampler(
-# likelihood=likelihood, priors=priors, sampler='emcee', nsteps=1000,
-# nwalkers=200, nburn=500, outdir=outdir, label=label)
-result.plot_corner()
\ No newline at end of file
+fig.savefig('{}/{}_corner.png'.format(outdir, label), dpi=300)
diff --git a/sampler_requirements.txt b/sampler_requirements.txt
index e81411d3..3204ab30 100644
--- a/sampler_requirements.txt
+++ b/sampler_requirements.txt
@@ -3,5 +3,5 @@ dynesty
emcee
nestle
ptemcee
-pymc3
+pymc3>=3.6
pymultinest
\ No newline at end of file
--
GitLab