diff --git a/GPy/examples/regression.py b/GPy/examples/regression.py index 43fa0147..b8838078 100644 --- a/GPy/examples/regression.py +++ b/GPy/examples/regression.py @@ -75,6 +75,74 @@ def silhouette(): print(m) return m +def coregionalisation_toy2(): + """ + A simple demonstration of coregionalisation on two sinusoidal functions + """ + X1 = np.random.rand(50,1)*8 + X2 = np.random.rand(30,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 + 2. + Y = np.vstack((Y1,Y2)) + + k1 = GPy.kern.rbf(1) + GPy.kern.bias(1) + k2 = GPy.kern.coregionalise(2,1) + k = k1.prod_orthogonal(k2) + m = GPy.models.GP_regression(X,Y,kernel=k) + m.constrain_fixed('rbf_var',1.) + m.constrain_positive('kappa') + 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 coregionalisation_toy(): + """ + A simple demonstration of coregionalisation on two sinusoidal functions + """ + X1 = np.random.rand(50,1)*8 + X2 = np.random.rand(30,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)) + + k1 = GPy.kern.rbf(1) + k2 = GPy.kern.coregionalise(2,1) + k = k1.prod_orthogonal(k2) + m = GPy.models.GP_regression(X,Y,kernel=k) + m.constrain_fixed('rbf_var',1.) + m.constrain_positive('kappa') + 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.""" @@ -160,3 +228,4 @@ def contour_data(data, length_scales, log_SNRs, signal_kernel_call=GPy.kern.rbf) length_scale_lls.append(model.log_likelihood()) lls.append(length_scale_lls) return np.array(lls) + diff --git a/GPy/kern/__init__.py b/GPy/kern/__init__.py index 3e20644b..625f6080 100644 --- a/GPy/kern/__init__.py +++ b/GPy/kern/__init__.py @@ -2,5 +2,5 @@ # Licensed under the BSD 3-clause license (see LICENSE.txt) -from constructors import rbf, Matern32, Matern52, exponential, linear, white, bias, finite_dimensional, spline, Brownian, rbf_sympy, sympykern, periodic_exponential, periodic_Matern32, periodic_Matern52, product, product_orthogonal, symmetric +from constructors import rbf, Matern32, Matern52, exponential, linear, white, bias, finite_dimensional, spline, Brownian, rbf_sympy, sympykern, periodic_exponential, periodic_Matern32, periodic_Matern52, product, product_orthogonal, symmetric, coregionalise from kern import kern diff --git a/GPy/kern/constructors.py b/GPy/kern/constructors.py index 5d665d93..9b58c282 100644 --- a/GPy/kern/constructors.py +++ b/GPy/kern/constructors.py @@ -21,6 +21,7 @@ from periodic_Matern52 import periodic_Matern52 as periodic_Matern52part from product import product as productpart from product_orthogonal import product_orthogonal as product_orthogonalpart from symmetric import symmetric as symmetric_part +from coregionalise import coregionalise as coregionalise_part #TODO these s=constructors are not as clean as we'd like. Tidy the code up #using meta-classes to make the objects construct properly wthout them. @@ -274,3 +275,8 @@ def symmetric(k): k_.parts = [symmetric_part(p) for p in k.parts] return k_ +def coregionalise(Nout,R=1, W=None, kappa=None): + p = coregionalise_part(Nout,R,W,kappa) + return kern(1,[p]) + + diff --git a/GPy/kern/coregionalise.py b/GPy/kern/coregionalise.py new file mode 100644 index 00000000..b8fdbf42 --- /dev/null +++ b/GPy/kern/coregionalise.py @@ -0,0 +1,83 @@ +# Copyright (c) 2012, James Hensman and Ricardo Andrade +# Licensed under the BSD 3-clause license (see LICENSE.txt) + +from kernpart import kernpart +import numpy as np +from GPy.util.linalg import mdot, pdinv + +class coregionalise(kernpart): + """ + Kernel for Intrisec Corregionalization Models + """ + def __init__(self,Nout,R=1, W=None, kappa=None): + self.D = 1 + self.name = 'coregion' + self.Nout = Nout + self.R = R + if W is None: + self.W = np.ones((self.Nout,self.R)) + else: + assert W.shape==(self.Nout,self.R) + self.W = W + if kappa is None: + kappa = np.ones(self.Nout) + else: + assert kappa.shape==(self.Nout,) + self.kappa = kappa + self.Nparam = self.Nout*(self.R + 1) + self._set_params(np.hstack([self.W.flatten(),self.kappa])) + + def _get_params(self): + return np.hstack([self.W.flatten(),self.kappa]) + + def _set_params(self,x): + assert x.size == self.Nparam + self.kappa = x[-self.Nout:] + self.W = x[:-self.Nout].reshape(self.Nout,self.R) + self.B = np.dot(self.W,self.W.T) + np.diag(self.kappa) + + def _get_param_names(self): + + return sum([['W%i_%i'%(i,j) for j in range(self.R)] for i in range(self.Nout)],[]) + ['kappa_%i'%i for i in range(self.Nout)] + + def K(self,index,index2,target): + index = np.asarray(index,dtype=np.int) + if index2 is None: + index2 = index + else: + index2 = np.asarray(index2,dtype=np.int) + ii,jj = np.meshgrid(index,index2) + target += self.B[ii,jj].T + + def Kdiag(self,index,target): + target += np.diag(self.B)[np.asarray(index,dtype=np.int).flatten()] + + def dK_dtheta(self,partial,index,index2,target): + index = np.asarray(index,dtype=np.int) + if index2 is None: + index2 = index + else: + index2 = np.asarray(index2,dtype=np.int) + ii,jj = np.meshgrid(index,index2) + 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 + 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 + + def dK_dX(self,partial,X,X2,target): + pass + + +