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

Fixed likelihood tests for new parameters structure

Browse source
This commit is contained in:
Alan Saul 2014-02-07 15:16:52 +00:00
parent c28f11f291
commit 186feb45a1
6 changed files with 207 additions and 168 deletions
Download patch file Download diff file Expand all files Collapse all files

10
GPy/inference/latent_function_inference/laplace.py
View file

@ -33,7 +33,6 @@ class LaplaceInference(object):
self._mode_finding_max_iter = 40 self._mode_finding_max_iter = 40
self.bad_fhat = True self.bad_fhat = True
def inference(self, kern, X, likelihood, Y, Y_metadata=None): def inference(self, kern, X, likelihood, Y, Y_metadata=None):
""" """
Returns a Posterior class containing essential quantities of the posterior Returns a Posterior class containing essential quantities of the posterior
@ -50,6 +49,7 @@ class LaplaceInference(object):
Ki_f_init = np.zeros_like(Y) Ki_f_init = np.zeros_like(Y)
else: else:
Ki_f_init = self._previous_Ki_fhat Ki_f_init = self._previous_Ki_fhat
f_hat, Ki_fhat = self.rasm_mode(K, Y, likelihood, Ki_f_init, Y_metadata=Y_metadata) f_hat, Ki_fhat = self.rasm_mode(K, Y, likelihood, Ki_f_init, Y_metadata=Y_metadata)
#Compute hessian and other variables at mode #Compute hessian and other variables at mode
@ -57,6 +57,7 @@ class LaplaceInference(object):
#likelihood.gradient = self.likelihood_gradients() #likelihood.gradient = self.likelihood_gradients()
kern.update_gradients_full(dL_dK, X) kern.update_gradients_full(dL_dK, X)
likelihood.update_gradients(np.ones(10))
self._previous_Ki_fhat = Ki_fhat.copy() self._previous_Ki_fhat = Ki_fhat.copy()
return Posterior(woodbury_vector=woodbury_vector, woodbury_inv = K_Wi_i, K=K), log_marginal, {'dL_dK':dL_dK} return Posterior(woodbury_vector=woodbury_vector, woodbury_inv = K_Wi_i, K=K), log_marginal, {'dL_dK':dL_dK}
@ -157,9 +158,12 @@ class LaplaceInference(object):
explicit_part = 0.5*(np.dot(Ki_f, Ki_f.T) - K_Wi_i) explicit_part = 0.5*(np.dot(Ki_f, Ki_f.T) - K_Wi_i)
#Implicit #Implicit
d3lik_d3fhat = likelihood.d3logpdf_df3(f_hat, Y, extra_data=Y_metadata) dW_df = likelihood.d3logpdf_df3(f_hat, Y, extra_data=Y_metadata) # d3lik_d3fhat
dL_dfhat = 0.5*(np.diag(Ki_W_i)[:, None]*d3lik_d3fhat) #why isn't this -0.5? s2 in R&W p126 line 9.
woodbury_vector = likelihood.dlogpdf_df(f_hat, Y, extra_data=Y_metadata) woodbury_vector = likelihood.dlogpdf_df(f_hat, Y, extra_data=Y_metadata)
dL_dfhat = 0.5*(np.diag(Ki_W_i)[:, None]*dW_df) #why isn't this -0.5? s2 in R&W p126 line 9.
#implicit_part = np.dot(woodbury_vector, dL_dfhat.T).dot(np.eye(Y.shape[0]) - np.dot(K, K_Wi_i))
BiK, _ = dpotrs(L, K, lower=1)
#dL_dfhat = 0.5*np.diag(BiK)[:, None]*dW_df
implicit_part = np.dot(woodbury_vector, dL_dfhat.T).dot(np.eye(Y.shape[0]) - np.dot(K, K_Wi_i)) implicit_part = np.dot(woodbury_vector, dL_dfhat.T).dot(np.eye(Y.shape[0]) - np.dot(K, K_Wi_i))
dL_dK = explicit_part + implicit_part dL_dK = explicit_part + implicit_part

17
GPy/likelihoods/gaussian.py
View file

