mirror of
https://github.com/SheffieldML/GPy.git
synced 2026-05-15 06:52:39 +02:00
Merge branch 'master' of github.com:SheffieldML/GPy
This commit is contained in:
Alan Saul
commit
8371804d56
13 changed files with 169 additions and 126 deletions
|
|
@ -6,5 +6,6 @@ import kern
|
|||
import models
|
||||
import inference
|
||||
import util
|
||||
import examples
|
||||
#import examples TODO: discuss!
|
||||
from core import priors
|
||||
|
|
|
|||
|
|
@ -80,19 +80,22 @@ class model(parameterised):
|
|||
for w in which:
|
||||
self.priors[w] = what
|
||||
|
||||
def get(self,name):
|
||||
def get(self,name, return_names=False):
|
||||
"""
|
||||
get a model parameter by name
|
||||
Get a model parameter by name. The name is applied as a regular expression and all parameters that match that regular expression are returned.
|
||||
"""
|
||||
matches = self.grep_param_names(name)
|
||||
if len(matches):
|
||||
if return_names:
|
||||
return self._get_params()[matches], np.asarray(self._get_param_names())[matches].tolist()
|
||||
else:
|
||||
return self._get_params()[matches]
|
||||
else:
|
||||
raise AttributeError, "no parameter matches %s"%name
|
||||
|
||||
def set(self,name,val):
|
||||
"""
|
||||
Set a model parameter by name
|
||||
Set model parameter(s) by name. The name is provided as a regular expression. All parameters matching that regular expression are set to ghe given value.
|
||||
"""
|
||||
matches = self.grep_param_names(name)
|
||||
if len(matches):
|
||||
|
|
@ -102,6 +105,20 @@ class model(parameterised):
|
|||
else:
|
||||
raise AttributeError, "no parameter matches %s"%name
|
||||
|
||||
def get_gradient(self,name, return_names=False):
|
||||
"""
|
||||
Get model gradient(s) by name. The name is applied as a regular expression and all parameters that match that regular expression are returned.
|
||||
"""
|
||||
matches = self.grep_param_names(name)
|
||||
if len(matches):
|
||||
if return_names:
|
||||
return self._log_likelihood_gradients()[matches], np.asarray(self._get_param_names())[matches].tolist()
|
||||
else:
|
||||
return self._log_likelihood_gradients()[matches]
|
||||
else:
|
||||
raise AttributeError, "no parameter matches %s"%name
|
||||
|
||||
|
||||
|
||||
|
||||
def log_prior(self):
|
||||
|
|
|
|||
|
|
@ -17,10 +17,8 @@ def toy_rbf_1d():
|
|||
# create simple GP model
|
||||
m = GPy.models.GP_regression(data['X'],data['Y'])
|
||||
|
||||
# contrain all parameters to be positive
|
||||
m.constrain_positive('')
|
||||
|
||||
# optimize
|
||||
m.ensure_default_constraints()
|
||||
m.optimize()
|
||||
|
||||
# plot
|
||||
|
|
@ -35,10 +33,8 @@ def rogers_girolami_olympics():
|
|||
# create simple GP model
|
||||
m = GPy.models.GP_regression(data['X'],data['Y'])
|
||||
|
||||
# contrain all parameters to be positive
|
||||
m.constrain_positive('')
|
||||
|
||||
# optimize
|
||||
m.ensure_default_constraints()
|
||||
m.optimize()
|
||||
|
||||
# plot
|
||||
|
|
@ -57,10 +53,8 @@ def toy_rbf_1d_50():
|
|||
# create simple GP model
|
||||
m = GPy.models.GP_regression(data['X'],data['Y'])
|
||||
|
||||
# contrain all parameters to be positive
|
||||
m.constrain_positive('')
|
||||
|
||||
# optimize
|
||||
m.ensure_default_constraints()
|
||||
m.optimize()
|
||||
|
||||
# plot
|
||||
|
|
@ -75,10 +69,8 @@ def silhouette():
|
|||
# create simple GP model
|
||||
m = GPy.models.GP_regression(data['X'],data['Y'])
|
||||
|
||||
# contrain all parameters to be positive
|
||||
m.constrain_positive('')
|
||||
|
||||
# optimize
|
||||
m.ensure_default_constraints()
|
||||
m.optimize()
|
||||
|
||||
print(m)
|
||||
|
|
@ -118,20 +110,15 @@ def multiple_optima(gene_number=937,resolution=80, model_restarts=10, seed=10000
|
|||
kern = GPy.kern.rbf(1, variance=np.random.exponential(1.), lengthscale=np.random.exponential(50.)) + GPy.kern.white(1,variance=np.random.exponential(1.))
|
||||
|
||||
m = GPy.models.GP_regression(data['X'],data['Y'], kernel=kern)
|
||||
params = m._get_params()
|
||||
optim_point_x[0] = params[1]
|
||||
optim_point_y[0] = np.log10(params[0]) - np.log10(params[2]);
|
||||
|
||||
# contrain all parameters to be positive
|
||||
m.constrain_positive('')
|
||||
optim_point_x[0] = m.get('rbf_lengthscale')
|
||||
optim_point_y[0] = np.log10(m.get('rbf_variance')) - np.log10(m.get('white_variance'));
|
||||
|
||||
# optimize
|
||||
m.ensure_default_constraints()
|
||||
m.optimize(xtol=1e-6,ftol=1e-6)
|
||||
|
||||
params = m._get_params()
|
||||
optim_point_x[1] = params[1]
|
||||
optim_point_y[1] = np.log10(params[0]) - np.log10(params[2]);
|
||||
print(m)
|
||||
optim_point_x[1] = m.get('rbf_lengthscale')
|
||||
optim_point_y[1] = np.log10(m.get('rbf_variance')) - np.log10(m.get('white_variance'));
|
||||
|
||||
pb.arrow(optim_point_x[0], optim_point_y[0], optim_point_x[1]-optim_point_x[0], optim_point_y[1]-optim_point_y[0], label=str(i), head_length=1, head_width=0.5, fc='k', ec='k')
|
||||
models.append(m)
|
||||
|
|
|
|||
|
|
@ -10,11 +10,11 @@ print "sparse GPLVM with RBF kernel"
|
|||
|
||||
N = 100
|
||||
M = 4
|
||||
Q = 1
|
||||
Q = 2
|
||||
D = 2
|
||||
#generate GPLVM-like data
|
||||
X = np.random.rand(N, Q)
|
||||
k = GPy.kern.rbf(Q, 1.0, 2.0) + GPy.kern.white(Q, 0.00001)
|
||||
k = GPy.kern.rbf(Q,1.,2*np.ones((1,))) + GPy.kern.white(Q, 0.00001)
|
||||
K = k.K(X)
|
||||
Y = np.random.multivariate_normal(np.zeros(N),K,D).T
|
||||
|
||||
|
|
|
|||
|
|
@ -20,43 +20,52 @@ class Matern32(kernpart):
|
|||
:type D: int
|
||||
:param variance: the variance :math:`\sigma^2`
|
||||
:type variance: float
|
||||
:param lengthscale: the lengthscales :math:`\ell_i`
|
||||
:param lengthscale: the lengthscale :math:`\ell_i`
|
||||
:type lengthscale: np.ndarray of size (D,)
|
||||
:rtype: kernel object
|
||||
|
||||
"""
|
||||
|
||||
def __init__(self,D,variance=1.,lengthscales=None):
|
||||
def __init__(self,D,variance=1.,lengthscale=None,ARD=False):
|
||||
self.D = D
|
||||
if lengthscales is not None:
|
||||
assert lengthscales.shape==(self.D,)
|
||||
else:
|
||||
lengthscales = np.ones(self.D)
|
||||
self.Nparam = self.D + 1
|
||||
self.ARD = ARD
|
||||
if ARD == False:
|
||||
self.Nparam = 2
|
||||
self.name = 'Mat32'
|
||||
self._set_params(np.hstack((variance,lengthscales)))
|
||||
if lengthscale is not None:
|
||||
assert lengthscale.shape == (1,)
|
||||
else:
|
||||
lengthscale = np.ones(1)
|
||||
else:
|
||||
self.Nparam = self.D + 1
|
||||
self.name = 'Mat32_ARD'
|
||||
if lengthscale is not None:
|
||||
assert lengthscale.shape == (self.D,)
|
||||
else:
|
||||
lengthscale = np.ones(self.D)
|
||||
self._set_params(np.hstack((variance,lengthscale)))
|
||||
|
||||
def _get_params(self):
|
||||
"""return the value of the parameters."""
|
||||
return np.hstack((self.variance,self.lengthscales))
|
||||
return np.hstack((self.variance,self.lengthscale))
|
||||
|
||||
def _set_params(self,x):
|
||||
"""set the value of the parameters."""
|
||||
assert x.size==(self.D+1)
|
||||
assert x.size == self.Nparam
|
||||
self.variance = x[0]
|
||||
self.lengthscales = x[1:]
|
||||
self.lengthscale = x[1:]
|
||||
|
||||
def _get_param_names(self):
|
||||
"""return parameter names."""
|
||||
if self.D==1:
|
||||
if self.Nparam == 2:
|
||||
return ['variance','lengthscale']
|
||||
else:
|
||||
return ['variance']+['lengthscale_%i'%i for i in range(self.lengthscales.size)]
|
||||
return ['variance']+['lengthscale_%i'%i for i in range(self.lengthscale.size)]
|
||||
|
||||
def K(self,X,X2,target):
|
||||
"""Compute the covariance matrix between X and X2."""
|
||||
if X2 is None: X2 = X
|
||||
dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscales),-1))
|
||||
dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscale),-1))
|
||||
np.add(self.variance*(1+np.sqrt(3.)*dist)*np.exp(-np.sqrt(3.)*dist), target,target)
|
||||
|
||||
def Kdiag(self,X,target):
|
||||
|
|
@ -66,13 +75,20 @@ class Matern32(kernpart):
|
|||
def dK_dtheta(self,partial,X,X2,target):
|
||||
"""derivative of the covariance matrix with respect to the parameters."""
|
||||
if X2 is None: X2 = X
|
||||
dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscales),-1))
|
||||
dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscale),-1))
|
||||
dvar = (1+np.sqrt(3.)*dist)*np.exp(-np.sqrt(3.)*dist)
|
||||
invdist = 1./np.where(dist!=0.,dist,np.inf)
|
||||
dist2M = np.square(X[:,None,:]-X2[None,:,:])/self.lengthscales**3
|
||||
dl = (self.variance* 3 * dist * np.exp(-np.sqrt(3.)*dist))[:,:,np.newaxis] * dist2M*invdist[:,:,np.newaxis]
|
||||
dist2M = np.square(X[:,None,:]-X2[None,:,:])/self.lengthscale**3
|
||||
#dl = (self.variance* 3 * dist * np.exp(-np.sqrt(3.)*dist))[:,:,np.newaxis] * dist2M*invdist[:,:,np.newaxis]
|
||||
target[0] += np.sum(dvar*partial)
|
||||
if self.ARD == True:
|
||||
dl = (self.variance* 3 * dist * np.exp(-np.sqrt(3.)*dist))[:,:,np.newaxis] * dist2M*invdist[:,:,np.newaxis]
|
||||
#dl = self.variance*dvar[:,:,None]*dist2M*invdist[:,:,None]
|
||||
target[1:] += (dl*partial[:,:,None]).sum(0).sum(0)
|
||||
else:
|
||||
dl = (self.variance* 3 * dist * np.exp(-np.sqrt(3.)*dist)) * dist2M.sum(-1)*invdist
|
||||
#dl = self.variance*dvar*dist2M.sum(-1)*invdist
|
||||
target[1] += np.sum(dl*partial)
|
||||
|
||||
def dKdiag_dtheta(self,partial,X,target):
|
||||
"""derivative of the diagonal of the covariance matrix with respect to the parameters."""
|
||||
|
|
@ -81,8 +97,8 @@ class Matern32(kernpart):
|
|||
def dK_dX(self,partial,X,X2,target):
|
||||
"""derivative of the covariance matrix with respect to X."""
|
||||
if X2 is None: X2 = X
|
||||
dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscales),-1))[:,:,None]
|
||||
ddist_dX = (X[:,None,:]-X2[None,:,:])/self.lengthscales**2/np.where(dist!=0.,dist,np.inf)
|
||||
dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscale),-1))[:,:,None]
|
||||
ddist_dX = (X[:,None,:]-X2[None,:,:])/self.lengthscale**2/np.where(dist!=0.,dist,np.inf)
|
||||
dK_dX = - np.transpose(3*self.variance*dist*np.exp(-np.sqrt(3)*dist)*ddist_dX,(1,0,2))
|
||||
target += np.sum(dK_dX*partial.T[:,:,None],0)
|
||||
|
||||
|
|
@ -104,7 +120,7 @@ class Matern32(kernpart):
|
|||
"""
|
||||
assert self.D == 1
|
||||
def L(x,i):
|
||||
return(3./self.lengthscales**2*F[i](x) + 2*np.sqrt(3)/self.lengthscales*F1[i](x) + F2[i](x))
|
||||
return(3./self.lengthscale**2*F[i](x) + 2*np.sqrt(3)/self.lengthscale*F1[i](x) + F2[i](x))
|
||||
n = F.shape[0]
|
||||
G = np.zeros((n,n))
|
||||
for i in range(n):
|
||||
|
|
@ -114,5 +130,5 @@ class Matern32(kernpart):
|
|||
F1lower = np.array([f(lower) for f in F1])[:,None]
|
||||
#print "OLD \n", np.dot(F1lower,F1lower.T), "\n \n"
|
||||
#return(G)
|
||||
return(self.lengthscales**3/(12.*np.sqrt(3)*self.variance) * G + 1./self.variance*np.dot(Flower,Flower.T) + self.lengthscales**2/(3.*self.variance)*np.dot(F1lower,F1lower.T))
|
||||
return(self.lengthscale**3/(12.*np.sqrt(3)*self.variance) * G + 1./self.variance*np.dot(Flower,Flower.T) + self.lengthscale**2/(3.*self.variance)*np.dot(F1lower,F1lower.T))
|
||||
|
||||
|
|
|
|||
|
|
@ -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, rbf_ARD, spline, Brownian, linear_ARD, rbf_sympy, sympykern
|
||||
from constructors import rbf, Matern32, Matern52, exponential, linear, white, bias, finite_dimensional, spline, Brownian, linear_ARD, rbf_sympy, sympykern
|
||||
from kern import kern
|
||||
|
|
|
|||
|
|
@ -22,7 +22,7 @@ from Brownian import Brownian as Brownianpart
|
|||
#using meta-classes to make the objects construct properly wthout them.
|
||||
|
||||
|
||||
def rbf(D,variance=1., lengthscale=1.):
|
||||
def rbf(D,variance=1., lengthscale=None,ARD=False):
|
||||
"""
|
||||
Construct an RBF kernel
|
||||
|
||||
|
|
@ -33,21 +33,7 @@ def rbf(D,variance=1., lengthscale=1.):
|
|||
:param lengthscale: the lengthscale of the kernel
|
||||
:type lengthscale: float
|
||||
"""
|
||||
part = rbfpart(D,variance,lengthscale)
|
||||
return kern(D, [part])
|
||||
|
||||
def rbf_ARD(D,variance=1., lengthscales=None):
|
||||
"""
|
||||
Construct an RBF kernel with Automatic Relevance Determination (ARD)
|
||||
|
||||
:param D: dimensionality of the kernel, obligatory
|
||||
:type D: int
|
||||
:param variance: the variance of the kernel
|
||||
:type variance: float
|
||||
:param lengthscales: the lengthscales of the kernel
|
||||
:type lengthscales: None|np.ndarray
|
||||
"""
|
||||
part = rbf_ARD_part(D,variance,lengthscales)
|
||||
part = rbfpart(D,variance,lengthscale,ARD)
|
||||
return kern(D, [part])
|
||||
|
||||
def linear(D,lengthscales=None):
|
||||
|
|
@ -86,7 +72,7 @@ def white(D,variance=1.):
|
|||
part = whitepart(D,variance)
|
||||
return kern(D, [part])
|
||||
|
||||
def exponential(D,variance=1., lengthscales=None):
|
||||
def exponential(D,variance=1., lengthscale=None, ARD=False):
|
||||
"""
|
||||
Construct a exponential kernel.
|
||||
|
||||
|
|
@ -96,10 +82,10 @@ def exponential(D,variance=1., lengthscales=None):
|
|||
variance (float)
|
||||
lengthscales (np.ndarray)
|
||||
"""
|
||||
part = exponentialpart(D,variance, lengthscales)
|
||||
part = exponentialpart(D,variance, lengthscale, ARD)
|
||||
return kern(D, [part])
|
||||
|
||||
def Matern32(D,variance=1., lengthscales=None):
|
||||
def Matern32(D,variance=1., lengthscale=None, ARD=False):
|
||||
"""
|
||||
Construct a Matern 3/2 kernel.
|
||||
|
||||
|
|
@ -109,7 +95,7 @@ def Matern32(D,variance=1., lengthscales=None):
|
|||
variance (float)
|
||||
lengthscales (np.ndarray)
|
||||
"""
|
||||
part = Matern32part(D,variance, lengthscales)
|
||||
part = Matern32part(D,variance, lengthscale, ARD)
|
||||
return kern(D, [part])
|
||||
|
||||
def Matern52(D,variance=1., lengthscales=None):
|
||||
|
|
|
|||
|
|
@ -24,37 +24,46 @@ class exponential(kernpart):
|
|||
:rtype: kernel object
|
||||
|
||||
"""
|
||||
def __init__(self,D,variance=1.,lengthscales=None):
|
||||
def __init__(self,D,variance=1.,lengthscale=None,ARD=False):
|
||||
self.D = D
|
||||
if lengthscales is not None:
|
||||
assert lengthscales.shape==(self.D,)
|
||||
else:
|
||||
lengthscales = np.ones(self.D)
|
||||
self.Nparam = self.D + 1
|
||||
self.ARD = ARD
|
||||
if ARD == False:
|
||||
self.Nparam = 2
|
||||
self.name = 'exp'
|
||||
self._set_params(np.hstack((variance,lengthscales)))
|
||||
if lengthscale is not None:
|
||||
assert lengthscale.shape == (1,)
|
||||
else:
|
||||
lengthscale = np.ones(1)
|
||||
else:
|
||||
self.Nparam = self.D + 1
|
||||
self.name = 'exp_ARD'
|
||||
if lengthscale is not None:
|
||||
assert lengthscale.shape == (self.D,)
|
||||
else:
|
||||
lengthscale = np.ones(self.D)
|
||||
self._set_params(np.hstack((variance,lengthscale)))
|
||||
|
||||
def _get_params(self):
|
||||
"""return the value of the parameters."""
|
||||
return np.hstack((self.variance,self.lengthscales))
|
||||
return np.hstack((self.variance,self.lengthscale))
|
||||
|
||||
def _set_params(self,x):
|
||||
"""set the value of the parameters."""
|
||||
assert x.size==(self.D+1)
|
||||
assert x.size == self.Nparam
|
||||
self.variance = x[0]
|
||||
self.lengthscales = x[1:]
|
||||
self.lengthscale = x[1:]
|
||||
|
||||
def _get_param_names(self):
|
||||
"""return parameter names."""
|
||||
if self.D==1:
|
||||
if self.Nparam == 2:
|
||||
return ['variance','lengthscale']
|
||||
else:
|
||||
return ['variance']+['lengthscale_%i'%i for i in range(self.lengthscales.size)]
|
||||
return ['variance']+['lengthscale_%i'%i for i in range(self.lengthscale.size)]
|
||||
|
||||
def K(self,X,X2,target):
|
||||
"""Compute the covariance matrix between X and X2."""
|
||||
if X2 is None: X2 = X
|
||||
dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscales),-1))
|
||||
dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscale),-1))
|
||||
np.add(self.variance*np.exp(-dist), target,target)
|
||||
|
||||
def Kdiag(self,X,target):
|
||||
|
|
@ -64,13 +73,17 @@ class exponential(kernpart):
|
|||
def dK_dtheta(self,partial,X,X2,target):
|
||||
"""derivative of the covariance matrix with respect to the parameters."""
|
||||
if X2 is None: X2 = X
|
||||
dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscales),-1))
|
||||
dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscale),-1))
|
||||
invdist = 1./np.where(dist!=0.,dist,np.inf)
|
||||
dist2M = np.square(X[:,None,:]-X2[None,:,:])/self.lengthscales**3
|
||||
dist2M = np.square(X[:,None,:]-X2[None,:,:])/self.lengthscale**3
|
||||
dvar = np.exp(-dist)
|
||||
dl = self.variance*dvar[:,:,None]*dist2M*invdist[:,:,None]
|
||||
target[0] += np.sum(dvar*partial)
|
||||
if self.ARD == True:
|
||||
dl = self.variance*dvar[:,:,None]*dist2M*invdist[:,:,None]
|
||||
target[1:] += (dl*partial[:,:,None]).sum(0).sum(0)
|
||||
else:
|
||||
dl = self.variance*dvar*dist2M.sum(-1)*invdist
|
||||
target[1] += np.sum(dl*partial)
|
||||
|
||||
def dKdiag_dtheta(self,partial,X,target):
|
||||
"""derivative of the diagonal of the covariance matrix with respect to the parameters."""
|
||||
|
|
@ -80,8 +93,8 @@ class exponential(kernpart):
|
|||
def dK_dX(self,partial,X,X2,target):
|
||||
"""derivative of the covariance matrix with respect to X."""
|
||||
if X2 is None: X2 = X
|
||||
dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscales),-1))[:,:,None]
|
||||
ddist_dX = (X[:,None,:]-X2[None,:,:])/self.lengthscales**2/np.where(dist!=0.,dist,np.inf)
|
||||
dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscale),-1))[:,:,None]
|
||||
ddist_dX = (X[:,None,:]-X2[None,:,:])/self.lengthscale**2/np.where(dist!=0.,dist,np.inf)
|
||||
dK_dX = - np.transpose(self.variance*np.exp(-dist)*ddist_dX,(1,0,2))
|
||||
target += np.sum(dK_dX*partial.T[:,:,None],0)
|
||||
|
||||
|
|
@ -101,14 +114,14 @@ class exponential(kernpart):
|
|||
"""
|
||||
assert self.D == 1
|
||||
def L(x,i):
|
||||
return(1./self.lengthscales*F[i](x) + F1[i](x))
|
||||
return(1./self.lengthscale*F[i](x) + F1[i](x))
|
||||
n = F.shape[0]
|
||||
G = np.zeros((n,n))
|
||||
for i in range(n):
|
||||
for j in range(i,n):
|
||||
G[i,j] = G[j,i] = integrate.quad(lambda x : L(x,i)*L(x,j),lower,upper)[0]
|
||||
Flower = np.array([f(lower) for f in F])[:,None]
|
||||
return(self.lengthscales/2./self.variance * G + 1./self.variance * np.dot(Flower,Flower.T))
|
||||
return(self.lengthscale/2./self.variance * G + 1./self.variance * np.dot(Flower,Flower.T))
|
||||
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -20,16 +20,32 @@ class rbf(kernpart):
|
|||
:type D: int
|
||||
:param variance: the variance of the kernel
|
||||
:type variance: float
|
||||
:param lengthscale: the lengthscale of the kernel
|
||||
:type lengthscale: float
|
||||
:param lengthscale: the vector of lengthscale of the kernel
|
||||
:type lengthscale: np.ndarray
|
||||
:param ARD: Auto Relevance Determination. If equal to "False", the kernel is isotropic (ie. one single lengthscale parameter \ell), otherwise there is one lengthscale parameter per dimension.
|
||||
:type ARD: Boolean
|
||||
|
||||
.. Note: for rbf with different lengthscale on each dimension, see rbf_ARD
|
||||
"""
|
||||
|
||||
def __init__(self,D,variance=1.,lengthscale=1.):
|
||||
def __init__(self,D,variance=1.,lengthscale=None,ARD=False):
|
||||
self.D = D
|
||||
self.ARD = ARD
|
||||
if ARD == False:
|
||||
self.Nparam = 2
|
||||
self.name = 'rbf'
|
||||
if lengthscale is not None:
|
||||
assert lengthscale.shape == (1,)
|
||||
else:
|
||||
lengthscale = np.ones(1)
|
||||
|
||||
else:
|
||||
self.Nparam = self.D + 1
|
||||
self.name = 'rbf_ARD'
|
||||
if lengthscale is not None:
|
||||
assert lengthscale.shape == (self.D,)
|
||||
else:
|
||||
lengthscale = np.ones(self.D)
|
||||
|
||||
self._set_params(np.hstack((variance,lengthscale)))
|
||||
|
||||
#initialize cache
|
||||
|
|
@ -40,14 +56,19 @@ class rbf(kernpart):
|
|||
return np.hstack((self.variance,self.lengthscale))
|
||||
|
||||
def _set_params(self,x):
|
||||
self.variance, self.lengthscale = x
|
||||
assert x.size==(self.Nparam)
|
||||
self.variance = x[0]
|
||||
self.lengthscale = x[1:]
|
||||
self.lengthscale2 = np.square(self.lengthscale)
|
||||
#reset cached results
|
||||
self._X, self._X2, self._params = np.empty(shape=(3,1))
|
||||
self._Z, self._mu, self._S = np.empty(shape=(3,1)) # cached versions of Z,mu,S
|
||||
|
||||
def _get_param_names(self):
|
||||
if self.Nparam == 2:
|
||||
return ['variance','lengthscale']
|
||||
else:
|
||||
return ['variance']+['lengthscale_%i'%i for i in range(self.lengthscale.size)]
|
||||
|
||||
def K(self,X,X2,target):
|
||||
if X2 is None:
|
||||
|
|
@ -61,7 +82,12 @@ class rbf(kernpart):
|
|||
def dK_dtheta(self,partial,X,X2,target):
|
||||
self._K_computations(X,X2)
|
||||
target[0] += np.sum(self._K_dvar*partial)
|
||||
target[1] += np.sum(self._K_dvar*self.variance*self._K_dist2/self.lengthscale*partial)
|
||||
if self.ARD == True:
|
||||
dl = self._K_dvar[:,:,None]*self.variance*self._K_dist2/self.lengthscale
|
||||
target[1:] += (dl*partial[:,:,None]).sum(0).sum(0)
|
||||
else:
|
||||
target[1] += np.sum(self._K_dvar*self.variance*(self._K_dist2.sum(-1))/self.lengthscale*partial)
|
||||
#np.sum(self._K_dvar*self.variance*self._K_dist2/self.lengthscale*partial)
|
||||
|
||||
def dKdiag_dtheta(self,partial,X,target):
|
||||
#NB: derivative of diagonal elements wrt lengthscale is 0
|
||||
|
|
@ -81,15 +107,13 @@ class rbf(kernpart):
|
|||
self._X = X
|
||||
self._X2 = X2
|
||||
if X2 is None: X2 = X
|
||||
XXT = np.dot(X,X2.T)
|
||||
if X is X2:
|
||||
self._K_dist2 = (-2.*XXT + np.diag(XXT)[:,np.newaxis] + np.diag(XXT)[np.newaxis,:])/self.lengthscale2
|
||||
else:
|
||||
self._K_dist2 = (-2.*XXT + np.sum(np.square(X),1)[:,None] + np.sum(np.square(X2),1)[None,:])/self.lengthscale2
|
||||
# TODO Remove comments if this is fine.
|
||||
# Commented out by Neil as doesn't seem to be used elsewhere.
|
||||
#self._K_exponent = -0.5*self._K_dist2
|
||||
self._K_dvar = np.exp(-0.5*self._K_dist2)
|
||||
self._K_dist = X[:,None,:]-X2[None,:,:] # this can be computationally heavy
|
||||
self._params = np.empty(shape=(1,0))#ensure the next section gets called
|
||||
if not np.all(self._params == self._get_params()):
|
||||
self._params == self._get_params()
|
||||
self._K_dist2 = np.square(self._K_dist/self.lengthscale)
|
||||
#self._K_exponent = -0.5*self._K_dist2.sum(-1) #ND: commented out because seems not to be used
|
||||
self._K_dvar = np.exp(-0.5*self._K_dist2.sum(-1))
|
||||
|
||||
def psi0(self,Z,mu,S,target):
|
||||
target += self.variance
|
||||
|
|
@ -132,7 +156,7 @@ class rbf(kernpart):
|
|||
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)
|
||||
target[1] += np.sum(d_length*partial)
|
||||
target[1:] += (d_length*partial[:,:,None]).sum(0).sum(0)
|
||||
|
||||
def dpsi2_dZ(self,partial,Z,mu,S,target):
|
||||
"""Returns shape N,M,M,Q"""
|
||||
|
|
@ -175,4 +199,3 @@ class rbf(kernpart):
|
|||
self._psi2 = np.square(self.variance)*np.exp(self._psi2_exponent) # N,M,M
|
||||
|
||||
self._Z, self._mu, self._S = Z, mu,S
|
||||
|
||||
|
|
|
|||
|
|
@ -63,10 +63,10 @@ class GP_regression(model):
|
|||
self._Ystd = np.ones((1,self.Y.shape[1]))
|
||||
|
||||
if self.D > self.N:
|
||||
# then it's more efficient to store Youter
|
||||
self.Youter = np.dot(self.Y, self.Y.T)
|
||||
# then it's more efficient to store YYT
|
||||
self.YYT = np.dot(self.Y, self.Y.T)
|
||||
else:
|
||||
self.Youter = None
|
||||
self.YYT = None
|
||||
|
||||
model.__init__(self)
|
||||
|
||||
|
|
@ -83,23 +83,23 @@ class GP_regression(model):
|
|||
|
||||
def _model_fit_term(self):
|
||||
"""
|
||||
Computes the model fit using Youter if it's available
|
||||
Computes the model fit using YYT if it's available
|
||||
"""
|
||||
if self.Youter is None:
|
||||
if self.YYT is None:
|
||||
return -0.5*np.sum(np.square(np.dot(self.Li,self.Y)))
|
||||
else:
|
||||
return -0.5*np.sum(np.multiply(self.Ki, self.Youter))
|
||||
return -0.5*np.sum(np.multiply(self.Ki, self.YYT))
|
||||
|
||||
def log_likelihood(self):
|
||||
complexity_term = -0.5*self.N*self.D*np.log(2.*np.pi) - 0.5*self.D*self.K_logdet
|
||||
return complexity_term + self._model_fit_term()
|
||||
|
||||
def dL_dK(self):
|
||||
if self.Youter is None:
|
||||
if self.YYT is None:
|
||||
alpha = np.dot(self.Ki,self.Y)
|
||||
dL_dK = 0.5*(np.dot(alpha,alpha.T)-self.D*self.Ki)
|
||||
else:
|
||||
dL_dK = 0.5*(mdot(self.Ki, self.Youter, self.Ki) - self.D*self.Ki)
|
||||
dL_dK = 0.5*(mdot(self.Ki, self.YYT, self.Ki) - self.D*self.Ki)
|
||||
|
||||
return dL_dK
|
||||
|
||||
|
|
|
|||
|
|
@ -91,9 +91,9 @@ class generalized_FITC(model):
|
|||
|
||||
def log_likelihood(self):
|
||||
self.posterior_param()
|
||||
self.Youter = np.dot(self.mu_tilde,self.mu_tilde.T)
|
||||
self.YYT = np.dot(self.mu_tilde,self.mu_tilde.T)
|
||||
A = -self.hld
|
||||
B = -.5*np.sum(self.Qi*self.Youter)
|
||||
B = -.5*np.sum(self.Qi*self.YYT)
|
||||
C = sum(np.log(self.ep_approx.Z_hat))
|
||||
D = .5*np.sum(np.log(1./self.ep_approx.tau_tilde + 1./self.ep_approx.tau_))
|
||||
E = .5*np.sum((self.ep_approx.v_/self.ep_approx.tau_ - self.mu_tilde.flatten())**2/(1./self.ep_approx.tau_ + 1./self.ep_approx.tau_tilde))
|
||||
|
|
|
|||
|
|
@ -48,9 +48,9 @@ class warpedGP(GP_regression):
|
|||
|
||||
# this supports the 'smart' behaviour in GP_regression
|
||||
if self.D > self.N:
|
||||
self.Youter = np.dot(self.Y, self.Y.T)
|
||||
self.YYT = np.dot(self.Y, self.Y.T)
|
||||
else:
|
||||
self.Youter = None
|
||||
self.YYT = None
|
||||
|
||||
return self.Y
|
||||
|
||||
|
|
|
|||
|
|
@ -121,7 +121,7 @@ class GradientTests(unittest.TestCase):
|
|||
""" Testing GPLVM with rbf + bias and white kernel """
|
||||
N, Q, D = 50, 1, 2
|
||||
X = np.random.rand(N, Q)
|
||||
k = GPy.kern.rbf(Q, 0.5, 0.9) + GPy.kern.bias(Q, 0.1) + GPy.kern.white(Q, 0.05)
|
||||
k = GPy.kern.rbf(Q, 0.5, 0.9*np.ones((1,))) + GPy.kern.bias(Q, 0.1) + GPy.kern.white(Q, 0.05)
|
||||
K = k.K(X)
|
||||
Y = np.random.multivariate_normal(np.zeros(N),K,D).T
|
||||
m = GPy.models.GPLVM(Y, Q, kernel = k)
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue