apunkt/GPy
1
0
Fork 0
mirror of https://github.com/SheffieldML/GPy.git synced 2026-06-11 15:15:15 +02:00
Code Issues Projects Releases Packages Wiki Activity Actions

Fixed a sign wrong, now gradchecks weirdly only above certain points

Browse source
This commit is contained in:
Alan Saul 2013-06-19 16:13:11 +01:00
parent de689fa8e9
commit e900509a7c
3 changed files with 64 additions and 51 deletions
Download patch file Download diff file Expand all files Collapse all files

61
GPy/examples/laplace_approximations.py
View file

@ -35,6 +35,54 @@ def timing():
print the_is
print np.mean(the_is)
def v_fail_test():
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
Y = Y/Y.max()
#Add student t random noise to datapoints
deg_free = 10
real_sd = np.sqrt(real_var)
print "Real noise std: ", real_sd
kernel1 = GPy.kern.white(X.shape[1]) #+ GPy.kern.white(X.shape[1])
edited_real_sd = 0.3#real_sd
edited_real_sd = 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, kernel1)
m.constrain_fixed('white', 1)
vs = 15
noises = 40
checkgrads = np.zeros((vs, noises))
vs_noises = np.zeros((vs, noises))
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)):
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.title('Checkgrads')
plt.imshow(checkgrads, interpolation='nearest')
plt.xlabel('noise')
plt.ylabel('v')
plt.figure(2)
plt.title('variance change')
plt.imshow(vs_noises, interpolation='nearest')
plt.xlabel('noise')
plt.ylabel('v')
print(m)
def debug_student_t_noise_approx():
plot = False
real_var = 0.1
@ -49,7 +97,7 @@ def debug_student_t_noise_approx():
Y = Y/Y.max()
#Add student t random noise to datapoints
deg_free = 1000
deg_free = 10
real_sd = np.sqrt(real_var)
print "Real noise std: ", real_sd
@ -60,7 +108,7 @@ def debug_student_t_noise_approx():
plt.close('all')
# Kernel object
kernel1 = GPy.kern.rbf(X.shape[1]) + GPy.kern.white(X.shape[1])
kernel1 = GPy.kern.rbf(X.shape[1]) #+ GPy.kern.white(X.shape[1])
kernel2 = kernel1.copy()
kernel3 = kernel1.copy()
kernel4 = kernel1.copy()
@ -90,12 +138,11 @@ def debug_student_t_noise_approx():
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['white'] = 1e-3
#m.constrain_positive('rbf')
#m.constrain_fixed('rbf_v', 1.0898)
#m.constrain_fixed('rbf_l', 1.8651)
#m['white'] = 1e-3
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('rbf')
m.constrain_positive('t_noise')
#m.constrain_fixed('t_noi', real_sd)
m.ensure_default_constraints()

47
GPy/likelihoods/Laplace.py
View file

@ -70,54 +70,38 @@ class Laplace(likelihood):
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)
Wi_K_i = self.W_12*self.Bi*self.W_12.T #same as rasms R
I_KW_i = np.eye(self.N) - np.dot(self.K, Wi_K_i)
return dL_dfhat, I_KW_i, Wi_K_i
I_KW_i = np.eye(self.N) - np.dot(self.K, self.Wi_K_i)
return dL_dfhat, I_KW_i
def _Kgradients(self, dK_dthetaK, X):
"""
Gradients with respect to prior kernel parameters
"""
dL_dfhat, I_KW_i, Wi_K_i = self._shared_gradients_components()
dL_dfhat, I_KW_i = self._shared_gradients_components()
dlp = self.likelihood_function.dlik_df(self.data, self.f_hat)
#Implicit
impl = mdot(dlp, dL_dfhat.T, I_KW_i)
expl_a = mdot(self.Ki_f, self.Ki_f.T)
expl_b = Wi_K_i
expl_b = self.Wi_K_i
expl = 0.5*expl_a + 0.5*expl_b
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)
dL_dthetaK = dL_dthetaK_imp + dL_dthetaK_exp
#dL_dthetaK = np.zeros(dK_dthetaK.shape)
#for thetaK_i, dK_dthetaK_i in enumerate(dK_dthetaK):
##Explicit
#f_Ki_dK_dtheta_Ki_f = mdot(self.Ki_f.T, dK_dthetaK_i, self.Ki_f)
#dL_dthetaK[thetaK_i] = 0.5*f_Ki_dK_dtheta_Ki_f - 0.5*np.trace(Wi_K_i*dK_dthetaK_i)
##Implicit
#df_hat_dthetaK = mdot(I_KW_i, dK_dthetaK_i, dlp)
#dL_dthetaK[thetaK_i] += np.dot(dL_dfhat.T, df_hat_dthetaK)
return dL_dthetaK
def _gradients(self, partial):
"""
Gradients with respect to likelihood parameters
"""
#return np.zeros(1)
dL_dfhat, I_KW_i, Wi_K_i = self._shared_gradients_components()
dL_dfhat, I_KW_i = self._shared_gradients_components()
dlik_dthetaL, dlik_grad_dthetaL, dlik_hess_dthetaL = self.likelihood_function._gradients(self.data, self.f_hat)
num_params = len(dlik_dthetaL)
dL_dthetaL = np.zeros(num_params) # make space for one derivative for each likelihood parameter
for thetaL_i in range(num_params):
#Explicit
#dL_dthetaL[thetaL_i] = np.sum(dlik_dthetaL[thetaL_i]) - 0.5*np.trace(np.dot(Ki_W_i.T, np.diagflat(dlik_hess_dthetaL[thetaL_i])))
#dL_dthetaL[thetaL_i] = np.sum(dlik_dthetaL[thetaL_i]) + 0.5*np.dot(Ki_W_i.T, dlik_hess_dthetaL[thetaL_i][:, None])
# might be +
dL_dthetaL_exp = np.sum(dlik_dthetaL[thetaL_i]) - 0.5*np.dot(np.diag(self.Ki_W_i), dlik_hess_dthetaL[thetaL_i])
#Implicit
df_hat_dthetaL = mdot(I_KW_i, self.K, dlik_grad_dthetaL[thetaL_i])
@ -165,34 +149,17 @@ class Laplace(likelihood):
Y_tilde = np.dot(self.Wi__Ki_W, self.f_hat)
ln_W_det = det_ln_diag(self.W)
yf_W_yf = mdot((Y_tilde - self.f_hat).T, np.diagflat(self.W), (Y_tilde - self.f_hat))
#Z_tilde = (+ self.NORMAL_CONST
#+ self.ln_z_hat
#+ 0.5*self.ln_I_KW_det
#- 0.5*ln_W_det
#+ 0.5*self.f_Ki_f
#+ 0.5*yf_W_yf
#)
self.Sigma_tilde = np.diagflat(1.0/self.W)
Ki, _, _, K_det = pdinv(self.K)
self.Wi_K_i = self.W_12*self.Bi*self.W_12.T #same as rasms R
ln_det_K_Wi__Bi = self.ln_I_KW_det + pddet(self.Sigma_tilde + self.K)
W = np.diagflat(self.W)
Wi = self.Sigma_tilde
W12i = np.sqrt(Wi)
#D = Ki - mdot((Ki + W), W12i, self.Bi, W12i, (Ki + W))
#fDf = mdot(self.f_hat.T, D, self.f_hat)
l = self.likelihood_function.link_function(self.data, self.f_hat, extra_data=self.extra_data)
#print "fDf:{} l:{} detKWiBi:{} W:{} Wi:{} Bi:{} Ki:{}".format(fDf, l, ln_det_K_Wi__Bi, W.sum(), Wi.sum(), self.Bi.sum(), Ki.sum())
y_Wi_Ki_i_y = mdot(Y_tilde.T, pdinv(self.K + Wi)[0], Y_tilde)
y_Wi_Ki_i_y = mdot(Y_tilde.T, self.Wi_K_i, Y_tilde)
Z_tilde = (+ self.NORMAL_CONST
+ l
+ 0.5*ln_det_K_Wi__Bi
#- 0.5*fDf
- 0.5*self.f_Ki_f
+ 0.5*y_Wi_Ki_i_y
)

7
GPy/likelihoods/likelihood_functions.py
View file

@ -194,10 +194,10 @@ class student_t(likelihood_function):
assert y.shape == f.shape
e = y - f
objective = (gammaln((self.v + 1) * 0.5)
objective = (+ gammaln((self.v + 1) * 0.5)
- gammaln(self.v * 0.5)
- np.log(self.sigma * np.sqrt(self.v * np.pi))
- (self.v + 1) * 0.5 * np.log(1 + ((e**2 / self.sigma**2) / self.v))
- (self.v + 1) * 0.5 * np.log(1 + (((e / self.sigma)**2) / self.v))
)
return np.sum(objective)
@ -234,7 +234,6 @@ class student_t(likelihood_function):
:returns: array which is diagonal of covariance matrix (second derivative of likelihood evaluated at points)
"""
assert y.shape == f.shape
e = y - f
hess = ((self.v + 1)*(e**2 - self.v*(self.sigma**2))) / ((((self.sigma**2)*self.v) + e**2)**2)
return hess
@ -247,7 +246,7 @@ class student_t(likelihood_function):
"""
assert y.shape == f.shape
e = y - f
d3lik_d3f = ( (2*(self.v + 1)*(-e)*(e**2 - 3*self.v*(self.sigma**2))) /
d3lik_d3f = ( -(2*(self.v + 1)*(-e)*(e**2 - 3*self.v*(self.sigma**2))) /
((e**2 + (self.sigma**2)*self.v)**3)
)
return d3lik_d3f