diff --git a/GPy/examples/regression.py b/GPy/examples/regression.py index b8838078..7d092c26 100644 --- a/GPy/examples/regression.py +++ b/GPy/examples/regression.py @@ -143,6 +143,43 @@ def coregionalisation_toy(): return m +def coregionalisation_sparse(): + """ + A simple demonstration of coregionalisation on two sinusoidal functions + """ + X1 = np.random.rand(500,1)*8 + X2 = np.random.rand(300,1)*5 + index = np.vstack((np.zeros_like(X1),np.ones_like(X2))) + X = np.hstack((np.vstack((X1,X2)),index)) + Y1 = np.sin(X1) + np.random.randn(*X1.shape)*0.05 + Y2 = -np.sin(X2) + np.random.randn(*X2.shape)*0.05 + Y = np.vstack((Y1,Y2)) + + Z = np.hstack((np.random.rand(25,1)*8,np.random.randint(0,2,25)[:,None])) + + k1 = GPy.kern.rbf(1) + k2 = GPy.kern.coregionalise(2,2) + k = k1.prod_orthogonal(k2) + GPy.kern.white(2,0.001) + + m = GPy.models.sparse_GP_regression(X,Y,kernel=k,Z=Z) + m.constrain_fixed('rbf_var',1.) + m.constrain_positive('kappa') + m.constrain_fixed('iip') + m.ensure_default_constraints() + #m.optimize() + + pb.figure() + Xtest1 = np.hstack((np.linspace(0,9,100)[:,None],np.zeros((100,1)))) + Xtest2 = np.hstack((np.linspace(0,9,100)[:,None],np.ones((100,1)))) + mean, var,low,up = m.predict(Xtest1) + GPy.util.plot.gpplot(Xtest1[:,0],mean,low,up) + mean, var,low,up = m.predict(Xtest2) + GPy.util.plot.gpplot(Xtest2[:,0],mean,low,up) + pb.plot(X1[:,0],Y1[:,0],'rx',mew=2) + pb.plot(X2[:,0],Y2[:,0],'gx',mew=2) + return m + + def multiple_optima(gene_number=937,resolution=80, model_restarts=10, seed=10000): """Show an example of a multimodal error surface for Gaussian process regression. Gene 939 has bimodal behaviour where the noisey mode is higher.""" diff --git a/GPy/kern/coregionalise.py b/GPy/kern/coregionalise.py index b8fdbf42..29073b3a 100644 --- a/GPy/kern/coregionalise.py +++ b/GPy/kern/coregionalise.py @@ -46,8 +46,8 @@ class coregionalise(kernpart): index2 = index else: index2 = np.asarray(index2,dtype=np.int) - ii,jj = np.meshgrid(index,index2) - target += self.B[ii,jj].T + ii,jj = np.meshgrid(index2,index) + target += self.B[ii,jj] def Kdiag(self,index,target): target += np.diag(self.B)[np.asarray(index,dtype=np.int).flatten()] @@ -58,26 +58,47 @@ class coregionalise(kernpart): index2 = index else: index2 = np.asarray(index2,dtype=np.int) - ii,jj = np.meshgrid(index,index2) + ii,jj = np.meshgrid(index2,index) PK = np.zeros((self.R,self.R)) - dkappa = np.zeros(self.Nout) partial_small = np.zeros_like(self.B) for i in range(self.Nout): for j in range(self.Nout): - partial_small[j,i] = np.sum(partial[(ii==i)*(jj==j)]) - #print partial_small + partial_small[i,j] = np.sum(partial[(ii==i)*(jj==j)]) dkappa = np.diag(partial_small) - ##target += (((X2[:, None, :] * self.variances)) * partial[:,:, None]).sum(0) dW = 2.*(self.W[:,None,:]*partial_small[:,:,None]).sum(0) target += np.hstack([dW.flatten(),dkappa]) def dKdiag_dtheta(self,partial,index,target): - raise NotImplementedError + index = np.asarray(index,dtype=np.int).flatten() + partial_small = np.zeros(self.Nout) + for i in range(self.Nout): + partial_small[i] += np.sum(partial[index==i]) + dW = 2.*self.W*partial_small[:,None] + dkappa = partial_small + target += np.hstack([dW.flatten(),dkappa]) def dK_dX(self,partial,X,X2,target): pass + def dKdiag_dthetai_(self,partial,index,target): + index = np.asarray(index,dtype=np.int) + index2 = index + ii,jj = np.meshgrid(index2,index) + PK = np.zeros((self.R,self.R)) + partial_small = np.zeros_like(self.B) + for i in range(self.Nout): + for j in range(self.Nout): + partial_small[j,i] = np.sum(partial[np.diag((ii==i)*(jj==j))]) + #print partial_small + dkappa = np.diag(partial_small) + + ##target += (((X2[:, None, :] * self.variances)) * partial[:,:, None]).sum(0) + partial_small = np.diag(np.diag(partial_small)) + #dW = 2.*(self.W[:,None,:]*partial_small[:,:,None]).sum(0) + dW = 2. + + target += np.hstack([dW.flatten(),dkappa]) diff --git a/GPy/kern/product_orthogonal.py b/GPy/kern/product_orthogonal.py index a729c126..6b02b868 100644 --- a/GPy/kern/product_orthogonal.py +++ b/GPy/kern/product_orthogonal.py @@ -46,14 +46,6 @@ class product_orthogonal(kernpart): self.k2.K(X[:,self.k1.D:],X2[:,self.k1.D:],target2) target += target1 * target2 - def Kdiag(self,X,target): - """Compute the diagonal of the covariance matrix associated to X.""" - target1 = np.zeros((X.shape[0],)) - target2 = np.zeros((X.shape[0],)) - self.k1.Kdiag(X[:,0:self.k1.D],target1) - self.k2.Kdiag(X[:,self.k1.D:],target2) - target += target1 * target2 - def dK_dtheta(self,partial,X,X2,target): """derivative of the covariance matrix with respect to the parameters.""" if X2 is None: X2 = X @@ -70,6 +62,22 @@ class product_orthogonal(kernpart): target[:self.k1.Nparam] += k1_target target[self.k1.Nparam:] += k2_target + def Kdiag(self,X,target): + """Compute the diagonal of the covariance matrix associated to X.""" + target1 = np.zeros((X.shape[0],)) + target2 = np.zeros((X.shape[0],)) + self.k1.Kdiag(X[:,:self.k1.D],target1) + self.k2.Kdiag(X[:,self.k1.D:],target2) + target += target1 * target2 + + def dKdiag_dtheta(self,partial,X,target): + K1 = np.zeros(X.shape[0]) + K2 = np.zeros(X.shape[0]) + self.k1.Kdiag(X[:,:self.k1.D],K1) + self.k2.Kdiag(X[:,self.k1.D:],K2) + self.k1.dKdiag_dtheta(partial*K2,X[:,:self.k1.D],target[:self.k1.Nparam]) + self.k2.dKdiag_dtheta(partial*K1,X[:,self.k1.D:],target[self.k1.Nparam:]) + def dK_dX(self,partial,X,X2,target): """derivative of the covariance matrix with respect to X.""" if X2 is None: X2 = X