Student t likelihood function checkgrads (summed gradients wrt to sigma2), maybe some numerical instability in laplace

GPy/likelihoods/laplace.py

@@ -127,7 +127,6 @@ class Laplace(likelihood):
                              #- 0.5*np.trace(mdot(self.Ki_W_i, (self.K, np.diagflat(dlik_hess_dthetaL[thetaL_i]))))
                              + np.dot(0.5*np.diag(self.Ki_W_i)[:,None].T, dlik_hess_dthetaL[thetaL_i])
                              )
-            import ipdb; ipdb.set_trace()  # XXX BREAKPOINT
 
             #Implicit
             df_hat_dthetaL = mdot(I_KW_i, self.K, dlik_grad_dthetaL[thetaL_i])
@@ -203,7 +202,7 @@ class Laplace(likelihood):
         #self.cC = 0.5*self.y_Wi_Ki_i_y
         #self.dD = -0.5*self.ln_B_det
         #print "Ztilde: {} lik: {} a: {} b: {} c: {} d:".format(Z_tilde, self.lik, self.aA, self.bB, self.cC, self.dD)
-        print "param value: {}".format(self.likelihood_function._get_params())
+        #print "param value: {}".format(self.likelihood_function._get_params())
 
         #Convert to float as its (1, 1) and Z must be a scalar
         self.Z = np.float64(Z_tilde)
@@ -330,7 +329,6 @@ class Laplace(likelihood):
         self.old_before_s = self.likelihood_function._get_params()
         #print "before: ", self.old_before_s
         #if self.old_before_s < 1e-5:
-            #import ipdb; ipdb.set_trace() ### XXX BREAKPOINT
 
         #old_a = np.zeros((self.N, 1))
         if self.old_a is None:
@@ -384,7 +382,6 @@ class Laplace(likelihood):
             new_obj = sp.optimize.minimize_scalar(i_o, method='brent', tol=1e-4, options={'maxiter':20, 'disp':True}).fun
             f = self.f.copy()
             a = self.a.copy()
-            #import ipdb; ipdb.set_trace() ### XXX BREAKPOINT
 
             #f_old = f.copy()
             #update_passed = False

GPy/likelihoods/likelihood_functions.py

@@ -218,16 +218,11 @@ class StudentT(LikelihoodFunction):
         """
         assert y.shape == f.shape
         e = y - f
-        #A = gammaln((self.v + 1) * 0.5)
-        #B = - gammaln(self.v * 0.5)
-        #C = - 0.5*np.log(self.sigma2 * self.v * np.pi)
-        #D = + (-(self.v + 1)*0.5)*np.log(1 + ((e**2)/self.sigma2)/np.float(self.v))
         objective = (+ gammaln((self.v + 1) * 0.5)
                      - gammaln(self.v * 0.5)
                      - 0.5*np.log(self.sigma2 * self.v * np.pi)
-                     + (-(self.v + 1)*0.5)*np.log(1 + ((e**2)/self.sigma2)/np.float(self.v))
+                     - 0.5*(self.v + 1)*np.log(1 + (1/np.float(self.v))*((e**2)/self.sigma2))
                     )
-        #print "C: {} D: {} obj: {}".format(C, np.sum(D), objective.sum())
         return np.sum(objective)
 
     def dlik_df(self, y, f, extra_data=None):
@@ -291,9 +286,13 @@ class StudentT(LikelihoodFunction):
         """
         assert y.shape == f.shape
         e = y - f
-        #FIXME: OUT BY SOME FUNCTION OF N
+        #FIXME: OUT BY SOME FUNCTION OF N, or the fact that we are summing over several things in the objective?
         dlik_dvar = self.v*(e**2 - self.sigma2)/(2*self.sigma2*(self.sigma2*self.v + e**2))
-        return dlik_dvar
+        #dlik_dvar = ( 0.5*(1/float(self.sigma2))
+                     #-0.5*(self.v + 1)*(-(1/float(self.v))*(e**2)/(1/(float(self.sigma2**2))))
+                     #/ (1 + (1/float(self.v))*((e**2)/float(self.sigma2)))
+                     #)
+        return np.sum(dlik_dvar) #May not want to sum over all dimensions if using many D?
 
     def dlik_df_dvar(self, y, f, extra_data=None):
         """
@@ -516,8 +515,7 @@ class Gaussian(LikelihoodFunction):
         e = y - f
         s_4 = 1.0/(self._variance**2)
         dlik_dsigma = -0.5*self.N/self._variance + 0.5*s_4*np.dot(e.T, e)
-        #dlik_dsigma = -0.5*self.N + 0.5*s_4*np.dot(e.T, e)
+        return np.sum(dlik_dsigma) # Sure about this sum?
-        return dlik_dsigma
 
     def dlik_df_dvar(self, y, f, extra_data=None):
         """

GPy/testing/laplace_tests.py

@@ -45,6 +45,7 @@ def dparam_checkgrad(func, dfunc, params, args, constrain_positive=True, randomi
         for fixed_val in range(dfnum):
             #dlik and dlik_dvar gives back 1 value for each
             f_ind = min(fnum, fixed_val+1) - 1
+            print "fnum: {} dfnum: {} f_ind: {} fixed_val: {}".format(fnum, dfnum, f_ind, fixed_val)
             grad = GradientChecker(lambda x: np.atleast_1d(partial_f(x))[f_ind],
                                    lambda x : np.atleast_1d(partial_df(x))[fixed_val],
                                    param, 'p')
@@ -63,9 +64,9 @@ def dparam_checkgrad(func, dfunc, params, args, constrain_positive=True, randomi
 
 class LaplaceTests(unittest.TestCase):
     def setUp(self):
-        self.N = 1
+        self.N = 5
-        self.D = 5
+        self.D = 1
-        self.X = np.linspace(0, 1, self.N)[:, None]
+        self.X = np.linspace(0, self.D, self.N)[:, None]
 
         self.real_std = 0.2
         noise = np.random.randn(*self.X.shape)*self.real_std
@@ -104,6 +105,27 @@ class LaplaceTests(unittest.TestCase):
         grad.checkgrad(verbose=1)
         self.assertTrue(grad.checkgrad())
 
+    def test_gaussian_d2lik_d2f_2(self):
+        print "\n{}".format(inspect.stack()[0][3])
+        self.Y = None
+        self.gauss = None
+
+        self.N = 2
+        self.D = 1
+        self.X = np.linspace(0, self.D, self.N)[:, None]
+        self.real_std = 0.2
+        noise = np.random.randn(*self.X.shape)*self.real_std
+        self.Y = np.sin(self.X*2*np.pi) + noise
+        self.f = np.random.rand(self.N, 1)
+        self.gauss = GPy.likelihoods.functions.Gaussian(self.var, self.D, self.N)
+
+        dlik_df = functools.partial(self.gauss.dlik_df, self.Y)
+        d2lik_d2f = functools.partial(self.gauss.d2lik_d2f, self.Y)
+        grad = GradientChecker(dlik_df, d2lik_d2f, self.f.copy(), 'f')
+        grad.randomize()
+        grad.checkgrad(verbose=1)
+        self.assertTrue(grad.checkgrad())
+
     def test_gaussian_d3lik_d3f(self):
         print "\n{}".format(inspect.stack()[0][3])
         d2lik_d2f = functools.partial(self.gauss.d2lik_d2f, self.Y)