diff --git a/GPy/inference/latent_function_inference/laplace.py b/GPy/inference/latent_function_inference/laplace.py index 6e252406..6fce94fc 100644 --- a/GPy/inference/latent_function_inference/laplace.py +++ b/GPy/inference/latent_function_inference/laplace.py @@ -33,7 +33,6 @@ class LaplaceInference(object): self._mode_finding_max_iter = 40 self.bad_fhat = True - def inference(self, kern, X, likelihood, Y, Y_metadata=None): """ Returns a Posterior class containing essential quantities of the posterior @@ -50,6 +49,7 @@ class LaplaceInference(object): Ki_f_init = np.zeros_like(Y) else: Ki_f_init = self._previous_Ki_fhat + f_hat, Ki_fhat = self.rasm_mode(K, Y, likelihood, Ki_f_init, Y_metadata=Y_metadata) #Compute hessian and other variables at mode @@ -57,6 +57,7 @@ class LaplaceInference(object): #likelihood.gradient = self.likelihood_gradients() kern.update_gradients_full(dL_dK, X) + likelihood.update_gradients(np.ones(10)) 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} @@ -157,9 +158,12 @@ class LaplaceInference(object): explicit_part = 0.5*(np.dot(Ki_f, Ki_f.T) - K_Wi_i) #Implicit - d3lik_d3fhat = likelihood.d3logpdf_df3(f_hat, Y, extra_data=Y_metadata) - 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. + dW_df = likelihood.d3logpdf_df3(f_hat, Y, extra_data=Y_metadata) # d3lik_d3fhat 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)) dL_dK = explicit_part + implicit_part @@ -219,7 +223,7 @@ class LaplaceInference(object): LiW12, _ = dtrtrs(L, np.diagflat(W_12), lower=1, trans=0) K_Wi_i = np.dot(LiW12.T, LiW12) # R = W12BiW12, in R&W p 126, eq 5.25 - #here's a better way to compute the required matrix. + #here's a better way to compute the required matrix. # you could do the model finding witha backsub, instead of a dot... #L2 = L/W_12 #K_Wi_i_2 , _= dpotri(L2) diff --git a/GPy/likelihoods/gaussian.py b/GPy/likelihoods/gaussian.py index c047e573..e6be2261 100644 --- a/GPy/likelihoods/gaussian.py +++ b/GPy/likelihoods/gaussian.py @@ -3,9 +3,9 @@ #TODO """ -A lot of this code assumes that the link functio nis the identity. +A lot of this code assumes that the link functio nis the identity. -I think laplace code is okay, but I'm quite sure that the EP moments will only work if the link is identity. +I think laplace code is okay, but I'm quite sure that the EP moments will only work if the link is identity. Furthermore, exact Guassian inference can only be done for the identity link, so we should be asserting so for all calls which relate to that. @@ -130,7 +130,10 @@ class Gaussian(Likelihood): :rtype: float """ 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): """ @@ -175,7 +178,8 @@ class Gaussian(Likelihood): (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 - 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 def d3logpdf_dlink3(self, link_f, y, extra_data=None): @@ -194,7 +198,8 @@ class Gaussian(Likelihood): :rtype: Nx1 array """ 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 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 e = y - link_f 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? 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 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 def dlogpdf_link_dtheta(self, f, y, extra_data=None): diff --git a/GPy/likelihoods/likelihood.py b/GPy/likelihoods/likelihood.py index b0ecfc37..d9a7e109 100644 --- a/GPy/likelihoods/likelihood.py +++ b/GPy/likelihoods/likelihood.py @@ -13,12 +13,12 @@ from ..core.parameterization import Parameterized class Likelihood(Parameterized): """ - Likelihood base class, used to defing p(y|f). + Likelihood base class, used to defing p(y|f). All instances use _inverse_ link functions, which can be swapped out. It is expected that inherriting classes define a default inverse link function - To use this class, inherrit and define missing functionality. + To use this class, inherrit and define missing functionality. Inherriting classes *must* implement: pdf_link : a bound method which turns the output of the link function into the pdf @@ -27,7 +27,7 @@ class Likelihood(Parameterized): To enable use with EP, inherriting classes *must* define: TODO: a suitable derivative function for any parameters of the class It is also desirable to define: - moments_match_ep : a function to compute the EP moments If this isn't defined, the moments will be computed using 1D quadrature. + moments_match_ep : a function to compute the EP moments If this isn't defined, the moments will be computed using 1D quadrature. To enable use with Laplace approximation, inherriting classes *must* define: Some derivative functions *AS TODO* @@ -36,7 +36,7 @@ class Likelihood(Parameterized): """ def __init__(self, gp_link, name): - super(Likelihood, self).__init__(name) + super(Likelihood, self).__init__(name) assert isinstance(gp_link,link_functions.GPTransformation), "gp_link is not a valid GPTransformation." self.gp_link = gp_link self.log_concave = False @@ -44,6 +44,10 @@ class Likelihood(Parameterized): def _gradients(self,partial): 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): """ 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 """ - if len(self._get_param_names()) > 0: + if self.size > 0: link_f = self.gp_link.transf(f) return self.dlogpdf_link_dtheta(link_f, y, extra_data=extra_data) else: @@ -314,7 +318,7 @@ class Likelihood(Parameterized): """ TODO: Doc strings """ - if len(self._get_param_names()) > 0: + if self.size > 0: link_f = self.gp_link.transf(f) dlink_df = self.gp_link.dtransf_df(f) dlogpdf_dlink_dtheta = self.dlogpdf_dlink_dtheta(link_f, y, extra_data=extra_data) @@ -327,7 +331,7 @@ class Likelihood(Parameterized): """ TODO: Doc strings """ - if len(self._get_param_names()) > 0: + if self.size > 0: link_f = self.gp_link.transf(f) dlink_df = self.gp_link.dtransf_df(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' # ensure we have gradients for every parameter we want to optimize - assert dlogpdf_dtheta.shape[1] == len(self._get_param_names()) - assert dlogpdf_df_dtheta.shape[1] == len(self._get_param_names()) - assert d2logpdf_df2_dtheta.shape[1] == len(self._get_param_names()) + assert dlogpdf_dtheta.shape[1] == self.size + assert dlogpdf_df_dtheta.shape[1] == self.size + assert d2logpdf_df2_dtheta.shape[1] == self.size return dlogpdf_dtheta, dlogpdf_df_dtheta, d2logpdf_df2_dtheta diff --git a/GPy/likelihoods/student_t.py b/GPy/likelihoods/student_t.py index 587e1b23..b4e0dfc3 100644 --- a/GPy/likelihoods/student_t.py +++ b/GPy/likelihoods/student_t.py @@ -8,6 +8,7 @@ import link_functions from scipy import stats, integrate from scipy.special import gammaln, gamma from likelihood import Likelihood +from ..core.parameterization import Param 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}} """ - def __init__(self,gp_link=None,analytical_mean=True,analytical_variance=True, deg_free=5, sigma2=2): - self.v = deg_free - self.sigma2 = sigma2 + def __init__(self,gp_link=None, deg_free=5, sigma2=2): + if gp_link is None: + 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 - def _get_params(self): - return np.asarray(self.sigma2) + def parameters_changed(self): + self.variance = (self.v / float(self.v - 2)) * self.sigma2 - def _get_param_names(self): - return ["t_noise_std2"] - - 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 update_gradients(self, partial): + self.sigma2.gradient = np.ones(1) #FIXME: Not done yet + self.v.gradient = np.ones(1) #FIXME: Not done yet def pdf_link(self, link_f, y, extra_data=None): """ @@ -82,10 +82,14 @@ class StudentT(Likelihood): """ assert np.atleast_1d(link_f).shape == np.atleast_1d(y).shape 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) - - gammaln(self.v * 0.5) - - 0.5*np.log(self.sigma2 * self.v * np.pi) - - 0.5*(self.v + 1)*np.log(1 + (1/np.float(self.v))*((e**2)/self.sigma2)) + - gammaln(self.v * 0.5) + - 0.5*np.log(self.sigma2 * self.v * np.pi) + - 0.5*(self.v + 1)*np.log(1 + (1/np.float(self.v))*((e**2)/self.sigma2)) ) return np.sum(objective) diff --git a/GPy/models/gradient_checker.py b/GPy/models/gradient_checker.py index 775334ac..b7c78449 100644 --- a/GPy/models/gradient_checker.py +++ b/GPy/models/gradient_checker.py @@ -1,9 +1,10 @@ # ## Copyright (c) 2012, GPy authors (see AUTHORS.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 numpy +from ..core.parameterization import Param def get_shape(x): if isinstance(x, numpy.ndarray): @@ -24,42 +25,42 @@ class GradientChecker(Model): """ :param f: Function to check gradient for :param df: Gradient of function to check - :param x0: + :param x0: Initial guess for inputs x (if it has a shape (a,b) this will be reflected in the parameter names). - Can be a list of arrays, if takes a list of arrays. This list will be passed + Can be a list of arrays, if takes a list of arrays. This list will be passed to f and df in the same order as given here. If only one argument, make sure not to pass a list!!! - + :type x0: [array-like] | array-like | float | int :param names: Names to print, when performing gradcheck. If a list was passed to x0 a list of names with the same length is expected. :param args: Arguments passed as f(x, *args, **kwargs) and df(x, *args, **kwargs) - + Examples: --------- from GPy.models import GradientChecker N, M, Q = 10, 5, 3 - + Sinusoid: - + X = numpy.random.rand(N, Q) grad = GradientChecker(numpy.sin,numpy.cos,X,'x') grad.checkgrad(verbose=1) - + Using GPy: - + X, Z = numpy.random.randn(N,Q), numpy.random.randn(M,Q) kern = GPy.kern.linear(Q, ARD=True) + GPy.kern.rbf(Q, ARD=True) - grad = GradientChecker(kern.K, + grad = GradientChecker(kern.K, lambda x: 2*kern.dK_dX(numpy.ones((1,1)), x), x0 = X.copy(), - names='X') + names='X') grad.checkgrad(verbose=1) 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: self.shapes = [get_shape(xi) for xi in x0] self.names = ['X{i}'.format(i=i) for i in range(len(x0))] @@ -72,8 +73,10 @@ class GradientChecker(Model): else: self.names = names self.shapes = [get_shape(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 = [] # 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)))))) @@ -93,20 +96,18 @@ class GradientChecker(Model): def _log_likelihood_gradients(self): 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): - return numpy.atleast_1d(numpy.hstack(map(lambda name: flatten_if_needed(self.__getattribute__(name)), self.names))) + #def _set_params(self, x): + #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 _set_params(self, x): - 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): - _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 + #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 diff --git a/GPy/testing/likelihood_tests.py b/GPy/testing/likelihood_tests.py index d14c9a41..c418a096 100644 --- a/GPy/testing/likelihood_tests.py +++ b/GPy/testing/likelihood_tests.py @@ -4,7 +4,8 @@ import GPy from GPy.models import GradientChecker import functools 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 #np.random.seed(300) 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 param """ - def param_func(param, inst_func, args): - inst_func.im_self._set_params(param) + def param_func(param_val, param_name, inst_func, args): + #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 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 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_df = dparam_partial(dfunc, *args) gradchecking = True - for param in params: - fnum = np.atleast_1d(partial_f(param)).shape[0] - dfnum = np.atleast_1d(partial_df(param)).shape[0] + zipped_params = zip(params, params_names) + for param_val, param_name in zipped_params: + 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): #dlik and dlik_dvar gives back 1 value for each f_ind = min(fnum, fixed_val+1) - 1 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 #The parameter is being set directly with __setattr__ - grad = GradientChecker(lambda x: np.atleast_1d(partial_f(x))[f_ind], - lambda x : np.atleast_1d(partial_df(x))[fixed_val], - param, 'p') + grad = GradientChecker(lambda p_val: np.atleast_1d(partial_f(p_val, param_name))[f_ind], + lambda p_val: np.atleast_1d(partial_df(p_val, param_name))[fixed_val], + param_val, [param_name]) #This is not general for more than one param... if constraints is not None: - for constraint in constraints: - constraint('p', grad) + for constrain_param, constraint in constraints: + if grad.grep_param_names(constrain_param): + constraint(constrain_param, grad) + else: + print "parameter didn't exist" + print constrain_param, " ", constraint if randomize: grad.randomize() if verbose: @@ -107,17 +115,20 @@ class TestNoiseModels(object): #################################################### # 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): - model.constrain_negative(regex) + model[regex].constrain_negative() def constrain_positive(regex, model): - model.constrain_positive(regex) + model[regex].constrain_positive() def constrain_bounded(regex, model, lower, upper): """ 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 @@ -134,71 +145,72 @@ class TestNoiseModels(object): } """ 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": { "names": ["t_noise"], "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 }, "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": { "names": ["t_noise"], "vals": [1.0], - "constraints": [constrain_positive] + "constraints": [("t_noise", constrain_positive), ("deg_free", constrain_positive)] }, "laplace": True }, "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": { "names": ["t_noise"], "vals": [0.01], - "constraints": [constrain_positive] + "constraints": [("t_noise", constrain_positive), ("deg_free", constrain_positive)] }, "laplace": True }, "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": { "names": ["t_noise"], "vals": [10.0], - "constraints": [constrain_positive] + "constraints": [("t_noise", constrain_positive), ("deg_free", constrain_positive)] }, "laplace": True }, "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": { "names": ["t_noise"], "vals": [self.var], - "constraints": [constrain_positive] + "constraints": [("t_noise", constrain_positive), ("deg_free", constrain_positive)] }, "laplace": True }, "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": { "names": ["t_noise"], "vals": [self.var], - "constraints": [constrain_positive] + "constraints": [("t_noise", constrain_positive), ("deg_free", constrain_positive)] }, "laplace": True }, "Gaussian_default": { - "model": GPy.likelihoods.gaussian(variance=self.var, D=self.D, N=self.N), + "model": GPy.likelihoods.Gaussian(variance=self.var), "grad_params": { - "names": ["noise_model_variance"], + "names": ["variance"], "vals": [self.var], - "constraints": [constrain_positive] + "constraints": [("variance", constrain_positive)] }, "laplace": True, "ep": True }, #"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": { #"names": ["noise_model_variance"], #"vals": [self.var], @@ -207,7 +219,7 @@ class TestNoiseModels(object): #"laplace": True #}, #"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": { #"names": ["noise_model_variance"], #"vals": [self.var], @@ -216,7 +228,7 @@ class TestNoiseModels(object): #"laplace": True #}, #"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": { #"names": ["noise_model_variance"], #"vals": [self.var], @@ -225,31 +237,31 @@ class TestNoiseModels(object): #"laplace": True #}, "Bernoulli_default": { - "model": GPy.likelihoods.bernoulli(), + "model": GPy.likelihoods.Bernoulli(), "link_f_constraints": [partial(constrain_bounded, lower=0, upper=1)], "laplace": True, "Y": self.binary_Y, "ep": True }, - "Exponential_default": { - "model": GPy.likelihoods.exponential(), - "link_f_constraints": [constrain_positive], - "Y": self.positive_Y, - "laplace": True, - }, - "Poisson_default": { - "model": GPy.likelihoods.poisson(), - "link_f_constraints": [constrain_positive], - "Y": self.integer_Y, - "laplace": True, - "ep": False #Should work though... - }, - "Gamma_default": { - "model": GPy.likelihoods.gamma(), - "link_f_constraints": [constrain_positive], - "Y": self.positive_Y, - "laplace": True - } + #"Exponential_default": { + #"model": GPy.likelihoods.exponential(), + #"link_f_constraints": [constrain_positive], + #"Y": self.positive_Y, + #"laplace": True, + #}, + #"Poisson_default": { + #"model": GPy.likelihoods.poisson(), + #"link_f_constraints": [constrain_positive], + #"Y": self.integer_Y, + #"laplace": True, + #"ep": False #Should work though... + #}, + #"Gamma_default": { + #"model": GPy.likelihoods.gamma(), + #"link_f_constraints": [constrain_positive], + #"Y": self.positive_Y, + #"laplace": True + #} } for name, attributes in noise_models.iteritems(): @@ -286,8 +298,8 @@ class TestNoiseModels(object): else: ep = False - if len(param_vals) > 1: - raise NotImplementedError("Cannot support multiple params in likelihood yet!") + #if len(param_vals) > 1: + #raise NotImplementedError("Cannot support multiple params in likelihood yet!") #Required by all #Normal derivatives @@ -302,13 +314,13 @@ class TestNoiseModels(object): yield self.t_d3logpdf_df3, model, Y, f yield self.t_d3logpdf_dlink3, model, Y, f, link_f_constraints #Params - yield self.t_dlogpdf_dparams, model, Y, f, param_vals, param_constraints - yield self.t_dlogpdf_df_dparams, model, Y, f, param_vals, param_constraints - yield self.t_d2logpdf2_df2_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_names, param_constraints + yield self.t_d2logpdf2_df2_dparams, model, Y, f, param_vals, param_names, param_constraints #Link params - yield self.t_dlogpdf_link_dparams, model, Y, f, param_vals, param_constraints - yield self.t_dlogpdf_dlink_dparams, model, Y, f, param_vals, param_constraints - yield self.t_d2logpdf2_dlink2_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_names, param_constraints + yield self.t_d2logpdf2_dlink2_dparams, model, Y, f, param_vals, param_names, param_constraints #laplace likelihood gradcheck 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 # ############## @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 model assert ( dparam_checkgrad(model.logpdf, model.dlogpdf_dtheta, - params, args=(f, Y), constraints=param_constraints, - randomize=True, verbose=True) + params, params_names, args=(f, Y), constraints=param_constraints, + randomize=False, verbose=True) ) @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 model assert ( dparam_checkgrad(model.dlogpdf_df, model.dlogpdf_df_dtheta, - params, args=(f, Y), constraints=param_constraints, - randomize=True, verbose=True) + params, params_names, args=(f, Y), constraints=param_constraints, + randomize=False, verbose=True) ) @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 model assert ( dparam_checkgrad(model.d2logpdf_df2, model.d2logpdf_df2_dtheta, - params, args=(f, Y), constraints=param_constraints, - randomize=True, verbose=True) + params, params_names, args=(f, Y), constraints=param_constraints, + randomize=False, verbose=True) ) ################ @@ -454,32 +466,32 @@ class TestNoiseModels(object): # dlink_dparams # ################# @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 model assert ( 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) ) @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 model assert ( 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) ) @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 model assert ( 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) ) @@ -493,18 +505,23 @@ class TestNoiseModels(object): Y = Y/Y.max() white_var = 1e-6 kernel = GPy.kern.rbf(X.shape[1]) + GPy.kern.white(X.shape[1]) - laplace_likelihood = GPy.likelihoods.Laplace(Y.copy(), model) - m = GPy.models.GPRegression(X.copy(), Y.copy(), kernel, likelihood=laplace_likelihood) + laplace_likelihood = GPy.inference.latent_function_inference.LaplaceInference() + m = GPy.core.GP(X.copy(), Y.copy(), kernel, likelihood=model, inference_method=laplace_likelihood) m.ensure_default_constraints() - m.constrain_fixed('white', white_var) + m['white'].constrain_fixed(white_var) - for param_num in range(len(param_names)): - name = param_names[param_num] - m[name] = param_vals[param_num] - constraints[param_num](name, m) + #Set constraints + for constrain_param, constraint in constraints: + constraint(constrain_param, m) print m 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) print m m.checkgrad(verbose=1, step=step) @@ -526,9 +543,9 @@ class TestNoiseModels(object): white_var = 1e-6 kernel = GPy.kern.rbf(X.shape[1]) + GPy.kern.white(X.shape[1]) 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.constrain_fixed('white', white_var) + m['white'].constrain_fixed(white_var) for param_num in range(len(param_names)): name = param_names[param_num] @@ -559,8 +576,8 @@ class LaplaceTests(unittest.TestCase): self.var = 0.2 self.var = np.random.rand(1) - self.stu_t = GPy.likelihoods.student_t(deg_free=5, sigma2=self.var) - self.gauss = GPy.likelihoods.gaussian(gp_transformations.Log(), variance=self.var, D=self.D, N=self.N) + self.stu_t = GPy.likelihoods.StudentT(deg_free=5, sigma2=self.var) + self.gauss = GPy.likelihoods.Gaussian(gp_link=link_functions.Log(), variance=self.var) #Make a bigger step as lower bound can be quite curved self.step = 1e-6 @@ -584,7 +601,7 @@ class LaplaceTests(unittest.TestCase): noise = np.random.randn(*self.X.shape)*self.real_std self.Y = np.sin(self.X*2*np.pi) + noise 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) 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]) kernel2 = kernel1.copy() - m1 = GPy.models.GPRegression(X, Y.copy(), kernel=kernel1) - m1.constrain_fixed('white', 1e-6) - m1['noise'] = initial_var_guess - m1.constrain_bounded('noise', 1e-4, 10) - m1.constrain_bounded('rbf', 1e-4, 10) + gauss_distr1 = GPy.likelihoods.Gaussian(variance=initial_var_guess) + exact_inf = GPy.inference.latent_function_inference.ExactGaussianInference() + m1 = GPy.core.GP(X, Y.copy(), kernel=kernel1, likelihood=gauss_distr1, inference_method=exact_inf) + m1['white'].constrain_fixed(1e-6) + m1['variance'] = initial_var_guess + m1['variance'].constrain_bounded(1e-4, 10) + m1['rbf'].constrain_bounded(1e-4, 10) m1.ensure_default_constraints() m1.randomize() - gauss_distr = GPy.likelihoods.gaussian(variance=initial_var_guess, D=1, N=Y.shape[0]) - laplace_likelihood = GPy.likelihoods.Laplace(Y.copy(), gauss_distr) - m2 = GPy.models.GPRegression(X, Y.copy(), kernel=kernel2, likelihood=laplace_likelihood) + gauss_distr2 = GPy.likelihoods.Gaussian(variance=initial_var_guess) + laplace_inf = GPy.inference.latent_function_inference.LaplaceInference() + m2 = GPy.core.GP(X, Y.copy(), kernel=kernel2, likelihood=gauss_distr2, inference_method=laplace_inf) m2.ensure_default_constraints() - m2.constrain_fixed('white', 1e-6) - m2.constrain_bounded('rbf', 1e-4, 10) - m2.constrain_bounded('noise', 1e-4, 10) + m2['white'].constrain_fixed(1e-6) + m2['rbf'].constrain_bounded(1e-4, 10) + m2['variance'].constrain_bounded(1e-4, 10) m2.randomize() if debug: