From 73f7b72708c8da7efe6e5ba8d9f7259342ec6276 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Tue, 19 Feb 2013 14:44:39 +0000 Subject: [PATCH 1/4] added a set_data method to the Gaussian likelihood --- GPy/likelihoods/Gaussian.py | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/GPy/likelihoods/Gaussian.py b/GPy/likelihoods/Gaussian.py index 9a992ac1..ed48caf0 100644 --- a/GPy/likelihoods/Gaussian.py +++ b/GPy/likelihoods/Gaussian.py @@ -5,20 +5,24 @@ class Gaussian(likelihood): def __init__(self,data,variance=1.,normalize=False): self.is_heteroscedastic = False self.Nparams = 1 - self.data = data - self.N,D = data.shape self.Z = 0. # a correction factor which accounts for the approximation made #normalisation if normalize: self._mean = data.mean(0)[None,:] self._std = data.std(0)[None,:] - self.Y = (self.data - self._mean)/self._std else: self._mean = np.zeros((1,D)) self._std = np.ones((1,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 + self.Y = (self.data - self._mean)/self._std if D > self.N: self.YYT = np.dot(self.Y,self.Y.T) self.trYYT = np.trace(self.YYT) @@ -26,9 +30,6 @@ class Gaussian(likelihood): self.YYT = None self.trYYT = None - self._set_params(np.asarray(variance)) - - def _get_params(self): return np.asarray(self._variance) From 70e2076cd4eb04cb96fc1a1005f8a22570176f16 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Tue, 19 Feb 2013 14:48:30 +0000 Subject: [PATCH 2/4] bugfixin' --- GPy/likelihoods/Gaussian.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/GPy/likelihoods/Gaussian.py b/GPy/likelihoods/Gaussian.py index ed48caf0..38e1d816 100644 --- a/GPy/likelihoods/Gaussian.py +++ b/GPy/likelihoods/Gaussian.py @@ -6,14 +6,15 @@ class Gaussian(likelihood): self.is_heteroscedastic = False self.Nparams = 1 self.Z = 0. # a correction factor which accounts for the approximation made + N, self.D = data.shape #normalisation if normalize: self._mean = data.mean(0)[None,:] self._std = data.std(0)[None,:] else: - self._mean = np.zeros((1,D)) - self._std = np.ones((1,D)) + self._mean = np.zeros((1,self.D)) + self._std = np.ones((1,self.D)) self.set_data(data) @@ -22,6 +23,7 @@ class Gaussian(likelihood): 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: self.YYT = np.dot(self.Y,self.Y.T) From 735a3403158efa023b4c0e25c597462438fdad05 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Tue, 19 Feb 2013 17:36:51 +0000 Subject: [PATCH 3/4] changes to the uncollapsed GP --- GPy/models/uncollapsed_sparse_GP.py | 42 ++++++++++++++--------------- 1 file changed, 21 insertions(+), 21 deletions(-) diff --git a/GPy/models/uncollapsed_sparse_GP.py b/GPy/models/uncollapsed_sparse_GP.py index 21d5e5ff..ddeccd18 100644 --- a/GPy/models/uncollapsed_sparse_GP.py +++ b/GPy/models/uncollapsed_sparse_GP.py @@ -27,10 +27,10 @@ class uncollapsed_sparse_GP(sparse_GP): """ def __init__(self, X, likelihood, kernel, Z, q_u=None, **kwargs): - self.D = Y.shape[1] - self.M = kwargs['Z'].shape[0] + self.M = Z.shape[0] 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) 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]) # 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_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 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 += 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) #Compute the gradient of the log likelihood wrt noise variance #TODO: suport heteroscedatic noise - dbeta = 0.5 * self.N*self.D/self.beta - dbeta += - 0.5 * self.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.N*self.likelihood.D/self.beta + dbeta += - 0.5 * self.likelihood.D * self.trace_K + 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 += np.sum(np.dot(self.Y.T,self.projected_mean)) 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 """ - A = -0.5*self.N*self.D*(np.log(2.*np.pi) - np.log(self.beta)) - B = -0.5*self.beta*self.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) + A = -0.5*self.N*self.likelihood.D*(np.log(2.*np.pi) - np.log(self.beta)) + B = -0.5*self.beta*self.likelihood.D*self.trace_K + 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 E = np.sum(np.dot(self.V.T,self.projected_mean)) 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) if full_cov: 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: 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 def set_vb_param(self,vb_param): """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_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) #TODO: computations now? @@ -133,8 +133,8 @@ 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) """ 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) - self.q_u_canonical[0] + 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_dmhSi = @@ -144,9 +144,9 @@ class uncollapsed_sparse_GP(sparse_GP): 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: 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') From 37f6ed02998178301003c275319dc6c85fddb28c Mon Sep 17 00:00:00 2001 From: James Hensman Date: Wed, 20 Feb 2013 09:12:00 +0000 Subject: [PATCH 4/4] made the name of the Gaussian noise variance noise_variance, for consistency --- GPy/likelihoods/Gaussian.py | 2 +- GPy/models/uncollapsed_sparse_GP.py | 1 - 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/GPy/likelihoods/Gaussian.py b/GPy/likelihoods/Gaussian.py index 38e1d816..4a32f066 100644 --- a/GPy/likelihoods/Gaussian.py +++ b/GPy/likelihoods/Gaussian.py @@ -36,7 +36,7 @@ class Gaussian(likelihood): return np.asarray(self._variance) def _get_param_names(self): - return ["noise variance"] + return ["noise_variance"] def _set_params(self,x): self._variance = float(x) diff --git a/GPy/models/uncollapsed_sparse_GP.py b/GPy/models/uncollapsed_sparse_GP.py index ddeccd18..43624e72 100644 --- a/GPy/models/uncollapsed_sparse_GP.py +++ b/GPy/models/uncollapsed_sparse_GP.py @@ -134,7 +134,6 @@ class uncollapsed_sparse_GP(sparse_GP): """ 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) - self.q_u_canonical[0] #dL_dSim = #dL_dmhSi =