From 77bca5547055bb76ef66b9ba132661bbdc631761 Mon Sep 17 00:00:00 2001 From: Alan Saul Date: Mon, 7 Oct 2013 15:28:40 +0100 Subject: [PATCH] Beginning to merge lik_functions and derivatives with richardos --- .../noise_models/gaussian_noise.py | 29 +++++++++++--- GPy/testing/laplace_tests.py | 39 ++++++++++++++++--- 2 files changed, 57 insertions(+), 11 deletions(-) diff --git a/GPy/likelihoods/noise_models/gaussian_noise.py b/GPy/likelihoods/noise_models/gaussian_noise.py index df351cf1..afd5d297 100644 --- a/GPy/likelihoods/noise_models/gaussian_noise.py +++ b/GPy/likelihoods/noise_models/gaussian_noise.py @@ -33,7 +33,8 @@ class Gaussian(NoiseDistribution): self.I = np.eye(self.N) self.covariance_matrix = self.I * self.variance self.Ki = self.I*(1.0 / self.variance) - self.ln_det_K = np.sum(np.log(np.diag(self.covariance_matrix))) + #self.ln_det_K = np.sum(np.log(np.diag(self.covariance_matrix))) + self.ln_det_K = self.N*np.log(self.variance) def _laplace_gradients(self, y, f, extra_data=None): #must be listed in same order as 'get_param_names' @@ -81,10 +82,26 @@ class Gaussian(NoiseDistribution): def _mass(self,gp,obs): #return std_norm_pdf( (self.gp_link.transf(gp)-obs)/np.sqrt(self.variance) ) - return stats.norm.pdf(obs,self.gp_link.transf(gp),np.sqrt(self.variance)) + #Assumes no covariance, exp, sum, log for numerical stability + return np.exp(np.sum(np.log(stats.norm.pdf(obs,self.gp_link.transf(gp),np.sqrt(self.variance))))) - def _nlog_mass(self,gp,obs): - return .5*((self.gp_link.transf(gp)-obs)**2/self.variance + np.log(2.*np.pi*self.variance)) + def _nlog_mass(self,gp,obs, extra_data=None): + """ + Negative Log likelihood function + + .. math:: + \\-ln p(y_{i}|f_{i}) = +\\frac{D \\ln 2\\pi}{2} + \\frac{\\ln |K|}{2} + \\frac{(y_{i} - f_{i})^{T}\\sigma^{-2}(y_{i} - f_{i})}{2} + + :param y: data + :type y: Nx1 array + :param f: latent variables f + :type f: Nx1 array + :param extra_data: extra_data which is not used in student t distribution - not used + :returns: likelihood evaluated for this point + :rtype: float + """ + assert gp.shape == obs.shape + return .5*(np.sum((self.gp_link.transf(gp)-obs)**2/self.variance) + self.ln_det_K + self.N*np.log(2.*np.pi)) def _dnlog_mass_dgp(self,gp,obs): return (self.gp_link.transf(gp)-obs)/self.variance * self.gp_link.dtransf_df(gp) @@ -139,7 +156,7 @@ class Gaussian(NoiseDistribution): """ assert y.shape == f.shape e = y - f - objective = (- 0.5*self.D*np.log(2*np.pi) + objective = (- 0.5*self.N*np.log(2*np.pi) - 0.5*self.ln_det_K - (0.5/self.variance)*np.sum(np.square(e)) # As long as K is diagonal ) @@ -206,7 +223,7 @@ class Gaussian(NoiseDistribution): :rtype: Nx1 array """ assert y.shape == f.shape - d3lik_d3f = np.diagonal(0*self.I)[:, None] # FIXME: CAREFUL THIS MAY NOT WORK WITH MULTIDIMENSIONS? + d3lik_d3f = np.diagonal(0*self.I)[:, None] return d3lik_d3f def dlik_dvar(self, y, f, extra_data=None): diff --git a/GPy/testing/laplace_tests.py b/GPy/testing/laplace_tests.py index e1876296..acd60b4a 100644 --- a/GPy/testing/laplace_tests.py +++ b/GPy/testing/laplace_tests.py @@ -64,18 +64,16 @@ def dparam_checkgrad(func, dfunc, params, args, constrain_positive=True, randomi class LaplaceTests(unittest.TestCase): def setUp(self): - self.N = 5 + self.N = 50 self.D = 3 self.X = np.random.rand(self.N, self.D)*10 self.real_std = 0.1 noise = np.random.randn(*self.X[:, 0].shape)*self.real_std self.Y = (np.sin(self.X[:, 0]*2*np.pi) + noise)[:, None] - #self.Y = np.array([[1.0]])#np.sin(self.X*2*np.pi) + noise - self.var = 0.2 - self.f = np.random.rand(self.N, 1) - #self.f = np.array([[3.0]])#np.sin(self.X*2*np.pi) + noise + + self.var = 0.2 self.var = np.random.rand(1) self.stu_t = GPy.likelihoods.student_t(deg_free=5, sigma2=self.var) @@ -91,6 +89,37 @@ class LaplaceTests(unittest.TestCase): self.f = None self.X = None + def test_lik_mass(self): + print "\n{}".format(inspect.stack()[0][3]) + np.testing.assert_almost_equal( + np.sum(self.gauss._nlog_mass(self.f.copy(), self.Y.copy())), + -self.gauss.lik_function(self.Y.copy(), self.f.copy())) + + def test_mass_nlog_mass(self): + print "\n{}".format(inspect.stack()[0][3]) + np.testing.assert_almost_equal( + -np.log(self.gauss._mass(self.f.copy(), self.Y.copy())), + self.gauss._nlog_mass(self.f.copy(), self.Y.copy())) + + def test_gaussian_dnlog_mass_dgp(self): + print "\n{}".format(inspect.stack()[0][3]) + link = functools.partial(self.gauss._nlog_mass, obs=self.Y) + dlik_df = functools.partial(self.gauss._dnlog_mass_dgp, obs=self.Y) + grad = GradientChecker(link, dlik_df, self.f.copy(), 'g') + grad.randomize() + grad.checkgrad(verbose=1) + self.assertTrue(grad.checkgrad()) + + def test_gaussian_d2nlog_mass_d2gp(self): + print "\n{}".format(inspect.stack()[0][3]) + link = functools.partial(self.gauss._dnlog_mass_dgp, obs=self.Y) + dlik_df = functools.partial(self.gauss._d2nlog_mass_dgp2, obs=self.Y) + grad = GradientChecker(link, dlik_df, self.f.copy(), 'g') + grad.randomize() + grad.checkgrad(verbose=1) + self.assertTrue(grad.checkgrad()) + + def test_gaussian_dlik_df(self): print "\n{}".format(inspect.stack()[0][3]) link = functools.partial(self.gauss.lik_function, self.Y)