@ -130,7 +130,10 @@ class Gaussian(Likelihood):
:rtype: float :rtype: float
""" """
assert np.asarray(link_f).shape == np.asarray(y).shape assert np.asarray(link_f).shape == np.asarray(y).shape
return -0.5*(np.sum((y-link_f)**2/self.variance) + self.ln_det_K + self.N*np.log(2.*np.pi)) N = y.shape[0]
ln_det_cov = N*np.log(self.variance)
return -0.5*(np.sum((y-link_f)**2/self.variance) + ln_det_cov + N*np.log(2.*np.pi))
def dlogpdf_dlink(self, link_f, y, extra_data=None): def dlogpdf_dlink(self, link_f, y, extra_data=None):
""" """
@ -175,7 +178,8 @@ class Gaussian(Likelihood):
(the distribution for y_i depends only on link(f_i) not on link(f_(j!=i)) (the distribution for y_i depends only on link(f_i) not on link(f_(j!=i))
""" """
assert np.asarray(link_f).shape == np.asarray(y).shape assert np.asarray(link_f).shape == np.asarray(y).shape
hess = -(1.0/self.variance)*np.ones((self.N, 1)) N = y.shape[0]
hess = -(1.0/self.variance)*np.ones((N, 1))
return hess return hess
def d3logpdf_dlink3(self, link_f, y, extra_data=None): def d3logpdf_dlink3(self, link_f, y, extra_data=None):
@ -194,7 +198,8 @@ class Gaussian(Likelihood):
:rtype: Nx1 array :rtype: Nx1 array
""" """
assert np.asarray(link_f).shape == np.asarray(y).shape assert np.asarray(link_f).shape == np.asarray(y).shape
d3logpdf_dlink3 = np.diagonal(0*self.I)[:, None] N = y.shape[0]
d3logpdf_dlink3 = np.zeros((N,1))
return d3logpdf_dlink3 return d3logpdf_dlink3
def dlogpdf_link_dvar(self, link_f, y, extra_data=None): def dlogpdf_link_dvar(self, link_f, y, extra_data=None):
@ -215,7 +220,8 @@ class Gaussian(Likelihood):
assert np.asarray(link_f).shape == np.asarray(y).shape assert np.asarray(link_f).shape == np.asarray(y).shape
e = y - link_f e = y - link_f
s_4 = 1.0/(self.variance**2) s_4 = 1.0/(self.variance**2)
dlik_dsigma = -0.5*self.N/self.variance + 0.5*s_4*np.sum(np.square(e)) N = y.shape[0]
dlik_dsigma = -0.5*N/self.variance + 0.5*s_4*np.sum(np.square(e))
return np.sum(dlik_dsigma) # Sure about this sum? return np.sum(dlik_dsigma) # Sure about this sum?
def dlogpdf_dlink_dvar(self, link_f, y, extra_data=None): def dlogpdf_dlink_dvar(self, link_f, y, extra_data=None):
@ -255,7 +261,8 @@ class Gaussian(Likelihood):
""" """
assert np.asarray(link_f).shape == np.asarray(y).shape assert np.asarray(link_f).shape == np.asarray(y).shape
s_4 = 1.0/(self.variance**2) s_4 = 1.0/(self.variance**2)
d2logpdf_dlink2_dvar = np.diag(s_4*self.I)[:, None] N = y.shape[0]
d2logpdf_dlink2_dvar = np.ones((N,1))*s_4
return d2logpdf_dlink2_dvar return d2logpdf_dlink2_dvar
def dlogpdf_link_dtheta(self, f, y, extra_data=None): def dlogpdf_link_dtheta(self, f, y, extra_data=None):

16
GPy/likelihoods/likelihood.py
View file

@ -44,6 +44,10 @@ class Likelihood(Parameterized):
def _gradients(self,partial): def _gradients(self,partial):
return np.zeros(0) return np.zeros(0)
def update_gradients(self, partial):
if self.size > 0:
raise NotImplementedError('Must be implemented for likelihoods with parameters to be optimized')
def _preprocess_values(self,Y): def _preprocess_values(self,Y):
""" """
In case it is needed, this function assess the output values or makes any pertinent transformation on them. In case it is needed, this function assess the output values or makes any pertinent transformation on them.
@ -303,7 +307,7 @@ class Likelihood(Parameterized):
""" """
TODO: Doc strings TODO: Doc strings
""" """
if len(self._get_param_names()) > 0: if self.size > 0:
link_f = self.gp_link.transf(f) link_f = self.gp_link.transf(f)
return self.dlogpdf_link_dtheta(link_f, y, extra_data=extra_data) return self.dlogpdf_link_dtheta(link_f, y, extra_data=extra_data)
else: else:
@ -314,7 +318,7 @@ class Likelihood(Parameterized):
""" """
TODO: Doc strings TODO: Doc strings
""" """
if len(self._get_param_names()) > 0: if self.size > 0:
link_f = self.gp_link.transf(f) link_f = self.gp_link.transf(f)
dlink_df = self.gp_link.dtransf_df(f) dlink_df = self.gp_link.dtransf_df(f)
dlogpdf_dlink_dtheta = self.dlogpdf_dlink_dtheta(link_f, y, extra_data=extra_data) dlogpdf_dlink_dtheta = self.dlogpdf_dlink_dtheta(link_f, y, extra_data=extra_data)
@ -327,7 +331,7 @@ class Likelihood(Parameterized):
""" """
TODO: Doc strings TODO: Doc strings
""" """
if len(self._get_param_names()) > 0: if self.size > 0:
link_f = self.gp_link.transf(f) link_f = self.gp_link.transf(f)
dlink_df = self.gp_link.dtransf_df(f) dlink_df = self.gp_link.dtransf_df(f)
d2link_df2 = self.gp_link.d2transf_df2(f) d2link_df2 = self.gp_link.d2transf_df2(f)
@ -345,9 +349,9 @@ class Likelihood(Parameterized):
#Parameters are stacked vertically. Must be listed in same order as 'get_param_names' #Parameters are stacked vertically. Must be listed in same order as 'get_param_names'
# ensure we have gradients for every parameter we want to optimize # ensure we have gradients for every parameter we want to optimize
assert dlogpdf_dtheta.shape[1] == len(self._get_param_names()) assert dlogpdf_dtheta.shape[1] == self.size
assert dlogpdf_df_dtheta.shape[1] == len(self._get_param_names()) assert dlogpdf_df_dtheta.shape[1] == self.size
assert d2logpdf_df2_dtheta.shape[1] == len(self._get_param_names()) assert d2logpdf_df2_dtheta.shape[1] == self.size
return dlogpdf_dtheta, dlogpdf_df_dtheta, d2logpdf_df2_dtheta return dlogpdf_dtheta, dlogpdf_df_dtheta, d2logpdf_df2_dtheta

36
GPy/likelihoods/student_t.py
View file

@ -8,6 +8,7 @@ import link_functions
from scipy import stats, integrate from scipy import stats, integrate
from scipy.special import gammaln, gamma from scipy.special import gammaln, gamma
from likelihood import Likelihood from likelihood import Likelihood
from ..core.parameterization import Param
class StudentT(Likelihood): class StudentT(Likelihood):
""" """
@ -19,26 +20,25 @@ class StudentT(Likelihood):
p(y_{i}|\\lambda(f_{i})) = \\frac{\\Gamma\\left(\\frac{v+1}{2}\\right)}{\\Gamma\\left(\\frac{v}{2}\\right)\\sqrt{v\\pi\\sigma^{2}}}\\left(1 + \\frac{1}{v}\\left(\\frac{(y_{i} - f_{i})^{2}}{\\sigma^{2}}\\right)\\right)^{\\frac{-v+1}{2}} p(y_{i}|\\lambda(f_{i})) = \\frac{\\Gamma\\left(\\frac{v+1}{2}\\right)}{\\Gamma\\left(\\frac{v}{2}\\right)\\sqrt{v\\pi\\sigma^{2}}}\\left(1 + \\frac{1}{v}\\left(\\frac{(y_{i} - f_{i})^{2}}{\\sigma^{2}}\\right)\\right)^{\\frac{-v+1}{2}}
""" """
def __init__(self,gp_link=None,analytical_mean=True,analytical_variance=True, deg_free=5, sigma2=2): def __init__(self,gp_link=None, deg_free=5, sigma2=2):
self.v = deg_free if gp_link is None:
self.sigma2 = sigma2 gp_link = link_functions.Identity()
super(StudentT, self).__init__(gp_link, name='Student_T')
self.sigma2 = Param('t_noise', float(sigma2))
self.v = Param('deg_free', float(deg_free))
self.add_parameter(self.sigma2)
self.add_parameter(self.v)
self._set_params(np.asarray(sigma2))
super(StudentT, self).__init__(gp_link,analytical_mean,analytical_variance)
self.log_concave = False self.log_concave = False
def _get_params(self): def parameters_changed(self):
return np.asarray(self.sigma2) self.variance = (self.v / float(self.v - 2)) * self.sigma2
def _get_param_names(self): def update_gradients(self, partial):
return ["t_noise_std2"] self.sigma2.gradient = np.ones(1) #FIXME: Not done yet
self.v.gradient = np.ones(1) #FIXME: Not done yet
def _set_params(self, x):
self.sigma2 = float(x)
@property
def variance(self, extra_data=None):
return (self.v / float(self.v - 2)) * self.sigma2
def pdf_link(self, link_f, y, extra_data=None): def pdf_link(self, link_f, y, extra_data=None):
""" """
@ -82,6 +82,10 @@ class StudentT(Likelihood):
""" """
assert np.atleast_1d(link_f).shape == np.atleast_1d(y).shape assert np.atleast_1d(link_f).shape == np.atleast_1d(y).shape
e = y - link_f e = y - link_f
#FIXME:
#Why does np.log(1 + (1/self.v)*((y-link_f)**2)/self.sigma2) suppress the divide by zero?!
#But np.log(1 + (1/float(self.v))*((y-link_f)**2)/self.sigma2) throws it correctly
#print - 0.5*(self.v + 1)*np.log(1 + (1/np.float(self.v))*((e**2)/self.sigma2))
objective = (+ gammaln((self.v + 1) * 0.5) objective = (+ gammaln((self.v + 1) * 0.5)
- gammaln(self.v * 0.5) - gammaln(self.v * 0.5)
- 0.5*np.log(self.sigma2 * self.v * np.pi) - 0.5*np.log(self.sigma2 * self.v * np.pi)

37
GPy/models/gradient_checker.py
View file

@ -1,9 +1,10 @@
# ## Copyright (c) 2012, GPy authors (see AUTHORS.txt). # ## Copyright (c) 2012, GPy authors (see AUTHORS.txt).
# Licensed under the BSD 3-clause license (see LICENSE.txt) # Licensed under the BSD 3-clause license (see LICENSE.txt)
from GPy.core.model import Model from ..core.model import Model
import itertools import itertools
import numpy import numpy
from ..core.parameterization import Param
def get_shape(x): def get_shape(x):
if isinstance(x, numpy.ndarray): if isinstance(x, numpy.ndarray):
@ -59,7 +60,7 @@ class GradientChecker(Model):
grad.randomize() grad.randomize()
grad.checkgrad(verbose=1) grad.checkgrad(verbose=1)
""" """
Model.__init__(self) Model.__init__(self, 'GradientChecker')
if isinstance(x0, (list, tuple)) and names is None: if isinstance(x0, (list, tuple)) and names is None:
self.shapes = [get_shape(xi) for xi in x0] self.shapes = [get_shape(xi) for xi in x0]
self.names = ['X{i}'.format(i=i) for i in range(len(x0))] self.names = ['X{i}'.format(i=i) for i in range(len(x0))]
@ -72,8 +73,10 @@ class GradientChecker(Model):
else: else:
self.names = names self.names = names
self.shapes = [get_shape(x0)] self.shapes = [get_shape(x0)]
for name, xi in zip(self.names, at_least_one_element(x0)): for name, xi in zip(self.names, at_least_one_element(x0)):
self.__setattr__(name, xi) self.__setattr__(name, Param(name, xi))
self.add_parameter(self.__getattribute__(name))
# self._param_names = [] # self._param_names = []
# for name, shape in zip(self.names, self.shapes): # for name, shape in zip(self.names, self.shapes):
# self._param_names.extend(map(lambda nameshape: ('_'.join(nameshape)).strip('_'), itertools.izip(itertools.repeat(name), itertools.imap(lambda t: '_'.join(map(str, t)), itertools.product(*map(lambda xi: range(xi), shape)))))) # self._param_names.extend(map(lambda nameshape: ('_'.join(nameshape)).strip('_'), itertools.izip(itertools.repeat(name), itertools.imap(lambda t: '_'.join(map(str, t)), itertools.product(*map(lambda xi: range(xi), shape))))))
@ -93,20 +96,18 @@ class GradientChecker(Model):
def _log_likelihood_gradients(self): def _log_likelihood_gradients(self):
return numpy.atleast_1d(self.df(*self._get_x(), **self.kwargs)).flatten() return numpy.atleast_1d(self.df(*self._get_x(), **self.kwargs)).flatten()
#def _get_params(self):
#return numpy.atleast_1d(numpy.hstack(map(lambda name: flatten_if_needed(self.__getattribute__(name)), self.names)))
def _get_params(self): #def _set_params(self, x):
return numpy.atleast_1d(numpy.hstack(map(lambda name: flatten_if_needed(self.__getattribute__(name)), self.names))) #current_index = 0
#for name, shape in zip(self.names, self.shapes):
#current_size = numpy.prod(shape)
#self.__setattr__(name, x[current_index:current_index + current_size].reshape(shape))
#current_index += current_size
#def _get_param_names(self):
def _set_params(self, x): #_param_names = []
current_index = 0 #for name, shape in zip(self.names, self.shapes):
for name, shape in zip(self.names, self.shapes): #_param_names.extend(map(lambda nameshape: ('_'.join(nameshape)).strip('_'), itertools.izip(itertools.repeat(name), itertools.imap(lambda t: '_'.join(map(str, t)), itertools.product(*map(lambda xi: range(xi), shape))))))
current_size = numpy.prod(shape) #return _param_names
self.__setattr__(name, x[current_index:current_index + current_size].reshape(shape))
current_index += current_size
def _get_param_names(self):
_param_names = []
for name, shape in zip(self.names, self.shapes):
_param_names.extend(map(lambda nameshape: ('_'.join(nameshape)).strip('_'), itertools.izip(itertools.repeat(name), itertools.imap(lambda t: '_'.join(map(str, t)), itertools.product(*map(lambda xi: range(xi), shape))))))
return _param_names

217
GPy/testing/likelihood_tests.py
View file

@ -4,7 +4,8 @@ import GPy
from GPy.models import GradientChecker from GPy.models import GradientChecker
import functools import functools
import inspect import inspect
from GPy.likelihoods.noise_models import gp_transformations from GPy.likelihoods import link_functions
from ..core.parameterization import Param
from functools import partial from functools import partial
#np.random.seed(300) #np.random.seed(300)
np.random.seed(7) np.random.seed(7)
@ -22,12 +23,14 @@ def dparam_partial(inst_func, *args):
the f or Y that are being used in the function whilst we tweak the the f or Y that are being used in the function whilst we tweak the
param param
""" """
def param_func(param, inst_func, args): def param_func(param_val, param_name, inst_func, args):
inst_func.im_self._set_params(param) #inst_func.im_self._set_params(param)
#inst_func.im_self.add_parameter(Param(param_name, param_val))
inst_func.im_self[param_name] = param_val
return inst_func(*args) return inst_func(*args)
return functools.partial(param_func, inst_func=inst_func, args=args) return functools.partial(param_func, inst_func=inst_func, args=args)
def dparam_checkgrad(func, dfunc, params, args, constraints=None, randomize=False, verbose=False): def dparam_checkgrad(func, dfunc, params, params_names, args, constraints=None, randomize=False, verbose=False):
""" """
checkgrad expects a f: R^N -> R^1 and df: R^N -> R^N checkgrad expects a f: R^N -> R^1 and df: R^N -> R^N
However if we are holding other parameters fixed and moving something else However if we are holding other parameters fixed and moving something else
@ -43,22 +46,27 @@ def dparam_checkgrad(func, dfunc, params, args, constraints=None, randomize=Fals
partial_f = dparam_partial(func, *args) partial_f = dparam_partial(func, *args)
partial_df = dparam_partial(dfunc, *args) partial_df = dparam_partial(dfunc, *args)
gradchecking = True gradchecking = True
for param in params: zipped_params = zip(params, params_names)
fnum = np.atleast_1d(partial_f(param)).shape[0] for param_val, param_name in zipped_params:
dfnum = np.atleast_1d(partial_df(param)).shape[0] fnum = np.atleast_1d(partial_f(param_val, param_name)).shape[0]
dfnum = np.atleast_1d(partial_df(param_val, param_name)).shape[0]
for fixed_val in range(dfnum): for fixed_val in range(dfnum):
#dlik and dlik_dvar gives back 1 value for each #dlik and dlik_dvar gives back 1 value for each
f_ind = min(fnum, fixed_val+1) - 1 f_ind = min(fnum, fixed_val+1) - 1
print "fnum: {} dfnum: {} f_ind: {} fixed_val: {}".format(fnum, dfnum, f_ind, fixed_val) print "fnum: {} dfnum: {} f_ind: {} fixed_val: {}".format(fnum, dfnum, f_ind, fixed_val)
#Make grad checker with this param moving, note that set_params is NOT being called #Make grad checker with this param moving, note that set_params is NOT being called
#The parameter is being set directly with __setattr__ #The parameter is being set directly with __setattr__
grad = GradientChecker(lambda x: np.atleast_1d(partial_f(x))[f_ind], grad = GradientChecker(lambda p_val: np.atleast_1d(partial_f(p_val, param_name))[f_ind],
lambda x : np.atleast_1d(partial_df(x))[fixed_val], lambda p_val: np.atleast_1d(partial_df(p_val, param_name))[fixed_val],
param, 'p') param_val, [param_name])
#This is not general for more than one param... #This is not general for more than one param...
if constraints is not None: if constraints is not None:
for constraint in constraints: for constrain_param, constraint in constraints:
constraint('p', grad) if grad.grep_param_names(constrain_param):
constraint(constrain_param, grad)
else:
print "parameter didn't exist"
print constrain_param, " ", constraint
if randomize: if randomize:
grad.randomize() grad.randomize()
if verbose: if verbose:
@ -107,17 +115,20 @@ class TestNoiseModels(object):
#################################################### ####################################################
# Constraint wrappers so we can just list them off # # Constraint wrappers so we can just list them off #
#################################################### ####################################################
def constrain_fixed(regex, model, value):
model[regex].constrain_fixed(value)
def constrain_negative(regex, model): def constrain_negative(regex, model):
model.constrain_negative(regex) model[regex].constrain_negative()
def constrain_positive(regex, model): def constrain_positive(regex, model):
model.constrain_positive(regex) model[regex].constrain_positive()
def constrain_bounded(regex, model, lower, upper): def constrain_bounded(regex, model, lower, upper):
""" """
Used like: partial(constrain_bounded, lower=0, upper=1) Used like: partial(constrain_bounded, lower=0, upper=1)
""" """
model.constrain_bounded(regex, lower, upper) model[regex].constrain_bounded(lower, upper)
""" """
Dictionary where we nest models we would like to check Dictionary where we nest models we would like to check
@ -134,71 +145,72 @@ class TestNoiseModels(object):
} }
""" """
noise_models = {"Student_t_default": { noise_models = {"Student_t_default": {
"model": GPy.likelihoods.student_t(deg_free=5, sigma2=self.var), "model": GPy.likelihoods.StudentT(deg_free=5, sigma2=self.var),
"grad_params": { "grad_params": {
"names": ["t_noise"], "names": ["t_noise"],
"vals": [self.var], "vals": [self.var],
"constraints": [constrain_positive] "constraints": [("t_noise", constrain_positive), ("deg_free", constrain_positive)]
#"constraints": [("t_noise", constrain_positive), ("deg_free", partial(constrain_fixed, value=5))]
}, },
"laplace": True "laplace": True
}, },
"Student_t_1_var": { "Student_t_1_var": {
"model": GPy.likelihoods.student_t(deg_free=5, sigma2=self.var), "model": GPy.likelihoods.StudentT(deg_free=5, sigma2=self.var),
"grad_params": { "grad_params": {
"names": ["t_noise"], "names": ["t_noise"],
"vals": [1.0], "vals": [1.0],
"constraints": [constrain_positive] "constraints": [("t_noise", constrain_positive), ("deg_free", constrain_positive)]
}, },
"laplace": True "laplace": True
}, },
"Student_t_small_var": { "Student_t_small_var": {
"model": GPy.likelihoods.student_t(deg_free=5, sigma2=self.var), "model": GPy.likelihoods.StudentT(deg_free=5, sigma2=self.var),
"grad_params": { "grad_params": {
"names": ["t_noise"], "names": ["t_noise"],
"vals": [0.01], "vals": [0.01],
"constraints": [constrain_positive] "constraints": [("t_noise", constrain_positive), ("deg_free", constrain_positive)]
}, },
"laplace": True "laplace": True
}, },
"Student_t_large_var": { "Student_t_large_var": {
"model": GPy.likelihoods.student_t(deg_free=5, sigma2=self.var), "model": GPy.likelihoods.StudentT(deg_free=5, sigma2=self.var),
"grad_params": { "grad_params": {
"names": ["t_noise"], "names": ["t_noise"],
"vals": [10.0], "vals": [10.0],
"constraints": [constrain_positive] "constraints": [("t_noise", constrain_positive), ("deg_free", constrain_positive)]
}, },
"laplace": True "laplace": True
}, },
"Student_t_approx_gauss": { "Student_t_approx_gauss": {
"model": GPy.likelihoods.student_t(deg_free=1000, sigma2=self.var), "model": GPy.likelihoods.StudentT(deg_free=1000, sigma2=self.var),
"grad_params": { "grad_params": {
"names": ["t_noise"], "names": ["t_noise"],
"vals": [self.var], "vals": [self.var],
"constraints": [constrain_positive] "constraints": [("t_noise", constrain_positive), ("deg_free", constrain_positive)]
}, },
"laplace": True "laplace": True
}, },
"Student_t_log": { "Student_t_log": {
"model": GPy.likelihoods.student_t(gp_link=gp_transformations.Log(), deg_free=5, sigma2=self.var), "model": GPy.likelihoods.StudentT(gp_link=link_functions.Log(), deg_free=5, sigma2=self.var),
"grad_params": { "grad_params": {
"names": ["t_noise"], "names": ["t_noise"],
"vals": [self.var], "vals": [self.var],
"constraints": [constrain_positive] "constraints": [("t_noise", constrain_positive), ("deg_free", constrain_positive)]
}, },
"laplace": True "laplace": True
}, },
"Gaussian_default": { "Gaussian_default": {
"model": GPy.likelihoods.gaussian(variance=self.var, D=self.D, N=self.N), "model": GPy.likelihoods.Gaussian(variance=self.var),
"grad_params": { "grad_params": {
"names": ["noise_model_variance"], "names": ["variance"],
"vals": [self.var], "vals": [self.var],
"constraints": [constrain_positive] "constraints": [("variance", constrain_positive)]
}, },
"laplace": True, "laplace": True,
"ep": True "ep": True
}, },
#"Gaussian_log": { #"Gaussian_log": {
#"model": GPy.likelihoods.gaussian(gp_link=gp_transformations.Log(), variance=self.var, D=self.D, N=self.N), #"model": GPy.likelihoods.gaussian(gp_link=link_functions.Log(), variance=self.var, D=self.D, N=self.N),
#"grad_params": { #"grad_params": {
#"names": ["noise_model_variance"], #"names": ["noise_model_variance"],
#"vals": [self.var], #"vals": [self.var],
@ -207,7 +219,7 @@ class TestNoiseModels(object):
#"laplace": True #"laplace": True
#}, #},
#"Gaussian_probit": { #"Gaussian_probit": {
#"model": GPy.likelihoods.gaussian(gp_link=gp_transformations.Probit(), variance=self.var, D=self.D, N=self.N), #"model": GPy.likelihoods.gaussian(gp_link=link_functions.Probit(), variance=self.var, D=self.D, N=self.N),
#"grad_params": { #"grad_params": {
#"names": ["noise_model_variance"], #"names": ["noise_model_variance"],
#"vals": [self.var], #"vals": [self.var],
@ -216,7 +228,7 @@ class TestNoiseModels(object):
#"laplace": True #"laplace": True
#}, #},
#"Gaussian_log_ex": { #"Gaussian_log_ex": {
#"model": GPy.likelihoods.gaussian(gp_link=gp_transformations.Log_ex_1(), variance=self.var, D=self.D, N=self.N), #"model": GPy.likelihoods.gaussian(gp_link=link_functions.Log_ex_1(), variance=self.var, D=self.D, N=self.N),
#"grad_params": { #"grad_params": {
#"names": ["noise_model_variance"], #"names": ["noise_model_variance"],
#"vals": [self.var], #"vals": [self.var],
@ -225,31 +237,31 @@ class TestNoiseModels(object):
#"laplace": True #"laplace": True
#}, #},
"Bernoulli_default": { "Bernoulli_default": {
"model": GPy.likelihoods.bernoulli(), "model": GPy.likelihoods.Bernoulli(),
"link_f_constraints": [partial(constrain_bounded, lower=0, upper=1)], "link_f_constraints": [partial(constrain_bounded, lower=0, upper=1)],
"laplace": True, "laplace": True,
"Y": self.binary_Y, "Y": self.binary_Y,
"ep": True "ep": True
}, },
"Exponential_default": { #"Exponential_default": {
"model": GPy.likelihoods.exponential(), #"model": GPy.likelihoods.exponential(),
"link_f_constraints": [constrain_positive], #"link_f_constraints": [constrain_positive],
"Y": self.positive_Y, #"Y": self.positive_Y,
"laplace": True, #"laplace": True,
}, #},
"Poisson_default": { #"Poisson_default": {
"model": GPy.likelihoods.poisson(), #"model": GPy.likelihoods.poisson(),
"link_f_constraints": [constrain_positive], #"link_f_constraints": [constrain_positive],
"Y": self.integer_Y, #"Y": self.integer_Y,
"laplace": True, #"laplace": True,
"ep": False #Should work though... #"ep": False #Should work though...
}, #},
"Gamma_default": { #"Gamma_default": {
"model": GPy.likelihoods.gamma(), #"model": GPy.likelihoods.gamma(),
"link_f_constraints": [constrain_positive], #"link_f_constraints": [constrain_positive],
"Y": self.positive_Y, #"Y": self.positive_Y,
"laplace": True #"laplace": True
} #}
} }
for name, attributes in noise_models.iteritems(): for name, attributes in noise_models.iteritems():
@ -286,8 +298,8 @@ class TestNoiseModels(object):
else: else:
ep = False ep = False
if len(param_vals) > 1: #if len(param_vals) > 1:
raise NotImplementedError("Cannot support multiple params in likelihood yet!") #raise NotImplementedError("Cannot support multiple params in likelihood yet!")
#Required by all #Required by all
#Normal derivatives #Normal derivatives
@ -302,13 +314,13 @@ class TestNoiseModels(object):
yield self.t_d3logpdf_df3, model, Y, f yield self.t_d3logpdf_df3, model, Y, f
yield self.t_d3logpdf_dlink3, model, Y, f, link_f_constraints yield self.t_d3logpdf_dlink3, model, Y, f, link_f_constraints
#Params #Params
yield self.t_dlogpdf_dparams, model, Y, f, param_vals, param_constraints yield self.t_dlogpdf_dparams, model, Y, f, param_vals, param_names, param_constraints
yield self.t_dlogpdf_df_dparams, model, Y, f, param_vals, param_constraints yield self.t_dlogpdf_df_dparams, model, Y, f, param_vals, param_names, param_constraints
yield self.t_d2logpdf2_df2_dparams, model, Y, f, param_vals, param_constraints yield self.t_d2logpdf2_df2_dparams, model, Y, f, param_vals, param_names, param_constraints
#Link params #Link params
yield self.t_dlogpdf_link_dparams, model, Y, f, param_vals, param_constraints yield self.t_dlogpdf_link_dparams, model, Y, f, param_vals, param_names, param_constraints
yield self.t_dlogpdf_dlink_dparams, model, Y, f, param_vals, param_constraints yield self.t_dlogpdf_dlink_dparams, model, Y, f, param_vals, param_names, param_constraints
yield self.t_d2logpdf2_dlink2_dparams, model, Y, f, param_vals, param_constraints yield self.t_d2logpdf2_dlink2_dparams, model, Y, f, param_vals, param_names, param_constraints
#laplace likelihood gradcheck #laplace likelihood gradcheck
yield self.t_laplace_fit_rbf_white, model, self.X, Y, f, self.step, param_vals, param_names, param_constraints yield self.t_laplace_fit_rbf_white, model, self.X, Y, f, self.step, param_vals, param_names, param_constraints
@ -370,33 +382,33 @@ class TestNoiseModels(object):
# df_dparams # # df_dparams #
############## ##############
@with_setup(setUp, tearDown) @with_setup(setUp, tearDown)
def t_dlogpdf_dparams(self, model, Y, f, params, param_constraints): def t_dlogpdf_dparams(self, model, Y, f, params, params_names, param_constraints):
print "\n{}".format(inspect.stack()[0][3]) print "\n{}".format(inspect.stack()[0][3])
print model print model
assert ( assert (
dparam_checkgrad(model.logpdf, model.dlogpdf_dtheta, dparam_checkgrad(model.logpdf, model.dlogpdf_dtheta,
params, args=(f, Y), constraints=param_constraints, params, params_names, args=(f, Y), constraints=param_constraints,
randomize=True, verbose=True) randomize=False, verbose=True)
) )
@with_setup(setUp, tearDown) @with_setup(setUp, tearDown)
def t_dlogpdf_df_dparams(self, model, Y, f, params, param_constraints): def t_dlogpdf_df_dparams(self, model, Y, f, params, params_names, param_constraints):
print "\n{}".format(inspect.stack()[0][3]) print "\n{}".format(inspect.stack()[0][3])
print model print model
assert ( assert (
dparam_checkgrad(model.dlogpdf_df, model.dlogpdf_df_dtheta, dparam_checkgrad(model.dlogpdf_df, model.dlogpdf_df_dtheta,
params, args=(f, Y), constraints=param_constraints, params, params_names, args=(f, Y), constraints=param_constraints,
randomize=True, verbose=True) randomize=False, verbose=True)
) )
@with_setup(setUp, tearDown) @with_setup(setUp, tearDown)
def t_d2logpdf2_df2_dparams(self, model, Y, f, params, param_constraints): def t_d2logpdf2_df2_dparams(self, model, Y, f, params, params_names, param_constraints):
print "\n{}".format(inspect.stack()[0][3]) print "\n{}".format(inspect.stack()[0][3])
print model print model
assert ( assert (
dparam_checkgrad(model.d2logpdf_df2, model.d2logpdf_df2_dtheta, dparam_checkgrad(model.d2logpdf_df2, model.d2logpdf_df2_dtheta,
params, args=(f, Y), constraints=param_constraints, params, params_names, args=(f, Y), constraints=param_constraints,
randomize=True, verbose=True) randomize=False, verbose=True)
) )
################ ################
@ -454,32 +466,32 @@ class TestNoiseModels(object):
# dlink_dparams # # dlink_dparams #
################# #################
@with_setup(setUp, tearDown) @with_setup(setUp, tearDown)
def t_dlogpdf_link_dparams(self, model, Y, f, params, param_constraints): def t_dlogpdf_link_dparams(self, model, Y, f, params, param_names, param_constraints):
print "\n{}".format(inspect.stack()[0][3]) print "\n{}".format(inspect.stack()[0][3])
print model print model
assert ( assert (
dparam_checkgrad(model.logpdf_link, model.dlogpdf_link_dtheta, dparam_checkgrad(model.logpdf_link, model.dlogpdf_link_dtheta,
params, args=(f, Y), constraints=param_constraints, params, param_names, args=(f, Y), constraints=param_constraints,
randomize=False, verbose=True) randomize=False, verbose=True)
) )
@with_setup(setUp, tearDown) @with_setup(setUp, tearDown)
def t_dlogpdf_dlink_dparams(self, model, Y, f, params, param_constraints): def t_dlogpdf_dlink_dparams(self, model, Y, f, params, param_names, param_constraints):
print "\n{}".format(inspect.stack()[0][3]) print "\n{}".format(inspect.stack()[0][3])
print model print model
assert ( assert (
dparam_checkgrad(model.dlogpdf_dlink, model.dlogpdf_dlink_dtheta, dparam_checkgrad(model.dlogpdf_dlink, model.dlogpdf_dlink_dtheta,
params, args=(f, Y), constraints=param_constraints, params, param_names, args=(f, Y), constraints=param_constraints,
randomize=False, verbose=True) randomize=False, verbose=True)
) )
@with_setup(setUp, tearDown) @with_setup(setUp, tearDown)
def t_d2logpdf2_dlink2_dparams(self, model, Y, f, params, param_constraints): def t_d2logpdf2_dlink2_dparams(self, model, Y, f, params, param_names, param_constraints):
print "\n{}".format(inspect.stack()[0][3]) print "\n{}".format(inspect.stack()[0][3])
print model print model
assert ( assert (
dparam_checkgrad(model.d2logpdf_dlink2, model.d2logpdf_dlink2_dtheta, dparam_checkgrad(model.d2logpdf_dlink2, model.d2logpdf_dlink2_dtheta,
params, args=(f, Y), constraints=param_constraints, params, param_names, args=(f, Y), constraints=param_constraints,
randomize=False, verbose=True) randomize=False, verbose=True)
) )
@ -493,18 +505,23 @@ class TestNoiseModels(object):
Y = Y/Y.max() Y = Y/Y.max()
white_var = 1e-6 white_var = 1e-6
kernel = GPy.kern.rbf(X.shape[1]) + GPy.kern.white(X.shape[1]) kernel = GPy.kern.rbf(X.shape[1]) + GPy.kern.white(X.shape[1])
laplace_likelihood = GPy.likelihoods.Laplace(Y.copy(), model) laplace_likelihood = GPy.inference.latent_function_inference.LaplaceInference()
m = GPy.models.GPRegression(X.copy(), Y.copy(), kernel, likelihood=laplace_likelihood) m = GPy.core.GP(X.copy(), Y.copy(), kernel, likelihood=model, inference_method=laplace_likelihood)
m.ensure_default_constraints() m.ensure_default_constraints()
m.constrain_fixed('white', white_var) m['white'].constrain_fixed(white_var)
for param_num in range(len(param_names)): #Set constraints
name = param_names[param_num] for constrain_param, constraint in constraints:
m[name] = param_vals[param_num] constraint(constrain_param, m)
constraints[param_num](name, m)
print m print m
m.randomize() m.randomize()
#Set params
for param_num in range(len(param_names)):
name = param_names[param_num]
m[name] = param_vals[param_num]
#m.optimize(max_iters=8) #m.optimize(max_iters=8)
print m print m
m.checkgrad(verbose=1, step=step) m.checkgrad(verbose=1, step=step)
@ -526,9 +543,9 @@ class TestNoiseModels(object):
white_var = 1e-6 white_var = 1e-6
kernel = GPy.kern.rbf(X.shape[1]) + GPy.kern.white(X.shape[1]) kernel = GPy.kern.rbf(X.shape[1]) + GPy.kern.white(X.shape[1])
ep_likelihood = GPy.likelihoods.EP(Y.copy(), model) ep_likelihood = GPy.likelihoods.EP(Y.copy(), model)
m = GPy.models.GPRegression(X.copy(), Y.copy(), kernel, likelihood=ep_likelihood) m = GPy.core.GP(X.copy(), Y.copy(), kernel, likelihood=ep_likelihood)
m.ensure_default_constraints() m.ensure_default_constraints()
m.constrain_fixed('white', white_var) m['white'].constrain_fixed(white_var)
for param_num in range(len(param_names)): for param_num in range(len(param_names)):
name = param_names[param_num] name = param_names[param_num]
@ -559,8 +576,8 @@ class LaplaceTests(unittest.TestCase):
self.var = 0.2 self.var = 0.2
self.var = np.random.rand(1) self.var = np.random.rand(1)
self.stu_t = GPy.likelihoods.student_t(deg_free=5, sigma2=self.var) self.stu_t = GPy.likelihoods.StudentT(deg_free=5, sigma2=self.var)
self.gauss = GPy.likelihoods.gaussian(gp_transformations.Log(), variance=self.var, D=self.D, N=self.N) self.gauss = GPy.likelihoods.Gaussian(gp_link=link_functions.Log(), variance=self.var)
#Make a bigger step as lower bound can be quite curved #Make a bigger step as lower bound can be quite curved
self.step = 1e-6 self.step = 1e-6
@ -584,7 +601,7 @@ class LaplaceTests(unittest.TestCase):
noise = np.random.randn(*self.X.shape)*self.real_std noise = np.random.randn(*self.X.shape)*self.real_std
self.Y = np.sin(self.X*2*np.pi) + noise self.Y = np.sin(self.X*2*np.pi) + noise
self.f = np.random.rand(self.N, 1) self.f = np.random.rand(self.N, 1)
self.gauss = GPy.likelihoods.gaussian(variance=self.var, D=self.D, N=self.N) self.gauss = GPy.likelihoods.Gaussian(variance=self.var)
dlogpdf_df = functools.partial(self.gauss.dlogpdf_df, y=self.Y) dlogpdf_df = functools.partial(self.gauss.dlogpdf_df, y=self.Y)
d2logpdf_df2 = functools.partial(self.gauss.d2logpdf_df2, y=self.Y) d2logpdf_df2 = functools.partial(self.gauss.d2logpdf_df2, y=self.Y)
@ -607,21 +624,23 @@ class LaplaceTests(unittest.TestCase):
kernel1 = GPy.kern.rbf(X.shape[1]) + GPy.kern.white(X.shape[1]) kernel1 = GPy.kern.rbf(X.shape[1]) + GPy.kern.white(X.shape[1])
kernel2 = kernel1.copy() kernel2 = kernel1.copy()
m1 = GPy.models.GPRegression(X, Y.copy(), kernel=kernel1) gauss_distr1 = GPy.likelihoods.Gaussian(variance=initial_var_guess)
m1.constrain_fixed('white', 1e-6) exact_inf = GPy.inference.latent_function_inference.ExactGaussianInference()
m1['noise'] = initial_var_guess m1 = GPy.core.GP(X, Y.copy(), kernel=kernel1, likelihood=gauss_distr1, inference_method=exact_inf)
m1.constrain_bounded('noise', 1e-4, 10) m1['white'].constrain_fixed(1e-6)
m1.constrain_bounded('rbf', 1e-4, 10) m1['variance'] = initial_var_guess
m1['variance'].constrain_bounded(1e-4, 10)
m1['rbf'].constrain_bounded(1e-4, 10)
m1.ensure_default_constraints() m1.ensure_default_constraints()
m1.randomize() m1.randomize()
gauss_distr = GPy.likelihoods.gaussian(variance=initial_var_guess, D=1, N=Y.shape[0]) gauss_distr2 = GPy.likelihoods.Gaussian(variance=initial_var_guess)
laplace_likelihood = GPy.likelihoods.Laplace(Y.copy(), gauss_distr) laplace_inf = GPy.inference.latent_function_inference.LaplaceInference()
m2 = GPy.models.GPRegression(X, Y.copy(), kernel=kernel2, likelihood=laplace_likelihood) m2 = GPy.core.GP(X, Y.copy(), kernel=kernel2, likelihood=gauss_distr2, inference_method=laplace_inf)
m2.ensure_default_constraints() m2.ensure_default_constraints()
m2.constrain_fixed('white', 1e-6) m2['white'].constrain_fixed(1e-6)
m2.constrain_bounded('rbf', 1e-4, 10) m2['rbf'].constrain_bounded(1e-4, 10)
m2.constrain_bounded('noise', 1e-4, 10) m2['variance'].constrain_bounded(1e-4, 10)
m2.randomize() m2.randomize()
if debug: if debug: