Skip to content
Snippets Groups Projects
Commit 07598631 authored by Matthew David Pitkin's avatar Matthew David Pitkin
Browse files

Update the multivariate Gaussian likelihood prior example

parent 7ecf9b8a
No related branches found
No related tags found
No related merge requests found
......@@ -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)
......@@ -3,5 +3,5 @@ dynesty
emcee
nestle
ptemcee
pymc3
pymc3>=3.6
pymultinest
\ No newline at end of file
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment