apunkt/GPy
1
0
Fork 0
mirror of https://github.com/SheffieldML/GPy.git synced 2026-06-08 15:05:15 +02:00
Code Issues Projects Releases Packages Wiki Activity Actions

kern params adapted: Nparams > num_params and fixes of input_dim

Browse source
This commit is contained in:
Max Zwiessele 2013-06-05 16:14:30 +01:00
parent aa7fd122ca
commit 0490861099
42 changed files with 480 additions and 502 deletions
Download patch file Download diff file Expand all files Collapse all files

12
GPy/core/fitc.py
View file

@ -57,8 +57,8 @@ class FITC(SparseGP):
def _set_params(self, p):
self.Z = p[:self.M * self.input_dim].reshape(self.M, self.input_dim)
self.kern._set_params(p[self.Z.size:self.Z.size + self.kern.Nparam])
self.likelihood._set_params(p[self.Z.size + self.kern.Nparam:])
self.kern._set_params(p[self.Z.size:self.Z.size + self.kern.num_params])
self.likelihood._set_params(p[self.Z.size + self.kern.num_params:])
self._compute_kernel_matrices()
self._computations()
@ -198,14 +198,14 @@ class FITC(SparseGP):
aux_1 = self.likelihood.Y.T * np.dot(self._LBi_Lmi_psi1V.T,LBiLmipsi1)
aux_2 = np.dot(LBiLmipsi1.T,self._LBi_Lmi_psi1V)
dA_dnoise = 0.5 * self.D * (dbstar_dnoise/self.beta_star).sum() - 0.5 * self.D * np.sum(self.likelihood.Y**2 * dbstar_dnoise)
dA_dnoise = 0.5 * self.input_dim * (dbstar_dnoise/self.beta_star).sum() - 0.5 * self.input_dim * np.sum(self.likelihood.Y**2 * dbstar_dnoise)
dC_dnoise = -0.5 * np.sum(mdot(self.LBi.T,self.LBi,Lmi_psi1) * Lmi_psi1 * dbstar_dnoise.T)
dC_dnoise = -0.5 * np.sum(mdot(self.LBi.T,self.LBi,Lmi_psi1) * Lmi_psi1 * dbstar_dnoise.T)
dD_dnoise_1 = mdot(self.V_star*LBiLmipsi1.T,LBiLmipsi1*dbstar_dnoise.T*self.likelihood.Y.T)
alpha = mdot(LBiLmipsi1,self.V_star)
alpha_ = mdot(LBiLmipsi1.T,alpha)
dD_dnoise_2 = -0.5 * self.D * np.sum(alpha_**2 * dbstar_dnoise )
dD_dnoise_2 = -0.5 * self.input_dim * np.sum(alpha_**2 * dbstar_dnoise )
dD_dnoise_1 = mdot(self.V_star.T,self.psi1.T,self.Lmi.T,self.LBi.T,self.LBi,self.Lmi,self.psi1,dbstar_dnoise*self.likelihood.Y)
dD_dnoise_2 = 0.5*mdot(self.V_star.T,self.psi1.T,Hi,self.psi1,dbstar_dnoise*self.psi1.T,Hi,self.psi1,self.V_star)
@ -270,8 +270,8 @@ class FITC(SparseGP):
# q(u|f) = N(u| R0i*mu_u*f, R0i*C*R0i.T)
# Ci = I + (RPT0)Di(RPT0).T
# C = I - [RPT0] * (D+[RPT0].T*[RPT0])^-1*[RPT0].T
# = I - [RPT0] * (D + self.Qnn)^-1 * [RPT0].T
# C = I - [RPT0] * (input_dim+[RPT0].T*[RPT0])^-1*[RPT0].T
# = I - [RPT0] * (input_dim + self.Qnn)^-1 * [RPT0].T
# = I - [RPT0] * (U*U.T)^-1 * [RPT0].T
# = I - V.T * V
U = np.linalg.cholesky(np.diag(self.Diag0) + self.Qnn)

24
GPy/core/gp_base.py
View file

@ -1,12 +1,12 @@
import numpy as np
import model
from .. import kern
from ..util.plot import gpplot, Tango, x_frame1D, x_frame2D
import pylab as pb
from GPy.core.model import Model
class GPBase(model.model):
class GPBase(Model):
"""
Gaussian Process model for holding shared behaviour between
Gaussian Process Model for holding shared behaviour between
sprase_GP and GP models
"""
@ -28,7 +28,7 @@ class GPBase(model.model):
self._Xmean = np.zeros((1, self.input_dim))
self._Xstd = np.ones((1, self.input_dim))
model.model.__init__(self)
Model.__init__(self)
# All leaf nodes should call self._set_params(self._get_params()) at
# the end
@ -84,8 +84,8 @@ class GPBase(model.model):
Xnew, xmin, xmax, xx, yy = x_frame2D(self.X, plot_limits, resolution)
m, v = self._raw_predict(Xnew, which_parts=which_parts)
m = m.reshape(resolution, resolution).T
ax.contour(xx, yy, m, vmin=m.min(), vmax=m.max(), cmap=pb.cm.jet)
ax.scatter(self.X[:, 0], self.X[:, 1], 40, self.likelihood.Y, linewidth=0, cmap=pb.cm.jet, vmin=m.min(), vmax=m.max())
ax.contour(xx, yy, m, vmin=m.min(), vmax=m.max(), cmap=pb.cm.jet) # @UndefinedVariable
ax.scatter(self.X[:, 0], self.X[:, 1], 40, self.likelihood.Y, linewidth=0, cmap=pb.cm.jet, vmin=m.min(), vmax=m.max()) # @UndefinedVariable
ax.set_xlim(xmin[0], xmax[0])
ax.set_ylim(xmin[1], xmax[1])
else:
@ -94,9 +94,9 @@ class GPBase(model.model):
def plot(self, plot_limits=None, which_data='all', which_parts='all', resolution=None, levels=20, samples=0, fignum=None, ax=None):
"""
TODO: Docstrings!
:param levels: for 2D plotting, the number of contour levels to use
is ax is None, create a new figure
"""
# TODO include samples
if which_data == 'all':
@ -111,7 +111,7 @@ class GPBase(model.model):
Xu = self.X * self._Xstd + self._Xmean # NOTE self.X are the normalized values now
Xnew, xmin, xmax = x_frame1D(Xu, plot_limits=plot_limits)
m, var, lower, upper = self.predict(Xnew, which_parts=which_parts)
m, _, lower, upper = self.predict(Xnew, which_parts=which_parts)
for d in range(m.shape[1]):
gpplot(Xnew, m[:, d], lower[:, d], upper[:, d], axes=ax)
ax.plot(Xu[which_data], self.likelihood.data[which_data, d], 'kx', mew=1.5)
@ -122,13 +122,13 @@ class GPBase(model.model):
elif self.X.shape[1] == 2: # FIXME
resolution = resolution or 50
Xnew, xx, yy, xmin, xmax = x_frame2D(self.X, plot_limits, resolution)
Xnew, _, _, xmin, xmax = x_frame2D(self.X, plot_limits, resolution)
x, y = np.linspace(xmin[0], xmax[0], resolution), np.linspace(xmin[1], xmax[1], resolution)
m, var, lower, upper = self.predict(Xnew, which_parts=which_parts)
m, _, lower, upper = self.predict(Xnew, which_parts=which_parts)
m = m.reshape(resolution, resolution).T
ax.contour(x, y, m, levels, vmin=m.min(), vmax=m.max(), cmap=pb.cm.jet)
ax.contour(x, y, m, levels, vmin=m.min(), vmax=m.max(), cmap=pb.cm.jet) # @UndefinedVariable
Yf = self.likelihood.Y.flatten()
ax.scatter(self.X[:, 0], self.X[:, 1], 40, Yf, cmap=pb.cm.jet, vmin=m.min(), vmax=m.max(), linewidth=0.)
ax.scatter(self.X[:, 0], self.X[:, 1], 40, Yf, cmap=pb.cm.jet, vmin=m.min(), vmax=m.max(), linewidth=0.) # @UndefinedVariable
ax.set_xlim(xmin[0], xmax[0])
ax.set_ylim(xmin[1], xmax[1])

55
GPy/core/model.py
View file

@ -6,37 +6,32 @@ from .. import likelihoods
from ..inference import optimization
from ..util.linalg import jitchol
from GPy.util.misc import opt_wrapper
from parameterised import parameterised
from scipy import optimize
from parameterised import Parameterised
import multiprocessing as mp
import numpy as np
import priors
import re
import sys
import pdb
from GPy.core.domains import POSITIVE, REAL
# import numdifftools as ndt
class model(parameterised):
class Model(Parameterised):
def __init__(self):
parameterised.__init__(self)
Parameterised.__init__(self)
self.priors = None
self.optimization_runs = []
self.sampling_runs = []
self.preferred_optimizer = 'scg'
# self._set_params(self._get_params()) has been taken out as it should only be called on leaf nodes
def _get_params(self):
raise NotImplementedError, "this needs to be implemented to use the model class"
raise NotImplementedError, "this needs to be implemented to use the Model class"
def _set_params(self, x):
raise NotImplementedError, "this needs to be implemented to use the model class"
raise NotImplementedError, "this needs to be implemented to use the Model class"
def log_likelihood(self):
raise NotImplementedError, "this needs to be implemented to use the model class"
raise NotImplementedError, "this needs to be implemented to use the Model class"
def _log_likelihood_gradients(self):
raise NotImplementedError, "this needs to be implemented to use the model class"
raise NotImplementedError, "this needs to be implemented to use the Model class"
def set_prior(self, regexp, what):
"""
Sets priors on the model parameters.
Sets priors on the Model parameters.
Arguments
---------
@ -95,7 +90,7 @@ class model(parameterised):
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.
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):
@ -135,7 +130,7 @@ class model(parameterised):
def randomize(self):
"""
Randomize the model.
Randomize the Model.
Make this draw from the Prior if one exists, else draw from N(0,1)
"""
# first take care of all parameters (from N(0,1))
@ -152,11 +147,11 @@ class model(parameterised):
def optimize_restarts(self, num_restarts=10, robust=False, verbose=True, parallel=False, num_processes=None, **kwargs):
"""
Perform random restarts of the model, and set the model to the best
Perform random restarts of the Model, and set the Model to the best
seen solution.
If the robust flag is set, exceptions raised during optimizations will
be handled silently. If _all_ runs fail, the model is reset to the
be handled silently. If _all_ runs fail, the Model is reset to the
existing parameter values.
Notes
@ -218,7 +213,7 @@ class model(parameterised):
Ensure that any variables which should clearly be positive have been constrained somehow.
"""
positive_strings = ['variance', 'lengthscale', 'precision', 'kappa']
param_names = self._get_param_names()
# param_names = self._get_param_names()
currently_constrained = self.all_constrained_indices()
to_make_positive = []
for s in positive_strings:
@ -251,7 +246,7 @@ class model(parameterised):
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.
Optimize the Model using self.log_likelihood and self.log_likelihood_gradient, as well as self.priors.
kwargs are passed to the optimizer. They can be:
:max_f_eval: maximum number of function evaluations
@ -274,7 +269,7 @@ class model(parameterised):
def optimize_SGD(self, momentum=0.1, learning_rate=0.01, iterations=20, **kwargs):
# assert self.Y.shape[1] > 1, "SGD only works with D > 1"
sgd = SGD.StochasticGD(self, iterations, learning_rate, momentum, **kwargs)
sgd = SGD.StochasticGD(self, iterations, learning_rate, momentum, **kwargs) # @UndefinedVariable
sgd.run()
self.optimization_runs.append(sgd)
@ -291,7 +286,7 @@ class model(parameterised):
def f(x):
self._set_params(x)
return self.log_likelihood()
h = ndt.Hessian(f)
h = ndt.Hessian(f) # @UndefinedVariable
A = -h(x)
self._set_params(x)
# check for almost zero components on the diagonal which screw up the cholesky
@ -300,7 +295,7 @@ class model(parameterised):
return A
def Laplace_evidence(self):
"""Returns an estiamte of the model evidence based on the Laplace approximation.
"""Returns an estiamte of the Model evidence based on the Laplace approximation.
Uses a numerical estimate of the hessian if none is available analytically"""
A = self.Laplace_covariance()
try:
@ -310,7 +305,7 @@ class model(parameterised):
return 0.5 * self._get_params().size * np.log(2 * np.pi) + self.log_likelihood() - hld
def __str__(self):
s = parameterised.__str__(self).split('\n')
s = Parameterised.__str__(self).split('\n')
# add priors to the string
if self.priors is not None:
strs = [str(p) if p is not None else '' for p in self.priors]
@ -336,7 +331,7 @@ class model(parameterised):
def checkgrad(self, target_param=None, verbose=False, step=1e-6, tolerance=1e-3):
"""
Check the gradient of the model by comparing to a numerical estimate.
Check the gradient of the Model by comparing to a numerical estimate.
If the verbose flag is passed, invividual components are tested (and printed)
:param verbose: If True, print a "full" checking of each parameter
@ -389,7 +384,7 @@ class model(parameterised):
param_list = range(len(x))
else:
param_list = self.grep_param_names(target_param, transformed=True, search=True)
if not param_list:
if not np.any(param_list):
print "No free parameters to check"
return
@ -419,15 +414,15 @@ class model(parameterised):
def input_sensitivity(self):
"""
return an array describing the sesitivity of the model to each input
return an array describing the sesitivity of the Model to each input
NB. Right now, we're basing this on the lengthscales (or
variances) of the kernel. TODO: proper sensitivity analysis
where we integrate across the model inputs and evaluate the
effect on the variance of the model output. """
where we integrate across the Model inputs and evaluate the
effect on the variance of the Model output. """
if not hasattr(self, 'kern'):
raise ValueError, "this model has no kernel"
raise ValueError, "this Model has no kernel"
k = [p for p in self.kern.parts if p.name in ['rbf', 'linear']]
if (not len(k) == 1) or (not k[0].ARD):
@ -475,7 +470,7 @@ class model(parameterised):
if ll_change < 0:
self.likelihood = last_approximation # restore previous likelihood approximation
self._set_params(last_params) # restore model parameters
self._set_params(last_params) # restore Model parameters
print "Log-likelihood decrement: %s \nLast likelihood update discarded." % ll_change
stop = True
else:

6
GPy/core/parameterised.py
View file

@ -9,7 +9,7 @@ import cPickle
import warnings
import transformations
class parameterised(object):
class Parameterised(object):
def __init__(self):
"""
This is the base class for model and kernel. Mostly just handles tieing and constraining of parameters
@ -34,7 +34,7 @@ class parameterised(object):
"""
Returns a **copy** of parameters in non transformed space
:see_also: :py:func:`GPy.core.parameterised.params_transformed`
:see_also: :py:func:`GPy.core.Parameterised.params_transformed`
"""
return self._get_params()
@ -47,7 +47,7 @@ class parameterised(object):
"""
Returns a **copy** of parameters in transformed space
:see_also: :py:func:`GPy.core.parameterised.params`
:see_also: :py:func:`GPy.core.Parameterised.params`
"""
return self._get_params_transformed()

2
GPy/core/sparse_gp.py
View file

@ -152,7 +152,7 @@ class SparseGP(GPBase):
return A + B + C + D + self.likelihood.Z
def _set_params(self, p):
self.Z = p[:self.num_inducing * self.output_dim].reshape(self.num_inducing, self.input_dim)
self.Z = p[:self.num_inducing * self.input_dim].reshape(self.num_inducing, self.input_dim)
self.kern._set_params(p[self.Z.size:self.Z.size + self.kern.num_params])
self.likelihood._set_params(p[self.Z.size + self.kern.num_params:])
self._compute_kernel_matrices()