Mods to regression.py now that we have get to get parameters. Moved Youter to YYT.

This commit is contained in:
Neil Lawrence 2013-01-18 15:31:20 +00:00
parent 11d088cf90
commit 99034d0fb0
6 changed files with 43 additions and 38 deletions

View file

@ -6,5 +6,6 @@ import kern
import models import models
import inference import inference
import util import util
import examples
#import examples TODO: discuss! #import examples TODO: discuss!
from core import priors from core import priors

View file

@ -80,19 +80,22 @@ class model(parameterised):
for w in which: for w in which:
self.priors[w] = what 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) matches = self.grep_param_names(name)
if len(matches): if len(matches):
return self._get_params()[matches] if return_names:
return self._get_params()[matches], np.asarray(self._get_param_names())[matches].tolist()
else:
return self._get_params()[matches]
else: else:
raise AttributeError, "no parameter matches %s"%name raise AttributeError, "no parameter matches %s"%name
def set(self,name,val): 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) matches = self.grep_param_names(name)
if len(matches): if len(matches):
@ -102,6 +105,20 @@ class model(parameterised):
else: else:
raise AttributeError, "no parameter matches %s"%name 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): def log_prior(self):

View file

@ -17,10 +17,8 @@ def toy_rbf_1d():
# create simple GP model # create simple GP model
m = GPy.models.GP_regression(data['X'],data['Y']) m = GPy.models.GP_regression(data['X'],data['Y'])
# contrain all parameters to be positive
m.constrain_positive('')
# optimize # optimize
m.ensure_default_constraints()
m.optimize() m.optimize()
# plot # plot
@ -35,10 +33,8 @@ def rogers_girolami_olympics():
# create simple GP model # create simple GP model
m = GPy.models.GP_regression(data['X'],data['Y']) m = GPy.models.GP_regression(data['X'],data['Y'])
# contrain all parameters to be positive
m.constrain_positive('')
# optimize # optimize
m.ensure_default_constraints()
m.optimize() m.optimize()
# plot # plot
@ -57,10 +53,8 @@ def toy_rbf_1d_50():
# create simple GP model # create simple GP model
m = GPy.models.GP_regression(data['X'],data['Y']) m = GPy.models.GP_regression(data['X'],data['Y'])
# contrain all parameters to be positive
m.constrain_positive('')
# optimize # optimize
m.ensure_default_constraints()
m.optimize() m.optimize()
# plot # plot
@ -75,10 +69,8 @@ def silhouette():
# create simple GP model # create simple GP model
m = GPy.models.GP_regression(data['X'],data['Y']) m = GPy.models.GP_regression(data['X'],data['Y'])
# contrain all parameters to be positive
m.constrain_positive('')
# optimize # optimize
m.ensure_default_constraints()
m.optimize() m.optimize()
print(m) 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.)) 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) m = GPy.models.GP_regression(data['X'],data['Y'], kernel=kern)
params = m._get_params() optim_point_x[0] = m.get('rbf_lengthscale')
optim_point_x[0] = params[1] optim_point_y[0] = np.log10(m.get('rbf_variance')) - np.log10(m.get('white_variance'));
optim_point_y[0] = np.log10(params[0]) - np.log10(params[2]);
# contrain all parameters to be positive
m.constrain_positive('')
# optimize # optimize
m.ensure_default_constraints()
m.optimize(xtol=1e-6,ftol=1e-6) m.optimize(xtol=1e-6,ftol=1e-6)
params = m._get_params() optim_point_x[1] = m.get('rbf_lengthscale')
optim_point_x[1] = params[1] optim_point_y[1] = np.log10(m.get('rbf_variance')) - np.log10(m.get('white_variance'));
optim_point_y[1] = np.log10(params[0]) - np.log10(params[2]);
print(m)
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') 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) models.append(m)

View file

@ -63,10 +63,10 @@ class GP_regression(model):
self._Ystd = np.ones((1,self.Y.shape[1])) self._Ystd = np.ones((1,self.Y.shape[1]))
if self.D > self.N: if self.D > self.N:
# then it's more efficient to store Youter # then it's more efficient to store YYT
self.Youter = np.dot(self.Y, self.Y.T) self.YYT = np.dot(self.Y, self.Y.T)
else: else:
self.Youter = None self.YYT = None
model.__init__(self) model.__init__(self)
@ -83,23 +83,23 @@ class GP_regression(model):
def _model_fit_term(self): 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))) return -0.5*np.sum(np.square(np.dot(self.Li,self.Y)))
else: 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): def log_likelihood(self):
complexity_term = -0.5*self.N*self.D*np.log(2.*np.pi) - 0.5*self.D*self.K_logdet 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() return complexity_term + self._model_fit_term()
def dL_dK(self): def dL_dK(self):
if self.Youter is None: if self.YYT is None:
alpha = np.dot(self.Ki,self.Y) alpha = np.dot(self.Ki,self.Y)
dL_dK = 0.5*(np.dot(alpha,alpha.T)-self.D*self.Ki) dL_dK = 0.5*(np.dot(alpha,alpha.T)-self.D*self.Ki)
else: 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 return dL_dK

View file

@ -91,9 +91,9 @@ class generalized_FITC(model):
def log_likelihood(self): def log_likelihood(self):
self.posterior_param() 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 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)) 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_)) 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)) 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))

View file

@ -48,9 +48,9 @@ class warpedGP(GP_regression):
# this supports the 'smart' behaviour in GP_regression # this supports the 'smart' behaviour in GP_regression
if self.D > self.N: if self.D > self.N:
self.Youter = np.dot(self.Y, self.Y.T) self.YYT = np.dot(self.Y, self.Y.T)
else: else:
self.Youter = None self.YYT = None
return self.Y return self.Y