From 5e88a885b127163a83336b3773894a2f76a924e9 Mon Sep 17 00:00:00 2001 From: Alan Saul Date: Fri, 13 Sep 2013 18:01:41 +0100 Subject: [PATCH] Student t likelihood function checkgrads (summed gradients wrt to sigma2), maybe some numerical instability in laplace --- GPy/likelihoods/laplace.py | 5 +---- GPy/likelihoods/likelihood_functions.py | 18 +++++++--------- GPy/testing/laplace_tests.py | 28 ++++++++++++++++++++++--- 3 files changed, 34 insertions(+), 17 deletions(-) diff --git a/GPy/likelihoods/laplace.py b/GPy/likelihoods/laplace.py index 2897e1de..7cc4834a 100644 --- a/GPy/likelihoods/laplace.py +++ b/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 diff --git a/GPy/likelihoods/likelihood_functions.py b/GPy/likelihoods/likelihood_functions.py index 39367734..b2f9ded7 100644 --- a/GPy/likelihoods/likelihood_functions.py +++ b/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 dlik_dsigma + return np.sum(dlik_dsigma) # Sure about this sum? def dlik_df_dvar(self, y, f, extra_data=None): """ diff --git a/GPy/testing/laplace_tests.py b/GPy/testing/laplace_tests.py index 7fc6f2f4..a52cc3cd 100644 --- a/GPy/testing/laplace_tests.py +++ b/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.D = 5 - self.X = np.linspace(0, 1, self.N)[:, None] + self.N = 5 + 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 @@ -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)