diff --git a/bilby/core/likelihood.py b/bilby/core/likelihood.py index e2b39142034cb6712fa12f919d49a4d589abf1c5..11fbd34cd649f0286cb3582ec31dd721f9bb83d2 100644 --- a/bilby/core/likelihood.py +++ b/bilby/core/likelihood.py @@ -3,6 +3,7 @@ import copy import numpy as np from scipy.special import gammaln +from scipy.stats import multivariate_normal from .utils import infer_parameters_from_function @@ -401,6 +402,69 @@ class StudentTLikelihood(Analytical1DLikelihood): self._nu = nu +class AnalyticalMultidimensionalCovariantGaussian(Likelihood): + """ + A multivariate Gaussian likelihood + with known analytic solution. + + Parameters + ---------- + mean: array_like + Array with the mean values of distribution + cov: array_like + The ndim*ndim covariance matrix + """ + + def __init__(self, mean, cov): + self.cov = np.atleast_2d(cov) + self.mean = np.atleast_1d(mean) + self.sigma = np.sqrt(np.diag(self.cov)) + self.pdf = multivariate_normal(mean=self.mean, cov=self.cov) + parameters = {"x{0}".format(i): 0 for i in range(self.dim)} + super(AnalyticalMultidimensionalCovariantGaussian, self).__init__(parameters=parameters) + + @property + def dim(self): + return len(self.cov[0]) + + def log_likelihood(self): + x = np.array([self.parameters["x{0}".format(i)] for i in range(self.dim)]) + return self.pdf.logpdf(x) + + +class AnalyticalMultidimensionalBimodalCovariantGaussian(Likelihood): + """ + A multivariate Gaussian likelihood + with known analytic solution. + + Parameters + ---------- + mean_1: array_like + Array with the mean value of the first mode + mean_2: array_like + Array with the mean value of the second mode + cov: array_like + """ + + def __init__(self, mean_1, mean_2, cov): + self.cov = np.atleast_2d(cov) + self.sigma = np.sqrt(np.diag(self.cov)) + self.mean_1 = np.atleast_1d(mean_1) + self.mean_2 = np.atleast_1d(mean_2) + self.pdf_1 = multivariate_normal(mean=self.mean_1, cov=self.cov) + self.pdf_2 = multivariate_normal(mean=self.mean_2, cov=self.cov) + parameters = {"x{0}".format(i): 0 for i in range(self.dim)} + super(AnalyticalMultidimensionalBimodalCovariantGaussian, self).__init__(parameters=parameters) + + @property + def dim(self): + return len(self.cov[0]) + + def log_likelihood(self): + x = np.array([self.parameters["x{0}".format(i)] for i in range(self.dim)]) + return -np.log(2) + np.logaddexp(self.pdf_1.logpdf(x), self.pdf_2.logpdf(x)) + + class JointLikelihood(Likelihood): def __init__(self, *likelihoods): """ diff --git a/examples/core_examples/15d_gaussian.py b/examples/core_examples/15d_gaussian.py index 8f588e3c4a7835465a5d5b48238154b346fcb349..036c66fdbd152b4bfd04514abb4fc463bc7a6730 100644 --- a/examples/core_examples/15d_gaussian.py +++ b/examples/core_examples/15d_gaussian.py @@ -1,8 +1,9 @@ -from scipy.stats import multivariate_normal - import bilby import numpy as np +from bilby.core.likelihood import AnalyticalMultidimensionalCovariantGaussian, \ + AnalyticalMultidimensionalBimodalCovariantGaussian + logger = bilby.core.utils.logger cov = [ @@ -58,68 +59,6 @@ cov = [ dim = 15 mean = np.zeros(dim) - -class AnalyticalMultidimensionalCovariantGaussian(bilby.Likelihood): - """ - A multivariate Gaussian likelihood - with known analytic solution. - - Parameters - ---------- - mean: array_like - Array with the mean value of distribution - cov: array_like - The ndim*ndim covariance matrix - """ - - def __init__(self, mean, cov): - super(AnalyticalMultidimensionalCovariantGaussian, self).__init__(parameters=dict()) - self.cov = np.array(cov) - self.mean = np.array(mean) - self.sigma = np.sqrt(np.diag(self.cov)) - self.pdf = multivariate_normal(mean=self.mean, cov=self.cov) - - @property - def dim(self): - return len(self.cov[0]) - - def log_likelihood(self): - x = np.array([self.parameters["x{0}".format(i)] for i in range(self.dim)]) - return self.pdf.logpdf(x) - - -class AnalyticalMultidimensionalBimodalCovariantGaussian(bilby.Likelihood): - """ - A multivariate Gaussian likelihood - with known analytic solution. - - Parameters - ---------- - mean_1: array_like - Array with the mean value of the first mode - mean_2: array_like - Array with the mean value of the second mode - cov: array_like - """ - - def __init__(self, mean_1, mean_2, cov): - super(AnalyticalMultidimensionalBimodalCovariantGaussian, self).__init__(parameters=dict()) - self.cov = np.array(cov) - self.mean_1 = np.array(mean_1) - self.mean_2 = np.array(mean_2) - self.sigma = np.sqrt(np.diag(self.cov)) - self.pdf_1 = multivariate_normal(mean=self.mean_1, cov=self.cov) - self.pdf_2 = multivariate_normal(mean=self.mean_2, cov=self.cov) - - @property - def dim(self): - return len(self.cov[0]) - - def log_likelihood(self): - x = np.array([self.parameters["x{0}".format(i)] for i in range(self.dim)]) - return -np.log(2) + np.logaddexp(self.pdf_1.logpdf(x), self.pdf_2.logpdf(x)) - - label = "multidim_gaussian_unimodal" outdir = "outdir" diff --git a/test/likelihood_test.py b/test/likelihood_test.py index cf878ad85d4407bcceabd6a0d641cace791a34b6..3bc3d0b97b541af86b7d3dd10bc8e2941a82801d 100644 --- a/test/likelihood_test.py +++ b/test/likelihood_test.py @@ -6,7 +6,9 @@ import mock import numpy as np from bilby.core.likelihood import ( Likelihood, GaussianLikelihood, PoissonLikelihood, StudentTLikelihood, - Analytical1DLikelihood, ExponentialLikelihood, JointLikelihood) + Analytical1DLikelihood, ExponentialLikelihood, + AnalyticalMultidimensionalCovariantGaussian, + AnalyticalMultidimensionalBimodalCovariantGaussian, JointLikelihood) class TestLikelihoodBase(unittest.TestCase): @@ -490,6 +492,83 @@ class TestExponentialLikelihood(unittest.TestCase): self.assertEqual(expected, repr(self.exponential_likelihood)) +class TestAnalyticalMultidimensionalCovariantGaussian(unittest.TestCase): + + def setUp(self): + self.cov = [[1, 0, 0], [0, 4, 0], [0, 0, 9]] + self.sigma = [1, 2, 3] + self.mean = [10, 11, 12] + self.likelihood = AnalyticalMultidimensionalCovariantGaussian( + mean=self.mean, + cov=self.cov) + + def tearDown(self): + del self.cov + del self.sigma + del self.mean + del self.likelihood + + def test_cov(self): + self.assertTrue(np.array_equal(self.cov, self.likelihood.cov)) + + def test_mean(self): + self.assertTrue(np.array_equal(self.mean, self.likelihood.mean)) + + def test_sigma(self): + self.assertTrue(np.array_equal(self.sigma, self.likelihood.sigma)) + + def test_parameters(self): + self.assertDictEqual(dict(x0=0, x1=0, x2=0), self.likelihood.parameters) + + def test_dim(self): + self.assertEqual(3, self.likelihood.dim) + + def test_log_likelihood(self): + likelihood = AnalyticalMultidimensionalCovariantGaussian(mean=[0], cov=[1]) + self.assertEqual(-np.log(2*np.pi)/2, likelihood.log_likelihood()) + + +class TestAnalyticalMultidimensionalBimodalCovariantGaussian(unittest.TestCase): + + def setUp(self): + self.cov = [[1, 0, 0], [0, 4, 0], [0, 0, 9]] + self.sigma = [1, 2, 3] + self.mean_1 = [10, 11, 12] + self.mean_2 = [20, 21, 22] + self.likelihood = AnalyticalMultidimensionalBimodalCovariantGaussian( + mean_1=self.mean_1, + mean_2=self.mean_2, + cov=self.cov) + + def tearDown(self): + del self.cov + del self.sigma + del self.mean_1 + del self.mean_2 + del self.likelihood + + def test_cov(self): + self.assertTrue(np.array_equal(self.cov, self.likelihood.cov)) + + def test_mean_1(self): + self.assertTrue(np.array_equal(self.mean_1, self.likelihood.mean_1)) + + def test_mean_2(self): + self.assertTrue(np.array_equal(self.mean_2, self.likelihood.mean_2)) + + def test_sigma(self): + self.assertTrue(np.array_equal(self.sigma, self.likelihood.sigma)) + + def test_parameters(self): + self.assertDictEqual(dict(x0=0, x1=0, x2=0), self.likelihood.parameters) + + def test_dim(self): + self.assertEqual(3, self.likelihood.dim) + + def test_log_likelihood(self): + likelihood = AnalyticalMultidimensionalBimodalCovariantGaussian(mean_1=[0], mean_2=[0], cov=[1]) + self.assertEqual(-np.log(2*np.pi)/2, likelihood.log_likelihood()) + class TestJointLikelihood(unittest.TestCase): def setUp(self):