apunkt/GPy
1
0
Fork 0
mirror of https://github.com/SheffieldML/GPy.git synced 2026-05-15 06:52:39 +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:
James Hensman 2013-02-19 14:05:55 +00:00
commit 389db915d4
8 changed files with 134 additions and 31 deletions
Download patch file Download diff file Expand all files Collapse all files

14
GPy/examples/BGPLVM_demo.py
View file

@ -8,21 +8,19 @@ np.random.seed(123344)
N = 10
M = 3
Q = 4
D = 5
Q = 2
D = 4
#generate GPLVM-like data
X = np.random.rand(N, Q)
k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001)
K = k.K(X)
Y = np.random.multivariate_normal(np.zeros(N),K,D).T
<<<<<<< HEAD
k = GPy.kern.rbf(Q) + GPy.kern.rbf(Q) + GPy.kern.white(Q)
# k = GPy.kern.linear(Q, ARD = True) + GPy.kern.white(Q, 0.00001)
=======
k = GPy.kern.linear(Q, ARD = True) + GPy.kern.white(Q)
# k = GPy.kern.rbf(Q) + GPy.kern.rbf(Q) + GPy.kern.white(Q)
# k = GPy.kern.rbf(Q) + GPy.kern.bias(Q) + GPy.kern.white(Q, 0.00001)
k = GPy.kern.rbf(Q, ARD = False) + GPy.kern.white(Q, 0.00001)
>>>>>>> master
# k = GPy.kern.rbf(Q, ARD = False) + GPy.kern.white(Q, 0.00001)
m = GPy.models.Bayesian_GPLVM(Y, Q, kernel = k, M=M)
m.constrain_positive('(rbf|bias|noise|white|S)')
# m.constrain_fixed('S', 1)

4
GPy/examples/oil_flow_demo.py
View file

@ -32,7 +32,7 @@ Y -= Y.mean(axis=0)
# Y /= Y.std(axis=0)
Q = 5
k = GPy.kern.linear(Q, ARD = False) + GPy.kern.white(Q)
k = GPy.kern.linear(Q, ARD = True) + GPy.kern.white(Q)
m = GPy.models.Bayesian_GPLVM(Y, Q, kernel = k, M = 20)
m.constrain_positive('(rbf|bias|S|linear|white|noise)')
@ -43,7 +43,7 @@ m.constrain_positive('(rbf|bias|S|linear|white|noise)')
# plot_oil(m.X, np.array([1,1]), labels, 'PCA initialization')
m.optimize(messages = True)
# m.optimize('tnc', messages = True)
plot_oil(m.X, m.kern.parts[0].lengthscale, labels, 'B-GPLVM')
# plot_oil(m.X, m.kern.parts[0].lengthscale, labels, 'B-GPLVM')
# # pb.figure()
# m.plot()
# pb.title('PCA initialisation')

18
GPy/kern/bias.py
View file

@ -47,6 +47,10 @@ class bias(kernpart):
def dKdiag_dX(self,partial,X,target):
pass
#---------------------------------------#
# PSI statistics #
#---------------------------------------#
def psi0(self, Z, mu, S, target):
target += self.variance
@ -59,27 +63,27 @@ class bias(kernpart):
def dpsi0_dtheta(self, partial, Z, mu, S, target):
target += partial.sum()
def dpsi1_dtheta(self, partial, Z, mu, S, target):
target += partial.sum()
def dpsi2_dtheta(self, partial, Z, mu, S, target):
target += 2.*self.variance*partial.sum()
def dpsi0_dZ(self, partial, Z, mu, S, target):
pass
def dpsi0_dmuS(self, partial, Z, mu, S, target_mu, target_S):
pass
def dpsi1_dtheta(self, partial, Z, mu, S, target):
target += partial.sum()
def dpsi1_dZ(self, partial, Z, mu, S, target):
pass
def dpsi1_dmuS(self, partial, Z, mu, S, target_mu, target_S):
pass
def dpsi2_dtheta(self, partial, Z, mu, S, target):
target += np.sum(2.*self.variance*partial)
def dpsi2_dZ(self, partial, Z, mu, S, target):
pass
def dpsi2_dmuS(self, partial, Z, mu, S, target_mu, target_S):
pass

59
GPy/kern/kern.py
View file

@ -3,10 +3,11 @@
import numpy as np
import pylab as pb
from ..core.parameterised import parameterised
from kernpart import kernpart
import itertools
from product_orthogonal import product_orthogonal
from product_orthogonal import product_orthogonal
class kern(parameterised):
def __init__(self,D,parts=[], input_slices=None):
@ -386,3 +387,59 @@ class kern(parameterised):
#TODO: there are some extra terms to compute here!
return target_mu, target_S
def plot(self, x = None, plot_limits=None,which_functions='all',resolution=None,*args,**kwargs):
if which_functions=='all':
which_functions = [True]*self.Nparts
if self.D == 1:
if x is None:
x = np.zeros((1,1))
else:
x = np.asarray(x)
assert x.size == 1, "The size of the fixed variable x is not 1"
x = x.reshape((1,1))
if plot_limits == None:
xmin, xmax = (x-5).flatten(), (x+5).flatten()
elif len(plot_limits) == 2:
xmin, xmax = plot_limits
else:
raise ValueError, "Bad limits for plotting"
Xnew = np.linspace(xmin,xmax,resolution or 201)[:,None]
Kx = self.K(Xnew,x,slices2=which_functions)
pb.plot(Xnew,Kx,*args,**kwargs)
pb.xlim(xmin,xmax)
pb.xlabel("x")
pb.ylabel("k(x,%0.1f)" %x)
elif self.D == 2:
if x is None:
x = np.zeros((1,2))
else:
x = np.asarray(x)
assert x.size == 2, "The size of the fixed variable x is not 2"
x = x.reshape((1,2))
if plot_limits == None:
xmin, xmax = (x-5).flatten(), (x+5).flatten()
elif len(plot_limits) == 2:
xmin, xmax = plot_limits
else:
raise ValueError, "Bad limits for plotting"
resolution = resolution or 51
xx,yy = np.mgrid[xmin[0]:xmax[0]:1j*resolution,xmin[1]:xmax[1]:1j*resolution]
xg = np.linspace(xmin[0],xmax[0],resolution)
yg = np.linspace(xmin[1],xmax[1],resolution)
Xnew = np.vstack((xx.flatten(),yy.flatten())).T
Kx = self.K(Xnew,x,slices2=which_functions)
Kx = Kx.reshape(resolution,resolution).T
pb.contour(xg,yg,Kx,vmin=Kx.min(),vmax=Kx.max(),cmap=pb.cm.jet)
pb.xlim(xmin[0],xmax[0])
pb.ylim(xmin[1],xmax[1])
pb.xlabel("x1")
pb.ylabel("x2")
pb.title("k(x1,x2 ; %0.1f,%0.1f)" %(x[0,0],x[0,1]) )
else:
raise NotImplementedError, "Cannot plot a kernel with more than two input dimensions"

13
GPy/kern/linear.py
View file

@ -90,22 +90,19 @@ class linear(kernpart):
def psi0(self,Z,mu,S,target):
self._psi_computations(Z,mu,S)
target += np.sum(self.variances*self.mu2_S)
target += np.sum(self.variances*self.mu2_S,1)
def dpsi0_dtheta(self,partial,Z,mu,S,target):
self._psi_computations(Z,mu,S)
tmp = (partial[:, None] * (np.sum(self.mu2_S,0)))
tmp = partial[:, None] * self.mu2_S
if self.ARD:
target += tmp.sum(0)
else:
target += tmp.sum()
def dpsi0_dmuS(self,partial, Z,mu,S,target_mu,target_S):
target_mu += np.sum(partial[:, None],0) * (2.0*mu*self.variances)
target_S += np.sum(partial[:, None] * self.variances, 0)
def dpsi0_dZ(self,Z,mu,S,target):
pass
target_mu += partial[:, None] * (2.0*mu*self.variances)
target_S += partial[:, None] * self.variances
def psi1(self,Z,mu,S,target):
"""the variance, it does nothing"""
@ -149,7 +146,7 @@ class linear(kernpart):
def dpsi2_dZ(self,partial,Z,mu,S,target):
self._psi_computations(Z,mu,S)
mu2_S = np.sum(self.mu2_S,0)# Q,
target += (partial[:,:,:,None]* (Z * mu2_S * np.square(self.variances))).sum(0).sum(1)
target += (partial[:,:,:,None] * (self.mu2_S[:,None,None,:]*(Z*np.square(self.variances)[None,:])[None,None,:,:])).sum(0).sum(1)
#---------------------------------------#
# Precomputations #

7
GPy/kern/rbf.py
View file

@ -155,21 +155,20 @@ class rbf(kernpart):
self._psi_computations(Z,mu,S)
d_var = 2.*self._psi2/self.variance
d_length = self._psi2[:,:,:,None]*(0.5*self._psi2_Zdist_sq*self._psi2_denom + 2.*self._psi2_mudist_sq + 2.*S[:,None,None,:]/self.lengthscale2)/(self.lengthscale*self._psi2_denom)
d_length = d_length.sum(0)
target[0] += np.sum(partial*d_var)
dpsi2_dlength = d_length*partial[:,:,:,None]
if not self.ARD:
target[1] += dpsi2_dlength.sum()
else:
target[1:] += dpsi2_dlength.sum(0).sum(0).sum(0)
def dpsi2_dZ(self,partial,Z,mu,S,target):
self._psi_computations(Z,mu,S)
term1 = 0.5*self._psi2_Zdist/self.lengthscale2 # M, M, Q
term2 = self._psi2_mudist/self._psi2_denom/self.lengthscale2 # N, M, M, Q
dZ = self._psi2[:,:,:,None] * (term1[None] + term2)
target += (partial[:,:,:,None]*dZ).sum(0).sum(0) # <----------------- TODO not sure about the first ':' here, should be a None (WAS a none in the debug branch)
target += (partial[:,:,:,None]*dZ).sum(0).sum(0)
def dpsi2_dmuS(self,partial,Z,mu,S,target_mu,target_S):
"""Think N,M,M,Q """

2
GPy/models/sparse_GP.py
View file

@ -34,7 +34,7 @@ class sparse_GP(GP):
"""
def __init__(self, X, likelihood, kernel, Z, X_uncertainty=None, Xslices=None,Zslices=None, normalize_X=False):
self.scale_factor = 1.0# a scaling factor to help keep the algorithm stable
self.scale_factor = 100.0# a scaling factor to help keep the algorithm stable
self.Z = Z
self.Zslices = Zslices

48
GPy/testing/bgplvm_tests.py Normal file
View file

@ -0,0 +1,48 @@
# Copyright (c) 2012, Nicolo Fusi
# Licensed under the BSD 3-clause license (see LICENSE.txt)
import unittest
import numpy as np
import GPy
class BGPLVMTests(unittest.TestCase):
def test_bias_kern(self):
N, M, Q, D = 10, 3, 2, 4
X = np.random.rand(N, Q)
k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001)
K = k.K(X)
Y = np.random.multivariate_normal(np.zeros(N),K,D).T
k = GPy.kern.bias(Q) + GPy.kern.white(Q, 0.00001)
m = GPy.models.Bayesian_GPLVM(Y, Q, kernel = k, M=M)
m.constrain_positive('(rbf|bias|noise|white|S)')
m.randomize()
self.assertTrue(m.checkgrad())
def test_linear_kern(self):
N, M, Q, D = 10, 3, 2, 4
X = np.random.rand(N, Q)
k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001)
K = k.K(X)
Y = np.random.multivariate_normal(np.zeros(N),K,D).T
k = GPy.kern.linear(Q) + GPy.kern.white(Q, 0.00001)
m = GPy.models.Bayesian_GPLVM(Y, Q, kernel = k, M=M)
m.constrain_positive('(linear|bias|noise|white|S)')
m.randomize()
self.assertTrue(m.checkgrad())
def test_rbf_kern(self):
N, M, Q, D = 10, 3, 2, 4
X = np.random.rand(N, Q)
k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001)
K = k.K(X)
Y = np.random.multivariate_normal(np.zeros(N),K,D).T
k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001)
m = GPy.models.Bayesian_GPLVM(Y, Q, kernel = k, M=M)
m.constrain_positive('(rbf|bias|noise|white|S)')
m.randomize()
self.assertTrue(m.checkgrad())
if __name__ == "__main__":
print "Running unit tests, please be (very) patient..."
unittest.main()