mirror of
https://github.com/SheffieldML/GPy.git
synced 2026-06-08 15:05:15 +02:00
Merge branch 'master' of github.com:SheffieldML/GPy
This commit is contained in:
Ricardo Andrade
commit
5011afda06
12 changed files with 289 additions and 97 deletions
|
|
@ -121,9 +121,6 @@ class model(parameterised):
|
|||
else:
|
||||
raise AttributeError, "no parameter matches %s"%name
|
||||
|
||||
|
||||
|
||||
|
||||
def log_prior(self):
|
||||
"""evaluate the prior"""
|
||||
return np.sum([p.lnpdf(x) for p, x in zip(self.priors,self._get_params()) if p is not None])
|
||||
|
|
@ -135,12 +132,11 @@ class model(parameterised):
|
|||
[np.put(ret,i,p.lnpdf_grad(xx)) for i,(p,xx) in enumerate(zip(self.priors,x)) if not p is None]
|
||||
return ret
|
||||
|
||||
def _log_likelihood_gradients_transformed(self):
|
||||
def _transform_gradients(self, g):
|
||||
"""
|
||||
Use self.log_likelihood_gradients and self.prior_gradients to get the gradients of the model.
|
||||
Adjust the gradient for constraints and ties, return.
|
||||
Takes a list of gradients and return an array of transformed gradients (positive/negative/tied/and so on)
|
||||
"""
|
||||
g = self._log_likelihood_gradients() + self._log_prior_gradients()
|
||||
|
||||
x = self._get_params()
|
||||
g[self.constrained_positive_indices] = g[self.constrained_positive_indices]*x[self.constrained_positive_indices]
|
||||
g[self.constrained_negative_indices] = g[self.constrained_negative_indices]*x[self.constrained_negative_indices]
|
||||
|
|
@ -152,6 +148,7 @@ class model(parameterised):
|
|||
else:
|
||||
return g
|
||||
|
||||
|
||||
def randomize(self):
|
||||
"""
|
||||
Randomize the model.
|
||||
|
|
@ -241,6 +238,27 @@ class model(parameterised):
|
|||
print "Warning! constraining %s postive"%name
|
||||
|
||||
|
||||
def objective_function(self, x):
|
||||
"""
|
||||
The objective function passed to the optimizer. It combines the likelihood and the priors.
|
||||
"""
|
||||
self._set_params_transformed(x)
|
||||
return -self.log_likelihood() - self.log_prior()
|
||||
|
||||
def objective_function_gradients(self, x):
|
||||
"""
|
||||
Gets the gradients from the likelihood and the priors.
|
||||
"""
|
||||
self._set_params_transformed(x)
|
||||
LL_gradients = self._transform_gradients(self._log_likelihood_gradients())
|
||||
prior_gradients = self._transform_gradients(self._log_prior_gradients())
|
||||
return -LL_gradients - prior_gradients
|
||||
|
||||
def objective_and_gradients(self, x):
|
||||
obj_f = self.objective_function(x)
|
||||
obj_grads = self.objective_function_gradients(x)
|
||||
return obj_f, obj_grads
|
||||
|
||||
def optimize(self, optimizer=None, start=None, **kwargs):
|
||||
"""
|
||||
Optimize the model using self.log_likelihood and self.log_likelihood_gradient, as well as self.priors.
|
||||
|
|
@ -254,22 +272,12 @@ class model(parameterised):
|
|||
if optimizer is None:
|
||||
optimizer = self.preferred_optimizer
|
||||
|
||||
def f(x):
|
||||
self._set_params_transformed(x)
|
||||
return -self.log_likelihood()-self.log_prior()
|
||||
def fp(x):
|
||||
self._set_params_transformed(x)
|
||||
return -self._log_likelihood_gradients_transformed()
|
||||
def f_fp(x):
|
||||
self._set_params_transformed(x)
|
||||
return -self.log_likelihood()-self.log_prior(),-self._log_likelihood_gradients_transformed()
|
||||
|
||||
if start == None:
|
||||
start = self._get_params_transformed()
|
||||
|
||||
optimizer = optimization.get_optimizer(optimizer)
|
||||
opt = optimizer(start, model = self, **kwargs)
|
||||
opt.run(f_fp=f_fp, f=f, fp=fp)
|
||||
opt.run(f_fp=self.objective_and_gradients, f=self.objective_function, fp=self.objective_function_gradients)
|
||||
self.optimization_runs.append(opt)
|
||||
|
||||
self._set_params_transformed(opt.x_opt)
|
||||
|
|
@ -357,12 +365,9 @@ class model(parameterised):
|
|||
dx = step*np.sign(np.random.uniform(-1,1,x.size))
|
||||
|
||||
#evaulate around the point x
|
||||
self._set_params_transformed(x+dx)
|
||||
f1,g1 = self.log_likelihood() + self.log_prior(), self._log_likelihood_gradients_transformed()
|
||||
self._set_params_transformed(x-dx)
|
||||
f2,g2 = self.log_likelihood() + self.log_prior(), self._log_likelihood_gradients_transformed()
|
||||
self._set_params_transformed(x)
|
||||
gradient = self._log_likelihood_gradients_transformed()
|
||||
f1, g1 = self.objective_and_gradients(x+dx)
|
||||
f2, g2 = self.objective_and_gradients(x-dx)
|
||||
gradient = self.objective_function_gradients(x)
|
||||
|
||||
numerical_gradient = (f1-f2)/(2*dx)
|
||||
global_ratio = (f1-f2)/(2*np.dot(dx,gradient))
|
||||
|
|
@ -398,14 +403,10 @@ class model(parameterised):
|
|||
for i in param_list:
|
||||
xx = x.copy()
|
||||
xx[i] += step
|
||||
self._set_params_transformed(xx)
|
||||
f1,g1 = self.log_likelihood() + self.log_prior(), self._log_likelihood_gradients_transformed()[i]
|
||||
f1, g1 = self.objective_and_gradients(xx)
|
||||
xx[i] -= 2.*step
|
||||
self._set_params_transformed(xx)
|
||||
f2,g2 = self.log_likelihood() + self.log_prior(), self._log_likelihood_gradients_transformed()[i]
|
||||
self._set_params_transformed(x)
|
||||
gradient = self._log_likelihood_gradients_transformed()[i]
|
||||
|
||||
f2, g2 = self.objective_and_gradients(xx)
|
||||
gradient = self.objective_function_gradients(x)[i]
|
||||
|
||||
numerical_gradient = (f1-f2)/(2*step)
|
||||
ratio = (f1-f2)/(2*step*gradient)
|
||||
|
|
|
|||
|
|
@ -41,7 +41,7 @@ m.constrain_positive('(rbf|bias|S|linear|white|noise)')
|
|||
# m.unconstrain('white')
|
||||
# m.constrain_bounded('white', 1e-6, 10.0)
|
||||
# plot_oil(m.X, np.array([1,1]), labels, 'PCA initialization')
|
||||
m.optimize(messages = True)
|
||||
#m.optimize(messages = True)
|
||||
# m.optimize('tnc', messages = True)
|
||||
# plot_oil(m.X, m.kern.parts[0].lengthscale, labels, 'B-GPLVM')
|
||||
# # pb.figure()
|
||||
|
|
|
|||
56
GPy/examples/tuto_GP_regression.py
Normal file
56
GPy/examples/tuto_GP_regression.py
Normal file
|
|
@ -0,0 +1,56 @@
|
|||
# The detailed explanations of the commands used in this file can be found in the tutorial section
|
||||
|
||||
import pylab as pb
|
||||
pb.ion()
|
||||
import numpy as np
|
||||
import GPy
|
||||
|
||||
X = np.random.uniform(-3.,3.,(20,1))
|
||||
Y = np.sin(X) + np.random.randn(20,1)*0.05
|
||||
|
||||
kernel = GPy.kern.rbf(D=1, variance=1., lengthscale=1.)
|
||||
|
||||
m = GPy.models.GP_regression(X,Y,kernel)
|
||||
|
||||
print m
|
||||
m.plot()
|
||||
|
||||
m.constrain_positive('')
|
||||
|
||||
m.unconstrain('') # Required to remove the previous constrains
|
||||
m.constrain_positive('rbf_variance')
|
||||
m.constrain_bounded('lengthscale',1.,10. )
|
||||
m.constrain_fixed('noise',0.0025)
|
||||
|
||||
m.optimize()
|
||||
|
||||
m.optimize_restarts(Nrestarts = 10)
|
||||
|
||||
###########################
|
||||
# 2-dimensional example #
|
||||
###########################
|
||||
|
||||
import pylab as pb
|
||||
pb.ion()
|
||||
import numpy as np
|
||||
import GPy
|
||||
|
||||
# sample inputs and outputs
|
||||
X = np.random.uniform(-3.,3.,(50,2))
|
||||
Y = np.sin(X[:,0:1]) * np.sin(X[:,1:2])+np.random.randn(50,1)*0.05
|
||||
|
||||
# define kernel
|
||||
ker = GPy.kern.Matern52(2,ARD=True) + GPy.kern.white(2)
|
||||
|
||||
# create simple GP model
|
||||
m = GPy.models.GP_regression(X,Y,ker)
|
||||
|
||||
# contrain all parameters to be positive
|
||||
m.constrain_positive('')
|
||||
|
||||
# optimize and plot
|
||||
pb.figure()
|
||||
m.optimize('tnc', max_f_eval = 1000)
|
||||
|
||||
m.plot()
|
||||
print(m)
|
||||
139
GPy/examples/tuto_kernel_overview.py
Normal file
139
GPy/examples/tuto_kernel_overview.py
Normal file
|
|
@ -0,0 +1,139 @@
|
|||
# The detailed explanations of the commands used in this file can be found in the tutorial section
|
||||
|
||||
import pylab as pb
|
||||
import numpy as np
|
||||
import GPy
|
||||
pb.ion()
|
||||
|
||||
ker1 = GPy.kern.rbf(1) # Equivalent to ker1 = GPy.kern.rbf(D=1, variance=1., lengthscale=1.)
|
||||
ker2 = GPy.kern.rbf(D=1, variance = .75, lengthscale=2.)
|
||||
ker3 = GPy.kern.rbf(1, .5, .5)
|
||||
|
||||
print ker2
|
||||
ker1.plot()
|
||||
ker2.plot()
|
||||
ker3.plot()
|
||||
|
||||
k1 = GPy.kern.rbf(1,1.,2.)
|
||||
k2 = GPy.kern.Matern32(1, 0.5, 0.2)
|
||||
|
||||
# Product of kernels
|
||||
k_prod = k1.prod(k2)
|
||||
k_prodorth = k1.prod_orthogonal(k2)
|
||||
|
||||
# Sum of kernels
|
||||
k_add = k1.add(k2)
|
||||
k_addorth = k1.add_orthogonal(k2)
|
||||
|
||||
pb.figure(figsize=(8,8))
|
||||
pb.subplot(2,2,1)
|
||||
k_prod.plot()
|
||||
pb.title('prod')
|
||||
pb.subplot(2,2,2)
|
||||
k_prodorth.plot()
|
||||
pb.title('prod_orthogonal')
|
||||
pb.subplot(2,2,3)
|
||||
k_add.plot()
|
||||
pb.title('add')
|
||||
pb.subplot(2,2,4)
|
||||
k_addorth.plot()
|
||||
pb.title('add_orthogonal')
|
||||
pb.subplots_adjust(wspace=0.3, hspace=0.3)
|
||||
|
||||
k1 = GPy.kern.rbf(1,1.,2)
|
||||
k2 = GPy.kern.periodic_Matern52(1,variance=1e3, lengthscale=1, period = 1.5, lower=-5., upper = 5)
|
||||
|
||||
k = k1 * k2 # equivalent to k = k1.prod(k2)
|
||||
print k
|
||||
|
||||
# Simulate sample paths
|
||||
X = np.linspace(-5,5,501)[:,None]
|
||||
Y = np.random.multivariate_normal(np.zeros(501),k.K(X),1)
|
||||
|
||||
# plot
|
||||
pb.figure(figsize=(10,4))
|
||||
pb.subplot(1,2,1)
|
||||
k.plot()
|
||||
pb.subplot(1,2,2)
|
||||
pb.plot(X,Y.T)
|
||||
pb.ylabel("Sample path")
|
||||
pb.subplots_adjust(wspace=0.3)
|
||||
|
||||
k = (k1+k2)*(k1+k2)
|
||||
print k.parts[0].name, '\n', k.parts[1].name, '\n', k.parts[2].name, '\n', k.parts[3].name
|
||||
|
||||
k1 = GPy.kern.rbf(1)
|
||||
k2 = GPy.kern.Matern32(1)
|
||||
k3 = GPy.kern.white(1)
|
||||
|
||||
k = k1 + k2 + k3
|
||||
print k
|
||||
|
||||
k.constrain_positive('var')
|
||||
k.constrain_fixed(np.array([1]),1.75)
|
||||
k.tie_param('len')
|
||||
k.unconstrain('white')
|
||||
k.constrain_bounded('white',lower=1e-5,upper=.5)
|
||||
print k
|
||||
|
||||
k_cst = GPy.kern.bias(1,variance=1.)
|
||||
k_mat = GPy.kern.Matern52(1,variance=1., lengthscale=3)
|
||||
Kanova = (k_cst + k_mat).prod_orthogonal(k_cst + k_mat)
|
||||
print Kanova
|
||||
|
||||
# sample inputs and outputs
|
||||
X = np.random.uniform(-3.,3.,(40,2))
|
||||
Y = 0.5*X[:,:1] + 0.5*X[:,1:] + 2*np.sin(X[:,:1]) * np.sin(X[:,1:])
|
||||
|
||||
# Create GP regression model
|
||||
m = GPy.models.GP_regression(X,Y,Kanova)
|
||||
pb.figure(figsize=(5,5))
|
||||
m.plot()
|
||||
|
||||
pb.figure(figsize=(20,3))
|
||||
pb.subplots_adjust(wspace=0.5)
|
||||
pb.subplot(1,5,1)
|
||||
m.plot()
|
||||
pb.subplot(1,5,2)
|
||||
pb.ylabel("= ",rotation='horizontal',fontsize='30')
|
||||
pb.subplot(1,5,3)
|
||||
m.plot(which_functions=[False,True,False,False])
|
||||
pb.ylabel("cst +",rotation='horizontal',fontsize='30')
|
||||
pb.subplot(1,5,4)
|
||||
m.plot(which_functions=[False,False,True,False])
|
||||
pb.ylabel("+ ",rotation='horizontal',fontsize='30')
|
||||
pb.subplot(1,5,5)
|
||||
pb.ylabel("+ ",rotation='horizontal',fontsize='30')
|
||||
m.plot(which_functions=[False,False,False,True])
|
||||
|
||||
import pylab as pb
|
||||
import numpy as np
|
||||
import GPy
|
||||
pb.ion()
|
||||
|
||||
ker1 = GPy.kern.rbf(D=1) # Equivalent to ker1 = GPy.kern.rbf(D=1, variance=1., lengthscale=1.)
|
||||
ker2 = GPy.kern.rbf(D=1, variance = .75, lengthscale=3.)
|
||||
ker3 = GPy.kern.rbf(1, .5, .25)
|
||||
|
||||
ker1.plot()
|
||||
ker2.plot()
|
||||
ker3.plot()
|
||||
#pb.savefig("Figures/tuto_kern_overview_basicdef.png")
|
||||
|
||||
kernels = [GPy.kern.rbf(1), GPy.kern.exponential(1), GPy.kern.Matern32(1), GPy.kern.Matern52(1), GPy.kern.Brownian(1), GPy.kern.bias(1), GPy.kern.linear(1), GPy.kern.spline(1), GPy.kern.periodic_exponential(1), GPy.kern.periodic_Matern32(1), GPy.kern.periodic_Matern52(1), GPy.kern.white(1)]
|
||||
kernel_names = ["GPy.kern.rbf", "GPy.kern.exponential", "GPy.kern.Matern32", "GPy.kern.Matern52", "GPy.kern.Brownian", "GPy.kern.bias", "GPy.kern.linear", "GPy.kern.spline", "GPy.kern.periodic_exponential", "GPy.kern.periodic_Matern32", "GPy.kern.periodic_Matern52", "GPy.kern.white"]
|
||||
|
||||
pb.figure(figsize=(16,12))
|
||||
pb.subplots_adjust(wspace=.5, hspace=.5)
|
||||
for i, kern in enumerate(kernels):
|
||||
pb.subplot(3,4,i+1)
|
||||
kern.plot(x=7.5,plot_limits=[0.00001,15.])
|
||||
pb.title(kernel_names[i]+ '\n')
|
||||
|
||||
# actual plot for the noise
|
||||
i = 11
|
||||
X = np.linspace(0.,15.,201)
|
||||
WN = 0*X
|
||||
WN[100] = 1.
|
||||
pb.subplot(3,4,i+1)
|
||||
pb.plot(X,WN,'b')
|
||||
|
|
@ -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, coregionalise
|
||||
from constructors import rbf, Matern32, Matern52, exponential, linear, white, bias, finite_dimensional, spline, Brownian, rbf_sympy, sympykern, periodic_exponential, periodic_Matern32, periodic_Matern52, prod, prod_orthogonal, symmetric, coregionalise
|
||||
from kern import kern
|
||||
|
|
|
|||
|
|
@ -18,8 +18,8 @@ from Brownian import Brownian as Brownianpart
|
|||
from periodic_exponential import periodic_exponential as periodic_exponentialpart
|
||||
from periodic_Matern32 import periodic_Matern32 as periodic_Matern32part
|
||||
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 prod import prod as prodpart
|
||||
from prod_orthogonal import prod_orthogonal as prod_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
|
||||
|
|
@ -245,7 +245,7 @@ def periodic_Matern52(D,variance=1., lengthscale=None, period=2*np.pi,n_freq=10,
|
|||
part = periodic_Matern52part(D,variance, lengthscale, period, n_freq, lower, upper)
|
||||
return kern(D, [part])
|
||||
|
||||
def product(k1,k2):
|
||||
def prod(k1,k2):
|
||||
"""
|
||||
Construct a product kernel over D from two kernels over D
|
||||
|
||||
|
|
@ -253,10 +253,10 @@ def product(k1,k2):
|
|||
:type k1, k2: kernpart
|
||||
:rtype: kernel object
|
||||
"""
|
||||
part = productpart(k1,k2)
|
||||
part = prodpart(k1,k2)
|
||||
return kern(k1.D, [part])
|
||||
|
||||
def product_orthogonal(k1,k2):
|
||||
def prod_orthogonal(k1,k2):
|
||||
"""
|
||||
Construct a product kernel over D1 x D2 from a kernel over D1 and another over D2.
|
||||
|
||||
|
|
@ -264,7 +264,7 @@ def product_orthogonal(k1,k2):
|
|||
:type k1, k2: kernpart
|
||||
:rtype: kernel object
|
||||
"""
|
||||
part = product_orthogonalpart(k1,k2)
|
||||
part = prod_orthogonalpart(k1,k2)
|
||||
return kern(k1.D+k2.D, [part])
|
||||
|
||||
def symmetric(k):
|
||||
|
|
|
|||
|
|
@ -7,8 +7,8 @@ import pylab as pb
|
|||
from ..core.parameterised import parameterised
|
||||
from kernpart import kernpart
|
||||
import itertools
|
||||
from product_orthogonal import product_orthogonal
|
||||
from product import product
|
||||
from prod_orthogonal import prod_orthogonal
|
||||
from prod import prod
|
||||
|
||||
class kern(parameterised):
|
||||
def __init__(self,D,parts=[], input_slices=None):
|
||||
|
|
@ -161,7 +161,7 @@ class kern(parameterised):
|
|||
K1 = self.copy()
|
||||
K2 = other.copy()
|
||||
|
||||
newkernparts = [product(k1,k2) for k1, k2 in itertools.product(K1.parts,K2.parts)]
|
||||
newkernparts = [prod(k1,k2) for k1, k2 in itertools.product(K1.parts,K2.parts)]
|
||||
|
||||
slices = []
|
||||
for sl1, sl2 in itertools.product(K1.input_slices,K2.input_slices):
|
||||
|
|
@ -183,7 +183,7 @@ class kern(parameterised):
|
|||
K1 = self.copy()
|
||||
K2 = other.copy()
|
||||
|
||||
newkernparts = [product_orthogonal(k1,k2) for k1, k2 in itertools.product(K1.parts,K2.parts)]
|
||||
newkernparts = [prod_orthogonal(k1,k2) for k1, k2 in itertools.product(K1.parts,K2.parts)]
|
||||
|
||||
slices = []
|
||||
for sl1, sl2 in itertools.product(K1.input_slices,K2.input_slices):
|
||||
|
|
@ -371,16 +371,17 @@ class kern(parameterised):
|
|||
|
||||
def psi2(self,Z,mu,S,slices1=None,slices2=None):
|
||||
"""
|
||||
:Z: np.ndarray of inducing inputs (M x Q)
|
||||
: mu, S: np.ndarrays of means and variacnes (each N x Q)
|
||||
:returns psi2: np.ndarray (N,M,M,Q) """
|
||||
:param Z: np.ndarray of inducing inputs (M x Q)
|
||||
:param mu, S: np.ndarrays of means and variances (each N x Q)
|
||||
:returns psi2: np.ndarray (N,M,M)
|
||||
"""
|
||||
target = np.zeros((mu.shape[0],Z.shape[0],Z.shape[0]))
|
||||
slices1, slices2 = self._process_slices(slices1,slices2)
|
||||
[p.psi2(Z[s2,i_s],mu[s1,i_s],S[s1,i_s],target[s1,s2,s2]) for p,i_s,s1,s2 in zip(self.parts,self.input_slices,slices1,slices2)]
|
||||
|
||||
#compute the "cross" terms
|
||||
for p1, p2 in itertools.combinations(self.parts,2):
|
||||
#white doesn;t compine with anything
|
||||
#white doesn;t combine with anything
|
||||
if p1.name=='white' or p2.name=='white':
|
||||
pass
|
||||
#rbf X bias
|
||||
|
|
@ -396,28 +397,9 @@ class kern(parameterised):
|
|||
else:
|
||||
raise NotImplementedError, "psi2 cannot be computed for this kernel"
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
# "crossterms". Here we are recomputing psi1 for white (we don't need to), but it's
|
||||
# not really expensive, since it's just a matrix of zeroes.
|
||||
# psi1_matrices = [np.zeros((mu.shape[0], Z.shape[0])) for p in self.parts]
|
||||
# [p.psi1(Z[s2],mu[s1],S[s1],psi1_target[s1,s2]) for p,s1,s2,psi1_target in zip(self.parts,slices1,slices2, psi1_matrices)]
|
||||
|
||||
crossterms = 0.0
|
||||
# for 3 kernels this returns something like
|
||||
# [(0,1), (0,2), (1,2)]
|
||||
# in theory, we should also account for (1,0), (2,0) and so on, but
|
||||
# the transpose deals exactly with that
|
||||
# for a,b in itertools.combinations(psi1_matrices, 2):
|
||||
# tmp = np.multiply(a,b)
|
||||
# crossterms += tmp[:,None,:] + tmp[:, :,None]
|
||||
|
||||
return target + crossterms
|
||||
return target
|
||||
|
||||
def dpsi2_dtheta(self,partial,partial1,Z,mu,S,slices1=None,slices2=None):
|
||||
"""Returns shape (N,M,M,Ntheta)"""
|
||||
slices1, slices2 = self._process_slices(slices1,slices2)
|
||||
target = np.zeros(self.Nparam)
|
||||
[p.dpsi2_dtheta(partial[s1,s2,s2],Z[s2,i_s],mu[s1,i_s],S[s1,i_s],target[ps]) for p,i_s,s1,s2,ps in zip(self.parts,self.input_slices,slices1,slices2,self.param_slices)]
|
||||
|
|
@ -429,7 +411,7 @@ class kern(parameterised):
|
|||
ipsl1, ipsl2 = self.input_slices[i1], self.input_slices[i2]
|
||||
ps1, ps2 = self.param_slices[i1], self.param_slices[i2]
|
||||
|
||||
#white doesn;t compine with anything
|
||||
#white doesn;t combine with anything
|
||||
if p1.name=='white' or p2.name=='white':
|
||||
pass
|
||||
#rbf X bias
|
||||
|
|
@ -447,26 +429,6 @@ class kern(parameterised):
|
|||
else:
|
||||
raise NotImplementedError, "psi2 cannot be computed for this kernel"
|
||||
|
||||
# # "crossterms"
|
||||
# # 1. get all the psi1 statistics
|
||||
# psi1_matrices = [np.zeros((mu.shape[0], Z.shape[0])) for p in self.parts]
|
||||
# [p.psi1(Z[s2],mu[s1],S[s1],psi1_target[s1,s2]) for p,s1,s2,psi1_target in zip(self.parts,slices1,slices2, psi1_matrices)]
|
||||
|
||||
# partial1 = np.ones_like(partial1)
|
||||
# # 2. get all the dpsi1/dtheta gradients
|
||||
# psi1_gradients = [np.zeros(self.Nparam) for p in self.parts]
|
||||
# [p.dpsi1_dtheta(partial1[s2,s1],Z[s2,i_s],mu[s1,i_s],S[s1,i_s],psi1g_target[ps]) for p,ps,s1,s2,i_s,psi1g_target in zip(self.parts, self.param_slices,slices1,slices2,self.input_slices,psi1_gradients)]
|
||||
|
||||
|
||||
# # 3. multiply them somehow
|
||||
# for a,b in itertools.combinations(range(len(psi1_matrices)), 2):
|
||||
|
||||
# tmp = (psi1_gradients[a][None, None] * psi1_matrices[b][:,:, None])
|
||||
# # target += (tmp[None] + tmp[:,None]).sum(0).sum(0).sum(0)
|
||||
# # gne = (psi1_gradients[a].sum()*psi1_matrices[b].sum())
|
||||
# # target += gne
|
||||
# #target += (gne[None] + gne[:, None]).sum(0)
|
||||
# target += (partial.sum(0)[:,:,None] * (tmp[:, None] + tmp[:,:,None]).sum(0)).sum(0).sum(0)
|
||||
return self._transform_gradients(target)
|
||||
|
||||
def dpsi2_dZ(self,partial,Z,mu,S,slices1=None,slices2=None):
|
||||
|
|
@ -475,16 +437,15 @@ class kern(parameterised):
|
|||
[p.dpsi2_dZ(partial[s1,s2,s2],Z[s2,i_s],mu[s1,i_s],S[s1,i_s],target[s2,i_s]) for p,i_s,s1,s2 in zip(self.parts,self.input_slices,slices1,slices2)]
|
||||
|
||||
#compute the "cross" terms
|
||||
#TODO: slices (need to iterate around the input slices also...)
|
||||
for p1, p2 in itertools.combinations(self.parts,2):
|
||||
#white doesn;t compine with anything
|
||||
#white doesn;t combine with anything
|
||||
if p1.name=='white' or p2.name=='white':
|
||||
pass
|
||||
#rbf X bias
|
||||
elif p1.name=='bias' and p2.name=='rbf':
|
||||
target += p2.dpsi1_dX(partial.sum(1)*p1.variance,Z,mu,S)
|
||||
target += p2.dpsi1_dX(partial.sum(1)*p1.variance,Z,mu,S,target)
|
||||
elif p2.name=='bias' and p1.name=='rbf':
|
||||
target += p1.dpsi1_dZ(partial.sum(2)*p2.variance,Z,mu,S)
|
||||
target += p1.dpsi1_dZ(partial.sum(2)*p2.variance,Z,mu,S,target)
|
||||
#rbf X linear
|
||||
elif p1.name=='linear' and p2.name=='rbf':
|
||||
raise NotImplementedError #TODO
|
||||
|
|
@ -502,7 +463,24 @@ class kern(parameterised):
|
|||
target_mu, target_S = np.zeros((2,mu.shape[0],mu.shape[1]))
|
||||
[p.dpsi2_dmuS(partial[s1,s2,s2],Z[s2,i_s],mu[s1,i_s],S[s1,i_s],target_mu[s1,i_s],target_S[s1,i_s]) for p,i_s,s1,s2 in zip(self.parts,self.input_slices,slices1,slices2)]
|
||||
|
||||
#TODO: there are some extra terms to compute here!
|
||||
#compute the "cross" terms
|
||||
for p1, p2 in itertools.combinations(self.parts,2):
|
||||
#white doesn;t combine with anything
|
||||
if p1.name=='white' or p2.name=='white':
|
||||
pass
|
||||
#rbf X bias
|
||||
elif p1.name=='bias' and p2.name=='rbf':
|
||||
target += p2.dpsi1_dmuS(partial.sum(1)*p1.variance,Z,mu,S,target_mu,target_S)
|
||||
elif p2.name=='bias' and p1.name=='rbf':
|
||||
target += p1.dpsi1_dmuS(partial.sum(2)*p2.variance,Z,mu,S,target_mu,target_S)
|
||||
#rbf X linear
|
||||
elif p1.name=='linear' and p2.name=='rbf':
|
||||
raise NotImplementedError #TODO
|
||||
elif p2.name=='linear' and p1.name=='rbf':
|
||||
raise NotImplementedError #TODO
|
||||
else:
|
||||
raise NotImplementedError, "psi2 cannot be computed for this kernel"
|
||||
|
||||
return target_mu, target_S
|
||||
|
||||
def plot(self, x = None, plot_limits=None,which_functions='all',resolution=None,*args,**kwargs):
|
||||
|
|
|
|||
|
|
@ -6,7 +6,7 @@ import numpy as np
|
|||
import hashlib
|
||||
#from scipy import integrate # This may not be necessary (Nicolas, 20th Feb)
|
||||
|
||||
class product(kernpart):
|
||||
class prod(kernpart):
|
||||
"""
|
||||
Computes the product of 2 kernels that are defined on the same space
|
||||
|
||||
|
|
@ -6,7 +6,7 @@ import numpy as np
|
|||
import hashlib
|
||||
#from scipy import integrate # This may not be necessary (Nicolas, 20th Feb)
|
||||
|
||||
class product_orthogonal(kernpart):
|
||||
class prod_orthogonal(kernpart):
|
||||
"""
|
||||
Computes the product of 2 kernels
|
||||
|
||||
17
doc/kernel_implementation.rst
Normal file
17
doc/kernel_implementation.rst
Normal file
|
|
@ -0,0 +1,17 @@
|
|||
|
||||
***************************
|
||||
List of implemented kernels
|
||||
***************************
|
||||
|
||||
The :math:`\checkmark` symbol represents the functions that have been implemented for each kernel.
|
||||
|
||||
.. |tick|
|
||||
|
||||
.. |tick| image:: tick.png
|
||||
|
||||
|
||||
====== =========== === ======= =========== =============== ======= =========== ====== ====== =======
|
||||
NAME get/set K Kdiag dK_dtheta dKdiag_dtheta dK_dX dKdiag_dX psi0 psi1 psi2
|
||||
====== =========== === ======= =========== =============== ======= =========== ====== ====== =======
|
||||
rbf \\checkmark y
|
||||
====== =========== === ======= =========== =============== ======= =========== ====== ====== =======
|
||||
|
|
@ -2,7 +2,7 @@
|
|||
Gaussian process regression tutorial
|
||||
*************************************
|
||||
|
||||
We will see in this tutorial the basics for building a 1 dimensional and a 2 dimensional Gaussian process regression model, also known as a kriging model.
|
||||
We will see in this tutorial the basics for building a 1 dimensional and a 2 dimensional Gaussian process regression model, also known as a kriging model. The code shown in this tutorial can be found without the comments at GPy/examples/tuto_GP_regression.py.
|
||||
|
||||
We first import the libraries we will need: ::
|
||||
|
||||
|
|
|
|||
|
|
@ -2,6 +2,7 @@
|
|||
****************************
|
||||
tutorial : A kernel overview
|
||||
****************************
|
||||
The aim of this tutorial is to give a better understanding of the kernel objects in GPy and to list the ones that are already implemented. The code shown in this tutorial can be found without the comments at GPy/examples/tuto_kernel_overview.py.
|
||||
|
||||
First we import the libraries we will need ::
|
||||
|
||||
|
|
@ -38,7 +39,7 @@ return::
|
|||
Implemented kernels
|
||||
===================
|
||||
|
||||
Many kernels are already implemented in GPy. Here is a summary of most of them:
|
||||
Many kernels are already implemented in GPy. A comprehensive list can be found `here <kernel_implementation.html>`_ . The following figure gives a summary of most of them:
|
||||
|
||||
.. figure:: Figures/tuto_kern_overview_allkern.png
|
||||
:align: center
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue