From a52c20f47008233495e20d96b4ab50be8eb7d4a3 Mon Sep 17 00:00:00 2001 From: Alan Saul Date: Tue, 7 May 2013 13:35:47 +0100 Subject: [PATCH] Added a debug examples --- GPy/examples/laplace_approximations.py | 84 +++++++++++++++++++++++++- GPy/likelihoods/Laplace.py | 23 +++++-- GPy/models/GP.py | 6 +- 3 files changed, 104 insertions(+), 9 deletions(-) diff --git a/GPy/examples/laplace_approximations.py b/GPy/examples/laplace_approximations.py index 5d1c1224..7e5c55bf 100644 --- a/GPy/examples/laplace_approximations.py +++ b/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) diff --git a/GPy/likelihoods/Laplace.py b/GPy/likelihoods/Laplace.py index 5e28212e..02f2c93f 100644 --- a/GPy/likelihoods/Laplace.py +++ b/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): diff --git a/GPy/models/GP.py b/GPy/models/GP.py index 96ec6582..07c7a708 100644 --- a/GPy/models/GP.py +++ b/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))))