From 9b2bc907e4ba0a892e9d5236cc228d1370966933 Mon Sep 17 00:00:00 2001 From: Neil Lawrence Date: Mon, 31 Mar 2014 09:20:03 +0100 Subject: [PATCH] Adding of symbolic likelihoods (not yet fully funcitonal). --- GPy/likelihoods/.DS_Store | Bin 12292 -> 0 bytes GPy/likelihoods/symbolic.py | 243 ++++++++++++++++++++++++++++++++++++ 2 files changed, 243 insertions(+) delete mode 100644 GPy/likelihoods/.DS_Store create mode 100644 GPy/likelihoods/symbolic.py diff --git a/GPy/likelihoods/.DS_Store b/GPy/likelihoods/.DS_Store deleted file mode 100644 index 8228ae90a5144c14ccc985f4a2c9f9cfc9145102..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 12292 zcmeHM!EVz)5S?ucggR8T2NVe;eB^}0A4F9_;tCgTCU&4$BzDxqQn=+a_$WSvzX)$; zw{*Q;ClUvQ!fdoVjy-Q?=goSZafrxVHq%!`BO)50vwWOjh%~;=z0;1|ndguOc%qV4 zl;QJ|79DLpp$sSk%78MU3@8Kt0t0wva~t$L_w7|3l>ueozhr>-hZLRV*vXBRvUOmP z5&-f9r={T-b%0<}C&x~1tQ070+UmhLsBtBRv2a{ZSsilhn}h-er`Q4!#NbZ|L+ zq8oUtsi7G@rMxc){e`%Ma%eJY4mkmi-!kk&9a8}>t6)wX%Xcr9dzjH8YELww{b0Mr zNOb>BGmO-{!*$?mv3P%-iyE*o~4d@!a=J-r{BlPtX(QgpXQoK}ssMZXJaa4y3tQ)iTB0KNk z_GElCYhf)&eSal0xa7@vnEh6&k| zx%#L_K$Y|rK3IOiEbIPt+~WPBABlX|bOU?WkvaeLa(A}FMzGJw%zbR!e}-LTa(5>M z{w4FC|8wzQA^wa$V;`yHjw&(x3N~W<4_om0djied%w}MYbG(w+bKaNmeYY(?qNgoe zb0DXZmkhJAR@ibQb7X$wGvu}8TYdItR>PG;uzU2u<8IhvcNuUVuI*mjbXm;ux}Nbb zz*8)KOY2fXkEEirAcjJHBc$;Lm!U^+eEYWokv-rLhhrPmi0l8?{qO$|F?5Yd8Bhis z1LoxF{naJDxZ1iy4A+0^7kKxI8LF;n6Z-^D{%(xufGV; s=f6Jx@zX?~|8o23pbRJj%78MU3@8K2fHI&ACjCZU6uP diff --git a/GPy/likelihoods/symbolic.py b/GPy/likelihoods/symbolic.py new file mode 100644 index 00000000..5eaafb2a --- /dev/null +++ b/GPy/likelihoods/symbolic.py @@ -0,0 +1,243 @@ +# Copyright (c) 2014 GPy Authors +# Licensed under the BSD 3-clause license (see LICENSE.txt) + +import numpy as np +import sympy as sp +from sympy.utilities.lambdify import lambdify +import link_functions +from scipy import stats, integrate +from scipy.special import gammaln, gamma, erf +from likelihood import Likelihood +from ..core.parameterization import Param +from ..core.parameterization.transformations import Logexp + +class Symbolic(Likelihood): + """ + Symbolic likelihood. + + Likelihood where the form of the likelihood is provided by a sympy expression. + + """ + def __init__(self, likelihood=None, log_likelihood=None, cdf=None, logZ=None, gp_link=None, name='symbolic', log_concave=False, param=None): + if gp_link is None: + gp_link = link_functions.Identity() + + if likelihood is None and log_likelihood is None and cdf is None: + raise ValueError, "You must provide an argument for the likelihood or the log likelihood." + + super(Symbolic, self).__init__(gp_link, name=name) + + if likelihood is None and log_likelihood: + self._sp_likelihood = sp.exp(log_likelihood).simplify() + self._sp_log_likelihood = log_likelihood + + if log_likelihood is None and likelihood: + self._sp_likelihood = likelihood + self._sp_log_likelihood = sp.log(likelihood).simplify() + + # TODO: build likelihood and log likelihood from CDF or + # compute CDF given likelihood/log-likelihood. Also check log + # likelihood, likelihood and CDF are consistent. + + # pull the variable names out of the symbolic likelihood + sp_vars = [e for e in self._sp_likelihood.atoms() if e.is_Symbol] + self._sp_f = [e for e in sp_vars if e.name=='f'] + if not self._sp_f: + raise ValueError('No variable f in likelihood or log likelihood.') + self._sp_y = [e for e in sp_vars if e.name=='y'] + if not self._sp_f: + raise ValueError('No variable y in likelihood or log likelihood.') + self._sp_theta = sorted([e for e in sp_vars if not (e.name=='f' or e.name=='y')],key=lambda e:e.name) + + # These are all the arguments need to compute likelihoods. + self.arg_list = self._sp_y + self._sp_f + self._sp_theta + + # these are arguments for computing derivatives. + derivative_arguments = self._sp_f + self._sp_theta + + # Do symbolic work to compute derivatives. + self._log_likelihood_derivatives = {theta.name : sp.diff(self._sp_log_likelihood,theta).simplify() for theta in derivative_arguments} + self._log_likelihood_second_derivatives = {theta.name : sp.diff(self._log_likelihood_derivatives['f'],theta).simplify() for theta in derivative_arguments} + self._log_likelihood_third_derivatives = {theta.name : sp.diff(self._log_likelihood_second_derivatives['f'],theta).simplify() for theta in derivative_arguments} + + # Add parameters to the model. + for theta in self._sp_theta: + val = 1.0 + # TODO: need to decide how to handle user passing values for the se parameter vectors. + if param is not None: + if param.has_key(theta): + val = param[theta] + setattr(self, theta.name, Param(theta.name, val, None)) + self.add_parameters(getattr(self, theta.name)) + + + # Is there some way to check whether the likelihood is log + # concave? For the moment, need user to specify. + self.log_concave = log_concave + + # initialise code arguments + self._arguments = {} + + # generate the code for the likelihood and derivatives + self._gen_code() + + def _gen_code(self): + """Generate the code from the symbolic parts that will be used for likleihod computation.""" + # TODO: Check here whether theano is available and set up + # functions accordingly. + self._likelihood_function = lambdify(self.arg_list, self._sp_likelihood, 'numpy') + self._log_likelihood_function = lambdify(self.arg_list, self._sp_log_likelihood, 'numpy') + + # compute code for derivatives (for implicit likelihood terms + # we need up to 3rd derivatives) + setattr(self, '_first_derivative_code', {key: lambdify(self.arg_list, self._log_likelihood_derivatives[key], 'numpy') for key in self._log_likelihood_derivatives.keys()}) + setattr(self, '_second_derivative_code', {key: lambdify(self.arg_list, self._log_likelihood_second_derivatives[key], 'numpy') for key in self._log_likelihood_second_derivatives.keys()}) + setattr(self, '_third_derivative_code', {key: lambdify(self.arg_list, self._log_likelihood_third_derivatives[key], 'numpy') for key in self._log_likelihood_third_derivatives.keys()}) + + # TODO: compute EP code parts based on logZ. We need dlogZ/dmu, d2logZ/dmu2 and dlogZ/dtheta + + def parameters_changed(self): + pass + + def update_gradients(self, grads): + """ + Pull out the gradients, be careful as the order must match the order + in which the parameters are added + """ + # The way the Laplace approximation is run requires the + # covariance function to compute the true gradient (because it + # is dependent on the mode). This means we actually compute + # the gradient outside this object. This function would + # normally ask the object to update its gradients internally, + # but here it provides them externally, because they are + # computed in the inference code. TODO: Thought: How does this + # effect EP? Shouldn't this be done by a separate + # Laplace-approximation specific call? + for grad, theta in zip(grads, self._sp_theta): + parameter = getattr(self, theta.name) + setattr(parameter, 'gradient', grad) + + def _arguments_update(self, f, y): + """Set up argument lists for the derivatives.""" + # If we do make use of Theano, then at this point we would + # need to do a lot of precomputation to ensure that the + # likelihoods and gradients are computed together, then check + # for parameter changes before updating. + for i, fvar in enumerate(self._sp_f): + self._arguments[fvar.name] = f + for i, yvar in enumerate(self._sp_y): + self._arguments[yvar.name] = y + for theta in self._sp_theta: + self._arguments[theta.name] = np.asarray(getattr(self, theta.name)) + + def pdf_link(self, inv_link_f, y, Y_metadata=None): + """ + Likelihood function given inverse link of f. + + :param inv_link_f: inverse link of latent variables. + :type inv_link_f: Nx1 array + :param y: data + :type y: Nx1 array + :param Y_metadata: Y_metadata which is not used in student t distribution + :returns: likelihood evaluated for this point + :rtype: float + """ + assert np.atleast_1d(inv_link_f).shape == np.atleast_1d(y).shape + self._arguments_update(inv_link_f, y) + l = self._likelihood_function(**self._arguments) + return np.prod(l) + + def logpdf_link(self, inv_link_f, y, Y_metadata=None): + """ + Log Likelihood Function given inverse link of latent variables. + + :param inv_inv_link_f: latent variables (inverse link of f) + :type inv_inv_link_f: Nx1 array + :param y: data + :type y: Nx1 array + :param Y_metadata: Y_metadata + :returns: likelihood evaluated for this point + :rtype: float + + """ + assert np.atleast_1d(inv_link_f).shape == np.atleast_1d(y).shape + self._arguments_update(inv_link_f, y) + ll = self._log_likelihood_function(**self._arguments) + return np.sum(ll) + + def dlogpdf_dlink(self, inv_link_f, y, Y_metadata=None): + """ + Gradient of log likelihood with respect to the inverse link function. + + :param inv_inv_link_f: latent variables (inverse link of f) + :type inv_inv_link_f: Nx1 array + :param y: data + :type y: Nx1 array + :param Y_metadata: Y_metadata + :returns: gradient of likelihood with respect to each point. + :rtype: Nx1 array + + """ + assert np.atleast_1d(inv_link_f).shape == np.atleast_1d(y).shape + self._arguments_update(inv_link_f, y) + return self._first_derivative_code['f'](**self._arguments) + + def d2logpdf_dlink2(self, inv_link_f, y, Y_metadata=None): + """ + Hessian of log likelihood given inverse link of latent variables with respect to that inverse link. + i.e. second derivative logpdf at y given inv_link(f_i) and inv_link(f_j) w.r.t inv_link(f_i) and inv_link(f_j). + + + :param inv_link_f: inverse link of the latent variables. + :type inv_link_f: Nx1 array + :param y: data + :type y: Nx1 array + :param Y_metadata: Y_metadata which is not used in student t distribution + :returns: Diagonal of Hessian matrix (second derivative of likelihood evaluated at points f) + :rtype: Nx1 array + + .. Note:: + Returns diagonal of Hessian, since every where else it is + 0, as the likelihood factorizes over cases (the + distribution for y_i depends only on link(f_i) not on + link(f_(j!=i)) + """ + assert np.atleast_1d(inv_link_f).shape == np.atleast_1d(y).shape + self._arguments_update(inv_link_f, y) + return self._second_derivative_code['f'](**self._arguments) + + def d3logpdf_dlink3(self, inv_link_f, y, Y_metadata=None): + assert np.atleast_1d(inv_link_f).shape == np.atleast_1d(y).shape + self._arguments_update(inv_link_f, y) + return self._third_derivative_code['f'](**self._arguments) + raise NotImplementedError + + def dlogpdf_link_dtheta(self, inv_link_f, y, Y_metadata=None): + assert np.atleast_1d(inv_link_f).shape == np.atleast_1d(y).shape + self._arguments_update(inv_link_f, y) + return np.asarray([self._first_derivative_code[theta.name](**self._arguments).sum() for theta in self._sp_theta]) + + def dlogpdf_dlink_dtheta(self, inv_link_f, y, Y_metadata=None): + assert np.atleast_1d(inv_link_f).shape == np.atleast_1d(y).shape + self._arguments_update(inv_link_f, y) + return np.asarray([self._second_derivative_code[theta.name](**self._arguments).sum() for theta in self._sp_theta]) + + def d2logpdf_dlink2_dtheta(self, inv_link_f, y, Y_metadata=None): + assert np.atleast_1d(inv_link_f).shape == np.atleast_1d(y).shape + self._arguments_update(inv_link_f, y) + return np.asarray([self._third_derivative_code[theta.name](**self._arguments).sum() for theta in self._sp_theta]) + + def predictive_mean(self, mu, sigma, Y_metadata=None): + raise NotImplementedError + + def predictive_variance(self, mu,variance, predictive_mean=None, Y_metadata=None): + raise NotImplementedError + + def conditional_mean(self, gp): + raise NotImplementedError + + def conditional_variance(self, gp): + raise NotImplementedError + + def samples(self, gp, Y_metadata=None): + raise NotImplementedError