apunkt/GPy
1
0
Fork 0
mirror of https://github.com/SheffieldML/GPy.git synced 2026-05-18 13:55:14 +02:00
Code Issues Projects Releases Packages Wiki Activity Actions

Merge branch 'master' of github.com:SheffieldML/GPy

Browse source
This commit is contained in:
Nicolò Fusi 2013-02-20 10:08:09 +00:00
commit 0f9eac5a25
2 changed files with 33 additions and 31 deletions
Download patch file Download diff file Expand all files Collapse all files

23
GPy/likelihoods/Gaussian.py
View file

@ -5,20 +5,26 @@ class Gaussian(likelihood):
def __init__(self,data,variance=1.,normalize=False): def __init__(self,data,variance=1.,normalize=False):
self.is_heteroscedastic = False self.is_heteroscedastic = False
self.Nparams = 1 self.Nparams = 1
self.data = data
self.N,D = data.shape
self.Z = 0. # a correction factor which accounts for the approximation made self.Z = 0. # a correction factor which accounts for the approximation made
N, self.D = data.shape
#normalisation #normalisation
if normalize: if normalize:
self._mean = data.mean(0)[None,:] self._mean = data.mean(0)[None,:]
self._std = data.std(0)[None,:] self._std = data.std(0)[None,:]
self.Y = (self.data - self._mean)/self._std
else: else:
self._mean = np.zeros((1,D)) self._mean = np.zeros((1,self.D))
self._std = np.ones((1,D)) self._std = np.ones((1,self.D))
self.Y = self.data
self.set_data(data)
self._set_params(np.asarray(variance))
def set_data(self,data):
self.data = data
self.N,D = data.shape
assert D == self.D
self.Y = (self.data - self._mean)/self._std
if D > self.N: if D > self.N:
self.YYT = np.dot(self.Y,self.Y.T) self.YYT = np.dot(self.Y,self.Y.T)
self.trYYT = np.trace(self.YYT) self.trYYT = np.trace(self.YYT)
@ -26,14 +32,11 @@ class Gaussian(likelihood):
self.YYT = None self.YYT = None
self.trYYT = None self.trYYT = None
self._set_params(np.asarray(variance))
def _get_params(self): def _get_params(self):
return np.asarray(self._variance) return np.asarray(self._variance)
def _get_param_names(self): def _get_param_names(self):
return ["noise variance"] return ["noise_variance"]
def _set_params(self,x): def _set_params(self,x):
self._variance = float(x) self._variance = float(x)

41
GPy/models/uncollapsed_sparse_GP.py
View file

@ -27,10 +27,10 @@ class uncollapsed_sparse_GP(sparse_GP):
""" """
def __init__(self, X, likelihood, kernel, Z, q_u=None, **kwargs): def __init__(self, X, likelihood, kernel, Z, q_u=None, **kwargs):
self.D = Y.shape[1] self.M = Z.shape[0]
self.M = kwargs['Z'].shape[0]
if q_u is None: if q_u is None:
q_u = np.hstack((np.random.randn(self.M*self.D),-0.5*np.eye(self.M).flatten())) q_u = np.hstack((np.random.randn(self.M*likelihood.D),-0.5*np.eye(self.M).flatten()))
self.likelihood = likelihood
self.set_vb_param(q_u) self.set_vb_param(q_u)
sparse_GP.__init__(self, X, likelihood, kernel, Z, **kwargs) sparse_GP.__init__(self, X, likelihood, kernel, Z, **kwargs)
@ -62,21 +62,21 @@ class uncollapsed_sparse_GP(sparse_GP):
self.projected_mean = mdot(self.psi1.T,self.Kmmi,self.q_u_expectation[0]) self.projected_mean = mdot(self.psi1.T,self.Kmmi,self.q_u_expectation[0])
# Compute dL_dpsi # Compute dL_dpsi
self.dL_dpsi0 = - 0.5 * self.D * self.beta * np.ones(self.N) self.dL_dpsi0 = - 0.5 * self.likelihood.D * self.beta * np.ones(self.N)
self.dL_dpsi1 = np.dot(self.VmT,self.Kmmi).T # This is the correct term for E I think... self.dL_dpsi1 = np.dot(self.VmT,self.Kmmi).T # This is the correct term for E I think...
self.dL_dpsi2 = 0.5 * self.beta * self.D * (self.Kmmi - mdot(self.Kmmi,self.q_u_expectation[1],self.Kmmi)) self.dL_dpsi2 = 0.5 * self.beta * self.likelihood.D * (self.Kmmi - mdot(self.Kmmi,self.q_u_expectation[1],self.Kmmi))
# Compute dL_dKmm # Compute dL_dKmm
tmp = self.beta*mdot(self.psi2,self.Kmmi,self.q_u_expectation[1]) -np.dot(self.q_u_expectation[0],self.psi1V.T) tmp = self.beta*mdot(self.psi2,self.Kmmi,self.q_u_expectation[1]) -np.dot(self.q_u_expectation[0],self.psi1V.T)
tmp += tmp.T tmp += tmp.T
tmp += self.D*(-self.beta*self.psi2 - self.Kmm + self.q_u_expectation[1]) tmp += self.likelihood.D*(-self.beta*self.psi2 - self.Kmm + self.q_u_expectation[1])
self.dL_dKmm = 0.5*mdot(self.Kmmi,tmp,self.Kmmi) self.dL_dKmm = 0.5*mdot(self.Kmmi,tmp,self.Kmmi)
#Compute the gradient of the log likelihood wrt noise variance #Compute the gradient of the log likelihood wrt noise variance
#TODO: suport heteroscedatic noise #TODO: suport heteroscedatic noise
dbeta = 0.5 * self.N*self.D/self.beta dbeta = 0.5 * self.N*self.likelihood.D/self.beta
dbeta += - 0.5 * self.D * self.trace_K dbeta += - 0.5 * self.likelihood.D * self.trace_K
dbeta += - 0.5 * self.D * np.sum(self.q_u_expectation[1]*mdot(self.Kmmi,self.psi2,self.Kmmi)) dbeta += - 0.5 * self.likelihood.D * np.sum(self.q_u_expectation[1]*mdot(self.Kmmi,self.psi2,self.Kmmi))
dbeta += - 0.5 * self.trYYT dbeta += - 0.5 * self.trYYT
dbeta += np.sum(np.dot(self.Y.T,self.projected_mean)) dbeta += np.sum(np.dot(self.Y.T,self.projected_mean))
self.partial_for_likelihood = -dbeta*self.likelihood.precision**2 self.partial_for_likelihood = -dbeta*self.likelihood.precision**2
@ -85,9 +85,9 @@ class uncollapsed_sparse_GP(sparse_GP):
""" """
Compute the (lower bound on the) log marginal likelihood Compute the (lower bound on the) log marginal likelihood
""" """
A = -0.5*self.N*self.D*(np.log(2.*np.pi) - np.log(self.beta)) A = -0.5*self.N*self.likelihood.D*(np.log(2.*np.pi) - np.log(self.beta))
B = -0.5*self.beta*self.D*self.trace_K B = -0.5*self.beta*self.likelihood.D*self.trace_K
C = -0.5*self.D *(self.Kmm_logdet-self.q_u_logdet + np.sum(self.Lambda * self.q_u_expectation[1]) - self.M) C = -0.5*self.likelihood.D *(self.Kmm_logdet-self.q_u_logdet + np.sum(self.Lambda * self.q_u_expectation[1]) - self.M)
D = -0.5*self.beta*self.trYYT D = -0.5*self.beta*self.trYYT
E = np.sum(np.dot(self.V.T,self.projected_mean)) E = np.sum(np.dot(self.V.T,self.projected_mean))
return A+B+C+D+E return A+B+C+D+E
@ -100,21 +100,21 @@ class uncollapsed_sparse_GP(sparse_GP):
tmp = self.Kmmi- mdot(self.Kmmi,self.q_u_cov,self.Kmmi) tmp = self.Kmmi- mdot(self.Kmmi,self.q_u_cov,self.Kmmi)
if full_cov: if full_cov:
Kxx = self.kern.K(Xnew) Kxx = self.kern.K(Xnew)
var = Kxx - mdot(Kx,tmp,Kx.T) + np.eye(Xnew.shape[0])/self.beta var = Kxx - mdot(Kx,tmp,Kx.T)
else: else:
Kxx = self.kern.Kdiag(Xnew) Kxx = self.kern.Kdiag(Xnew)
var = Kxx - np.sum(Kx*np.dot(Kx,tmp),1) + 1./self.beta var = (Kxx - np.sum(Kx*np.dot(Kx,tmp),1))[:,None]
return mu,var return mu,var
def set_vb_param(self,vb_param): def set_vb_param(self,vb_param):
"""set the distribution q(u) from the canonical parameters""" """set the distribution q(u) from the canonical parameters"""
self.q_u_prec = -2.*vb_param[self.M*self.D:].reshape(self.M,self.M) self.q_u_prec = -2.*vb_param[-self.M**2:].reshape(self.M, self.M)
self.q_u_cov, q_u_Li, q_u_L, tmp = pdinv(self.q_u_prec) self.q_u_cov, q_u_Li, q_u_L, tmp = pdinv(self.q_u_prec)
self.q_u_logdet = -tmp self.q_u_logdet = -tmp
self.q_u_mean = np.dot(self.q_u_cov,vb_param[:self.M*self.D].reshape(self.M,self.D)) self.q_u_mean = np.dot(self.q_u_cov,vb_param[:self.M*self.likelihood.D].reshape(self.M,self.likelihood.D))
self.q_u_expectation = (self.q_u_mean, np.dot(self.q_u_mean,self.q_u_mean.T)+self.q_u_cov) self.q_u_expectation = (self.q_u_mean, np.dot(self.q_u_mean,self.q_u_mean.T)+self.q_u_cov*self.likelihood.D)
self.q_u_canonical = (np.dot(self.q_u_prec, self.q_u_mean),-0.5*self.q_u_prec) self.q_u_canonical = (np.dot(self.q_u_prec, self.q_u_mean),-0.5*self.q_u_prec)
#TODO: computations now? #TODO: computations now?
@ -133,8 +133,7 @@ class uncollapsed_sparse_GP(sparse_GP):
Note that the natural gradient in either is given by the gradient in the other (See Hensman et al 2012 Fast Variational inference in the conjugate exponential Family) Note that the natural gradient in either is given by the gradient in the other (See Hensman et al 2012 Fast Variational inference in the conjugate exponential Family)
""" """
dL_dmmT_S = -0.5*self.Lambda-self.q_u_canonical[1] dL_dmmT_S = -0.5*self.Lambda-self.q_u_canonical[1]
#dL_dm = np.dot(self.Kmmi,self.psi1V) - np.dot(self.Lambda,self.q_u_mean) dL_dm = np.dot(self.Kmmi,self.psi1V) - np.dot(self.Lambda,self.q_u_mean)
dL_dm = np.dot(self.Kmmi,self.psi1V) - self.q_u_canonical[0]
#dL_dSim = #dL_dSim =
#dL_dmhSi = #dL_dmhSi =
@ -144,9 +143,9 @@ class uncollapsed_sparse_GP(sparse_GP):
def plot(self, *args, **kwargs): def plot(self, *args, **kwargs):
""" """
add the distribution q(u) to the plot from sparse_GP_regression add the distribution q(u) to the plot from sparse_GP
""" """
sparse_GP_regression.plot(self,*args,**kwargs) sparse_GP.plot(self,*args,**kwargs)
if self.Q==1: if self.Q==1:
pb.errorbar(self.Z[:,0],self.q_u_expectation[0][:,0],yerr=2.*np.sqrt(np.diag(self.q_u_cov)),fmt=None,ecolor='b') pb.errorbar(self.Z[:,0],self.q_u_expectation[0][:,0],yerr=2.*np.sqrt(np.diag(self.q_u_cov)),fmt=None,ecolor='b')