Starting to fiddle with mode finding code

GPy/examples/laplace_approximations.py

@@ -36,7 +36,7 @@ def timing():
     print np.mean(the_is)
 
 def v_fail_test():
-    plt.close('all')
+    #plt.close('all')
     real_var = 0.1
     X = np.linspace(0.0, 10.0, 50)[:, None]
     Y = np.sin(X) + np.random.randn(*X.shape)*real_var
@@ -57,6 +57,7 @@ def v_fail_test():
     stu_t_likelihood = GPy.likelihoods.Laplace(Y.copy(), t_distribution, rasm=True)
     m = GPy.models.GP(X, stu_t_likelihood, kernel1)
     m.constrain_fixed('white', 1)
+    m.constrain_positive('t_noise')
     vs = 15
     noises = 40
     checkgrads = np.zeros((vs, noises))
@@ -64,23 +65,24 @@ def v_fail_test():
     for v_ind, v in enumerate(np.linspace(1, 20, vs)):
         m.likelihood.likelihood_function.v = v
         print v
-        for noise_ind, noise in enumerate(np.linspace(0.0000001, 1, noises)):
+        for noise_ind, noise in enumerate(np.linspace(0.0001, 1, noises)):
             m['t_noise'] = noise
             m.update_likelihood_approximation()
             checkgrads[v_ind, noise_ind] = m.checkgrad()
             vs_noises[v_ind, noise_ind] = (float(v)/(float(v) - 2))*(noise**2)
 
-    plt.figure(1)
+    plt.figure()
     plt.title('Checkgrads')
     plt.imshow(checkgrads, interpolation='nearest')
     plt.xlabel('noise')
     plt.ylabel('v')
 
-    plt.figure(2)
+    plt.figure()
     plt.title('variance change')
     plt.imshow(vs_noises, interpolation='nearest')
     plt.xlabel('noise')
     plt.ylabel('v')
+    import ipdb; ipdb.set_trace() ### XXX BREAKPOINT
     print(m)
 
 def debug_student_t_noise_approx():
@@ -139,13 +141,13 @@ def debug_student_t_noise_approx():
     stu_t_likelihood = GPy.likelihoods.Laplace(Y.copy(), t_distribution, rasm=True)
     m = GPy.models.GP(X, stu_t_likelihood, kernel6)
     #m['white'] = 1e-3
-    m.constrain_fixed('rbf_v', 1.0898)
-    m.constrain_fixed('rbf_l', 1.8651)
+    #m.constrain_fixed('rbf_v', 1.0898)
+    #m.constrain_fixed('rbf_l', 1.8651)
     #m.constrain_fixed('t_noise_variance', real_sd)
     #m.constrain_positive('rbf')
-    m.constrain_positive('t_noise')
+    #m.constrain_positive('t_noise')
+    m.constrain_positive('')
     #m.constrain_fixed('t_noi', real_sd)
-    m.ensure_default_constraints()
     m.update_likelihood_approximation()
     #m.optimize(messages=True)
     print(m)

GPy/likelihoods/Laplace.py

@@ -68,8 +68,7 @@ class Laplace(likelihood):
     def _shared_gradients_components(self):
         #FIXME: Careful of side effects! And make sure W and K are up to date!
         d3lik_d3fhat = self.likelihood_function.d3lik_d3f(self.data, self.f_hat)
-        dL_dfhat = -0.5*(np.diag(self.Ki_W_i)[:, None]*d3lik_d3fhat)
-
+        dL_dfhat = -0.5*(np.diag(self.Ki_W_i)[:, None]*d3lik_d3fhat).T
         I_KW_i = np.eye(self.N) - np.dot(self.K, self.Wi_K_i)
         return dL_dfhat, I_KW_i
 
@@ -81,10 +80,10 @@ class Laplace(likelihood):
         dlp = self.likelihood_function.dlik_df(self.data, self.f_hat)
 
         #Implicit
-        impl = mdot(dlp, dL_dfhat.T, I_KW_i)
+        impl = mdot(dlp, dL_dfhat, I_KW_i)
         expl_a = mdot(self.Ki_f, self.Ki_f.T)
         expl_b = self.Wi_K_i
-        expl = 0.5*expl_a + 0.5*expl_b
+        expl = 0.5*expl_a - 0.5*expl_b # Might need to be -?
         dL_dthetaK_exp = dK_dthetaK(expl, X)
         dL_dthetaK_imp = dK_dthetaK(impl, X)
         #print "dL_dthetaK_exp: {}     dL_dthetaK_implicit: {}".format(dL_dthetaK_exp, dL_dthetaK_imp)
@@ -103,10 +102,11 @@ class Laplace(likelihood):
         for thetaL_i in range(num_params):
             #Explicit
             dL_dthetaL_exp = np.sum(dlik_dthetaL[thetaL_i]) - 0.5*np.dot(np.diag(self.Ki_W_i), dlik_hess_dthetaL[thetaL_i])
+            #dL_dthetaL_exp = np.sum(dlik_dthetaL[thetaL_i]) - 0.5*np.trace(mdot(self.Bi, self.K, dlik_hess_dthetaL[thetaL_i]))
             #Implicit
             df_hat_dthetaL = mdot(I_KW_i, self.K, dlik_grad_dthetaL[thetaL_i])
-            dL_dthetaL_imp = np.dot(dL_dfhat.T, df_hat_dthetaL)
-            #print "dL_dthetaL_exp: {}     dL_dthetaL_implicit: {}".format(dL_dthetaL_exp, dL_dthetaL_imp)
+            dL_dthetaL_imp = np.dot(dL_dfhat, df_hat_dthetaL)
+            print "dL_dthetaL_exp: {}     dL_dthetaL_implicit: {}".format(dL_dthetaL_exp, dL_dthetaL_imp)
             dL_dthetaL[thetaL_i] = dL_dthetaL_imp + dL_dthetaL_exp
 
         return dL_dthetaL #should be array of length *params-being optimized*, for student t just optimising 1 parameter, this is (1,)

GPy/likelihoods/likelihood_functions.py

@@ -192,7 +192,6 @@ class student_t(likelihood_function):
 
         """
         assert y.shape == f.shape
-
         e = y - f
         objective = (+ gammaln((self.v + 1) * 0.5)
                      - gammaln(self.v * 0.5)