From d4bfd99c21c835e5cf7873e20295561c031d5221 Mon Sep 17 00:00:00 2001 From: Alan Saul Date: Thu, 20 Jun 2013 14:30:25 +0100 Subject: [PATCH] Starting to fiddle with mode finding code --- GPy/examples/laplace_approximations.py | 18 ++++++++++-------- GPy/likelihoods/Laplace.py | 12 ++++++------ GPy/likelihoods/likelihood_functions.py | 1 - 3 files changed, 16 insertions(+), 15 deletions(-) diff --git a/GPy/examples/laplace_approximations.py b/GPy/examples/laplace_approximations.py index 887e35ae..d300806f 100644 --- a/GPy/examples/laplace_approximations.py +++ b/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) diff --git a/GPy/likelihoods/Laplace.py b/GPy/likelihoods/Laplace.py index 70ec568a..ed3229a9 100644 --- a/GPy/likelihoods/Laplace.py +++ b/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,) diff --git a/GPy/likelihoods/likelihood_functions.py b/GPy/likelihoods/likelihood_functions.py index d6dbf55f..4d298122 100644 --- a/GPy/likelihoods/likelihood_functions.py +++ b/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)