Beginning to merge lik_functions and derivatives with richardos

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):

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)