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

This commit is contained in:
Ricardo Andrade 2013-03-07 14:02:06 +00:00
commit 6350ae7afc
2 changed files with 16 additions and 15 deletions

View file

@ -123,7 +123,7 @@ def coregionalisation_toy():
Y = np.vstack((Y1,Y2)) Y = np.vstack((Y1,Y2))
k1 = GPy.kern.rbf(1) k1 = GPy.kern.rbf(1)
k2 = GPy.kern.coregionalise(2,1) k2 = GPy.kern.coregionalise(2,2)
k = k1.prod_orthogonal(k2) k = k1.prod_orthogonal(k2)
m = GPy.models.GP_regression(X,Y,kernel=k) m = GPy.models.GP_regression(X,Y,kernel=k)
m.constrain_fixed('rbf_var',1.) m.constrain_fixed('rbf_var',1.)
@ -147,8 +147,8 @@ def coregionalisation_sparse():
""" """
A simple demonstration of coregionalisation on two sinusoidal functions A simple demonstration of coregionalisation on two sinusoidal functions
""" """
X1 = np.random.rand(500,1)*8 X1 = np.random.rand(50,1)*8
X2 = np.random.rand(300,1)*5 X2 = np.random.rand(30,1)*5
index = np.vstack((np.zeros_like(X1),np.ones_like(X2))) index = np.vstack((np.zeros_like(X1),np.ones_like(X2)))
X = np.hstack((np.vstack((X1,X2)),index)) X = np.hstack((np.vstack((X1,X2)),index))
Y1 = np.sin(X1) + np.random.randn(*X1.shape)*0.05 Y1 = np.sin(X1) + np.random.randn(*X1.shape)*0.05
@ -158,7 +158,7 @@ def coregionalisation_sparse():
Z = np.hstack((np.random.rand(25,1)*8,np.random.randint(0,2,25)[:,None])) Z = np.hstack((np.random.rand(25,1)*8,np.random.randint(0,2,25)[:,None]))
k1 = GPy.kern.rbf(1) k1 = GPy.kern.rbf(1)
k2 = GPy.kern.coregionalise(2,2) k2 = GPy.kern.coregionalise(2,1)
k = k1.prod_orthogonal(k2) + GPy.kern.white(2,0.001) k = k1.prod_orthogonal(k2) + GPy.kern.white(2,0.001)
m = GPy.models.sparse_GP_regression(X,Y,kernel=k,Z=Z) m = GPy.models.sparse_GP_regression(X,Y,kernel=k,Z=Z)
@ -180,7 +180,6 @@ def coregionalisation_sparse():
return m return m
def multiple_optima(gene_number=937,resolution=80, model_restarts=10, seed=10000): 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.""" """Show an example of a multimodal error surface for Gaussian process regression. Gene 939 has bimodal behaviour where the noisey mode is higher."""

View file

@ -4,6 +4,7 @@
from kernpart import kernpart from kernpart import kernpart
import numpy as np import numpy as np
from GPy.util.linalg import mdot, pdinv from GPy.util.linalg import mdot, pdinv
import pdb
class coregionalise(kernpart): class coregionalise(kernpart):
""" """
@ -37,7 +38,6 @@ class coregionalise(kernpart):
self.B = np.dot(self.W,self.W.T) + np.diag(self.kappa) self.B = np.dot(self.W,self.W.T) + np.diag(self.kappa)
def _get_param_names(self): 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)] 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): def K(self,index,index2,target):
@ -46,7 +46,8 @@ class coregionalise(kernpart):
index2 = index index2 = index
else: else:
index2 = np.asarray(index2,dtype=np.int) index2 = np.asarray(index2,dtype=np.int)
ii,jj = np.meshgrid(index2,index) ii,jj = np.meshgrid(index,index2)
ii,jj = ii.T, jj.T
target += self.B[ii,jj] target += self.B[ii,jj]
def Kdiag(self,index,target): def Kdiag(self,index,target):
@ -58,19 +59,22 @@ class coregionalise(kernpart):
index2 = index index2 = index
else: else:
index2 = np.asarray(index2,dtype=np.int) index2 = np.asarray(index2,dtype=np.int)
ii,jj = np.meshgrid(index2,index) ii,jj = np.meshgrid(index,index2)
PK = np.zeros((self.R,self.R)) ii,jj = ii.T, jj.T
partial_small = np.zeros_like(self.B) partial_small = np.zeros_like(self.B)
for i in range(self.Nout): for i in range(self.Nout):
for j in range(self.Nout): for j in range(i,self.Nout):
partial_small[i,j] = np.sum(partial[(ii==i)*(jj==j)]) tmp = np.sum(partial[(ii==i)*(jj==j)])
dkappa = np.diag(partial_small) partial_small[i,j] = tmp
partial_small[j,i] = tmp
dkappa = np.diag(partial_small)
dW = 2.*(self.W[:,None,:]*partial_small[:,:,None]).sum(0) dW = 2.*(self.W[:,None,:]*partial_small[:,:,None]).sum(0)
target += np.hstack([dW.flatten(),dkappa]) target += np.hstack([dW.flatten(),dkappa])
def dKdiag_dtheta_foo(self,partial,index,target): def dKdiag_dtheta(self,partial,index,target):
index = np.asarray(index,dtype=np.int).flatten() index = np.asarray(index,dtype=np.int).flatten()
partial_small = np.zeros(self.Nout) partial_small = np.zeros(self.Nout)
for i in range(self.Nout): for i in range(self.Nout):
@ -82,7 +86,5 @@ class coregionalise(kernpart):
def dK_dX(self,partial,X,X2,target): def dK_dX(self,partial,X,X2,target):
pass pass
def dKdiag_dtheta(self,partial,index,target):
self.dK_dtheta(np.diag(partial),index,index,target)