From ff82f12c3d321bfc3ce6615447fad25aea9de6bd Mon Sep 17 00:00:00 2001 From: monabf <26089777+monabf@users.noreply.github.com> Date: Sat, 15 May 2021 00:21:47 +0200 Subject: [PATCH] Corrected Multivariate Gaussian prior (#775) * Corrected MultivariateGaussian prior Some corrections including adapting to current version of pdinv, correcting the expressions of constant, pdf and its gradient, and adding the printing function. After some tests, seems to run as expected, similarly to the Gaussian prior which was already working. * Added test of MultivariateGaussian prior Simple unit test for creating a kernel with Multivariate Gaussian prior over the lengthscales, then performing GP regression. * Took care of case where x of shape (n, 1) for multivariate Gaussian prior * Got rid of unnecessary asserts in Multivariate Gaussian prior since loss of time Co-authored-by: and --- GPy/core/parameterization/priors.py | 42 +++++++++++++++++++---------- GPy/testing/model_tests.py | 28 +++++++++++++++++++ 2 files changed, 56 insertions(+), 14 deletions(-) diff --git a/GPy/core/parameterization/priors.py b/GPy/core/parameterization/priors.py index 906ad65f..6269436b 100644 --- a/GPy/core/parameterization/priors.py +++ b/GPy/core/parameterization/priors.py @@ -200,37 +200,50 @@ class MultivariateGaussian(Prior): def __new__(cls, mu=0, var=1): # Singleton: if cls._instances: - cls._instances[:] = [instance for instance in cls._instances if instance()] + cls._instances[:] = [instance for instance in cls._instances if + instance()] for instance in cls._instances: - if np.all(instance().mu == mu) and np.all(instance().var == var): + if np.all(instance().mu == mu) and np.all( + instance().var == var): return instance() - o = super(Prior, cls).__new__(cls, mu, var) + newfunc = super(Prior, cls).__new__ + if newfunc is object.__new__: + o = newfunc(cls) + else: + o = newfunc(cls, mu, var) cls._instances.append(weakref.ref(o)) return cls._instances[-1]() def __init__(self, mu, var): self.mu = np.array(mu).flatten() self.var = np.array(var) - assert len(self.var.shape) == 2 - assert self.var.shape[0] == self.var.shape[1] + assert len(self.var.shape) == 2, 'Covariance must be a matrix' + assert self.var.shape[0] == self.var.shape[1], \ + 'Covariance must be a square matrix' assert self.var.shape[0] == self.mu.size self.input_dim = self.mu.size - self.inv, self.hld = pdinv(self.var) - self.constant = -0.5 * self.input_dim * np.log(2 * np.pi) - self.hld + self.inv, _, self.hld, _ = pdinv(self.var) + self.constant = -0.5 * (self.input_dim * np.log(2 * np.pi) + self.hld) + + def __str__(self): + return 'MultiN(' + str(self.mu) + ', ' + str(np.diag(self.var)) + ')' def summary(self): raise NotImplementedError def pdf(self, x): + x = np.array(x).flatten() return np.exp(self.lnpdf(x)) def lnpdf(self, x): + x = np.array(x).flatten() d = x - self.mu - return self.constant - 0.5 * np.sum(d * np.dot(d, self.inv), 1) + return self.constant - 0.5 * np.dot(d.T, np.dot(self.inv, d)) def lnpdf_grad(self, x): + x = np.array(x).flatten() d = x - self.mu - return -np.dot(self.inv, d) + return - np.dot(self.inv, d) def rvs(self, n): return np.random.multivariate_normal(self.mu, self.var, n) @@ -247,14 +260,15 @@ class MultivariateGaussian(Prior): return self.mu, self.var def __setstate__(self, state): - self.mu = state[0] + self.mu = np.array(state[0]).flatten() self.var = state[1] - assert len(self.var.shape) == 2 - assert self.var.shape[0] == self.var.shape[1] + assert len(self.var.shape) == 2, 'Covariance must be a matrix' + assert self.var.shape[0] == self.var.shape[1], \ + 'Covariance must be a square matrix' assert self.var.shape[0] == self.mu.size self.input_dim = self.mu.size - self.inv, self.hld = pdinv(self.var) - self.constant = -0.5 * self.input_dim * np.log(2 * np.pi) - self.hld + self.inv, _, self.hld, _ = pdinv(self.var) + self.constant = -0.5 * (self.input_dim * np.log(2 * np.pi) + self.hld) def gamma_from_EV(E, V): warnings.warn("use Gamma.from_EV to create Gamma Prior", FutureWarning) diff --git a/GPy/testing/model_tests.py b/GPy/testing/model_tests.py index 9ee16f99..a241fc9c 100644 --- a/GPy/testing/model_tests.py +++ b/GPy/testing/model_tests.py @@ -936,6 +936,34 @@ class GradientTests(np.testing.TestCase): #import ipdb;ipdb.set_trace() #m.constrain_fixed('.*rbf_var', 1.) self.assertTrue(m.checkgrad()) + + def test_simple_MultivariateGaussian_prior(self): + X = np.random.multivariate_normal( + [1, 5], np.diag([0.5, 0.3]), (100, 1)).reshape(100, 2) + Y = X + np.random.randn(100, 2) * 0.05 + kernel = GPy.kern.RBF(input_dim=2, variance=1,lengthscale=1, ARD=True) + kernel.unconstrain() + kernel.variance.set_prior(GPy.priors.Gaussian(150, 5)) + kernel.lengthscale.set_prior(GPy.priors.MultivariateGaussian( + np.array([20, 20]), np.diag([5, 5]))) + m = GPy.models.GPRegression(X, Y, kernel=kernel) + m.optimize() + print(m.kern.variance) + print(m.kern.lengthscale) + + def test_simple_MultivariateGaussian_prior_matrixmean(self): + X = np.random.multivariate_normal( + [1, 5], np.diag([0.5, 0.3]), (100, 1)).reshape(100, 2) + Y = X + np.random.randn(100, 2) * 0.05 + kernel = GPy.kern.RBF(input_dim=2, variance=1,lengthscale=1, ARD=True) + kernel.unconstrain() + kernel.variance.set_prior(GPy.priors.Gaussian(150, 5)) + kernel.lengthscale.set_prior(GPy.priors.MultivariateGaussian( + np.array([[20, 20]]), np.diag([5, 5]))) + m = GPy.models.GPRegression(X, Y, kernel=kernel) + m.optimize() + print(m.kern.variance) + print(m.kern.lengthscale) def test_multioutput_sparse_regression_1D(self): X1 = np.random.rand(500, 1) * 8