Added a debug examples

GPy/examples/laplace_approximations.py

@@ -35,12 +35,86 @@ def timing():
     print the_is
     print np.mean(the_is)
 
+def debug_student_t_noise_approx():
+    real_var = 0.2
+    #Start a function, any function
+    X = np.linspace(0.0, 10.0, 30)[:, None]
+    Y = np.sin(X) + np.random.randn(*X.shape)*real_var
+
+    X_full = np.linspace(0.0, 10.0, 500)[:, None]
+    Y_full = np.sin(X_full)
+
+    #Y = Y/Y.max()
+
+    #Add student t random noise to datapoints
+    deg_free = 10000
+    real_sd = np.sqrt(real_var)
+    print "Real noise: ", real_sd
+
+    initial_var_guess = 0.01
+    #t_rv = t(deg_free, loc=0, scale=real_var)
+    #noise = t_rvrvs(size=Y.shape)
+    #Y += noise
+
+    plt.figure(1)
+    plt.suptitle('Gaussian likelihood')
+    # Kernel object
+    kernel1 = GPy.kern.rbf(X.shape[1])
+    kernel2 = kernel1.copy()
+    kernel3 = kernel1.copy()
+    kernel4 = kernel1.copy()
+    kernel5 = kernel1.copy()
+    kernel6 = kernel1.copy()
+
+    print "Clean Gaussian"
+    #A GP should completely break down due to the points as they get a lot of weight
+    # create simple GP model
+    m = GPy.models.GP_regression(X, Y, kernel=kernel1)
+    # optimize
+    m.ensure_default_constraints()
+    m.optimize()
+    # plot
+    plt.subplot(131)
+    m.plot()
+    plt.plot(X_full, Y_full)
+    print m
+
+    plt.suptitle('Student-t likelihood')
+    edited_real_sd = initial_var_guess #real_sd
+
+    print "Clean student t, rasm"
+    t_distribution = GPy.likelihoods.likelihood_functions.student_t(deg_free, sigma=edited_real_sd)
+    stu_t_likelihood = GPy.likelihoods.Laplace(Y.copy(), t_distribution, rasm=True)
+    m = GPy.models.GP(X, stu_t_likelihood, kernel6)
+    m.ensure_default_constraints()
+    m.update_likelihood_approximation()
+    m.optimize()
+    print(m)
+    plt.subplot(132)
+    m.plot()
+    plt.plot(X_full, Y_full)
+    plt.ylim(-2.5, 2.5)
+
+    print "Clean student t, ncg"
+    t_distribution = GPy.likelihoods.likelihood_functions.student_t(deg_free, sigma=edited_real_sd)
+    stu_t_likelihood = GPy.likelihoods.Laplace(Y, t_distribution, rasm=False)
+    m = GPy.models.GP(X, stu_t_likelihood, kernel3)
+    m.ensure_default_constraints()
+    m.update_likelihood_approximation()
+    m.optimize()
+    print(m)
+    plt.subplot(133)
+    m.plot()
+    plt.plot(X_full, Y_full)
+    plt.ylim(-2.5, 2.5)
+
+    plt.show()
 
 def student_t_approx():
     """
     Example of regressing with a student t likelihood
     """
-    real_var = 0.1
+    real_var = 0.2
     #Start a function, any function
     X = np.linspace(0.0, 10.0, 30)[:, None]
     Y = np.sin(X) + np.random.randn(*X.shape)*real_var
@@ -58,8 +132,11 @@ def student_t_approx():
     #Yc = Yc/Yc.max()
 
     #Add student t random noise to datapoints
-    deg_free = 10
+    deg_free = 1000000000000
     real_sd = np.sqrt(real_var)
+    print "Real noise: ", real_sd
+
+    initial_var_guess = 0.01
     #t_rv = t(deg_free, loc=0, scale=real_var)
     #noise = t_rvrvs(size=Y.shape)
     #Y += noise
@@ -73,6 +150,7 @@ def student_t_approx():
     #print corrupted_indices
     #noise = t_rv.rvs(size=(len(corrupted_indices), 1))
     #Y[corrupted_indices] += noise
+
     plt.figure(1)
     plt.suptitle('Gaussian likelihood')
     # Kernel object
@@ -108,7 +186,7 @@ def student_t_approx():
 
     plt.figure(2)
     plt.suptitle('Student-t likelihood')
-    edited_real_sd = real_sd
+    edited_real_sd = initial_var_guess #real_sd
 
     print "Clean student t, rasm"
     t_distribution = GPy.likelihoods.likelihood_functions.student_t(deg_free, sigma=edited_real_sd)

GPy/likelihoods/Laplace.py

@@ -5,7 +5,7 @@ from scipy.linalg import inv, cho_solve, det
 from numpy.linalg import cond
 from GPy.likelihoods.likelihood import likelihood
 from GPy.util.linalg import pdinv, mdot, jitchol, chol_inv, det_ln_diag, pddet
-from scipy.linalg.lapack import dtrtrs
+from scipy.linalg.flapack import dtrtrs
 import pylab as plt
 
 
@@ -63,6 +63,7 @@ class Laplace(likelihood):
         return self.likelihood_function._get_param_names()
 
     def _set_params(self, p):
+        print "Setting noise sd: ", p
         return self.likelihood_function._set_params(p)
 
     def both_gradients(self, dL_d_K_Sigma, dK_dthetaK):
@@ -79,7 +80,9 @@ class Laplace(likelihood):
 
     def _shared_gradients_components(self):
         dL_dytil = -np.dot(self.Y.T, (self.K+self.Sigma_tilde))
-        dytil_dfhat = self.Wi__Ki_W # np.dot(self.Sigma_tilde, self.Ki) + np.eye(self.N) # or self.Wi__Ki_W?
+        #dytil_dfhat = self.Wi__Ki_W # np.dot(self.Sigma_tilde, self.Ki) + np.eye(self.N) # or self.Wi__Ki_W?
+        Ki = inv(self.K)
+        dytil_dfhat = np.dot(self.Sigma_tilde, Ki) + np.eye(self.N) # or self.Wi__Ki_W?
         return dL_dytil, dytil_dfhat
 
     def _Kgradients(self, dL_d_K_Sigma, dK_dthetaK):
@@ -93,19 +96,26 @@ class Laplace(likelihood):
         dL_dytil, dytil_dfhat = self._shared_gradients_components()
 
         print "Computing K gradients"
+        print "dytil_dfhat: ", np.mean(dytil_dfhat)
         I = np.eye(self.N)
         C = np.dot(self.K, self.W)
         A = I + C
         #plt.imshow(A)
         #plt.show()
-        ki, _, _, _ = pdinv(self.K)
+
-        I_KW_i, _, _, _ = pdinv(A)
+        #FIXME: K ISNT SYMMETRIC SO NEITHER IS A AND IT MAKES IT NON-PD!
+        #ki, _, _, _ = pdinv(self.K)
+        #I_KW_i, _, _, _ = pdinv(A)
+
+        I_KW_i = inv(A)
+
 
         #FIXME: Careful dK_dthetaK is not the derivative with respect to the marginal just prior K!
         #Derivative for each f dimension, for each of K's hyper parameters
         dfhat_dthetaK = np.zeros((self.f_hat.shape[0], dK_dthetaK.shape[0]))
+        grad = self.likelihood_function.link_grad(self.data, self.f_hat, self.extra_data)
         for ind_j, thetaj in enumerate(dK_dthetaK):
-            dfhat_dthetaK[:, ind_j] = mdot(I_KW_i, thetaj, self.likelihood_function.link_grad(self.data, self.f_hat, self.extra_data))
+            dfhat_dthetaK[:, ind_j] = np.dot(I_KW_i, np.dot(thetaj, grad))
 
         dytil_dthetaK = np.dot(dytil_dfhat, dfhat_dthetaK) # should be (D,thetaK)
         #FIXME: Careful dL_dK = dL_d_K_Sigma
@@ -116,8 +126,11 @@ class Laplace(likelihood):
         dSigmai_dthetaK = 0 #+ np.sum(d3phi_d3fhat*dfhat_dthetaK) #FIXME: CAREFUL OF THIS SUM! SHOULD SUM OVER FHAT NOT THETAS
         dSigma_dthetaK = -mdot(self.Sigma_tilde, dSigmai_dthetaK, self.Sigma_tilde)
 
+        print "dL_dytil: ", np.mean(dL_dytil)
+        print "dytil_dthetaK: ", np.mean(dytil_dthetaK)
         dL_dthetaK_implicit = np.sum(np.dot(dL_dytil, dytil_dthetaK), axis=0)# + np.dot(dL_dSigma, dSigma_dthetaK)
         #dL_dthetaK_implicit = np.dot(dL_dytil.T, dytil_dthetaK.T)
+        import ipdb; ipdb.set_trace()  # XXX BREAKPOINT
         return np.squeeze(dL_dthetaK_implicit)
 
     def _gradients(self, partial):

GPy/models/GP.py

@@ -116,7 +116,6 @@ class GP(model):
         """
         return -0.5 * self.D * self.K_logdet + self._model_fit_term() + self.likelihood.Z
 
-
     def _log_likelihood_gradients(self):
         """
         The gradient of all parameters.
@@ -132,9 +131,14 @@ class GP(model):
 
             dL_dthetaK_implicit = self.likelihood._Kgradients(self.dL_dK, dK_dthetaK)
             dL_dthetaK = dL_dthetaK_explicit + dL_dthetaK_implicit
+
+            print "dL_dthetaK_explicit: {dldkx}     dL_dthetaK_implicit: {dldki}        dL_dthetaK: {dldk}".format(dldkx=dL_dthetaK_explicit, dldki=dL_dthetaK_implicit, dldk=dL_dthetaK)
+
             dL_dthetaL = self.likelihood._gradients(partial=self.dL_dK)
         else:
+            print "dL_dthetaK: ", dL_dthetaK
             dL_dthetaL = self.likelihood._gradients(partial=np.diag(self.dL_dK))
+        print "dL_dthetaL: ", dL_dthetaL
         return np.hstack((dL_dthetaK, dL_dthetaL))
         #return np.hstack((self.kern.dK_dtheta(dL_dK=self.dL_dK, X=self.X), self.likelihood._gradients(partial=np.diag(self.dL_dK))))