Commit 7877dd30 authored by Moritz Huebner's avatar Moritz Huebner Committed by Gregory Ashton
Browse files

Resolve "Add multivariate gaussian to core likelihood"

parent ce2a4047
......@@ -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):
"""
......
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"
......
......@@ -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):
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment