From 7877dd303febd99eb761660eaf029f12f1a591d4 Mon Sep 17 00:00:00 2001
From: Moritz Huebner <moritz.huebner@ligo.org>
Date: Thu, 23 Apr 2020 18:05:09 -0500
Subject: [PATCH] Resolve "Add multivariate gaussian to core likelihood"
---
bilby/core/likelihood.py | 64 ++++++++++++++++++++
examples/core_examples/15d_gaussian.py | 67 +--------------------
test/likelihood_test.py | 81 +++++++++++++++++++++++++-
3 files changed, 147 insertions(+), 65 deletions(-)
diff --git a/bilby/core/likelihood.py b/bilby/core/likelihood.py
index e2b391420..11fbd34cd 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 8f588e3c4..036c66fdb 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 cf878ad85..3bc3d0b97 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):
--
GitLab