From 61872fb3159b1a02f897c905556fbf2bd4a9f240 Mon Sep 17 00:00:00 2001 From: Akash Kumar Dhaka Date: Thu, 1 Jun 2017 15:03:19 +0300 Subject: [PATCH 1/3] fresh branch for these new likelihood functions along with test code- will work atleast with La[lace approximation ... --- GPy/likelihoods/__init__.py | 4 +- GPy/likelihoods/loggaussian.py | 304 ++++++++++++++++++++++++++++ GPy/likelihoods/loglogistic.py | 338 ++++++++++++++++++++++++++++++++ GPy/likelihoods/weibull.py | 224 +++++++++++++++++++++ GPy/testing/likelihood_tests.py | 32 +++ 5 files changed, 901 insertions(+), 1 deletion(-) create mode 100644 GPy/likelihoods/loggaussian.py create mode 100644 GPy/likelihoods/loglogistic.py create mode 100644 GPy/likelihoods/weibull.py diff --git a/GPy/likelihoods/__init__.py b/GPy/likelihoods/__init__.py index 627ef39f..aa9464b5 100644 --- a/GPy/likelihoods/__init__.py +++ b/GPy/likelihoods/__init__.py @@ -7,4 +7,6 @@ from .student_t import StudentT from .likelihood import Likelihood from .mixed_noise import MixedNoise from .binomial import Binomial - +from .weibull import Weibull +from .loglogistic import LogLogistic +from .loggaussian import LogGaussian diff --git a/GPy/likelihoods/loggaussian.py b/GPy/likelihoods/loggaussian.py new file mode 100644 index 00000000..58a53aed --- /dev/null +++ b/GPy/likelihoods/loggaussian.py @@ -0,0 +1,304 @@ +# Copyright (c) 2012 - 2014, GPy authors (see AUTHORS.txt). +# Licensed under the BSD 3-clause license (see LICENSE.txt) + +import numpy as np +from scipy import stats, special +from ..core.parameterization import Param +from ..core.parameterization.transformations import Logexp +from . import link_functions +from .likelihood import Likelihood + + + +class LogGaussian(Likelihood): + """ + .. math:: + $$ p(y_{i}|f_{i}, z_{i}) = \\prod_{i=1}^{n} (\\frac{ry^{r-1}}{\\exp{f(x_{i})}})^{1-z_i} (1 + (\\frac{y}{\\exp(f(x_{i}))})^{r})^{z_i-2} $$ + + .. note: + where z_{i} is the censoring indicator- 0 for non-censored data, and 1 for censored data. + + + """ + def __init__(self,gp_link=None, sigma=1.): + if gp_link is None: + gp_link = link_functions.Identity() + # gp_link = link_functions.Log() + + super(LogGaussian, self).__init__(gp_link, name='loggaussian') + + self.sigma = Param('sigma', sigma, Logexp()) + self.variance = Param('variance', sigma**2, Logexp()) + self.link_parameter(self.variance) + # self.link_parameter() + + def pdf_link(self, link_f, y, Y_metadata=None): + """ + :param link_f: latent variables link(f) + :type link_f: Nx1 array + :param y: data + :type y: Nx1 array + :param Y_metadata: includes censoring information in dictionary key 'censored' + :returns: likelihood evaluated for this point + :rtype: float + """ + return np.exp(self.logpdf_link(link_f, y, Y_metadata=Y_metadata)) + + def logpdf_link(self, link_f, y, Y_metadata=None): + """ + :param link_f: latent variables (link(f)) + :type link_f: Nx1 array + :param y: data + :type y: Nx1 array + :param Y_metadata: includes censoring information in dictionary key 'censored' + :returns: likelihood evaluated for this point + :rtype: float + """ + assert np.atleast_1d(link_f).shape == np.atleast_1d(y).shape + c = np.zeros_like(y) + if Y_metadata is not None and 'censored' in Y_metadata.keys(): + c = Y_metadata['censored'] + + uncensored = (1-c)* (-0.5*np.log(2*np.pi*self.variance) - np.log(y) - (np.log(y)-link_f)**2 /(2*self.variance) ) + censored = c*np.log( 1 - stats.norm.cdf((np.log(y) - link_f)/np.sqrt(self.variance)) ) + logpdf = uncensored + censored + return logpdf + + def dlogpdf_dlink(self, link_f, y, Y_metadata=None): + """ + derivative of logpdf wrt link_f param + .. math:: + + :param link_f: latent variables link(f) + :type link_f: Nx1 array + :param y: data + :type y: Nx1 array + :param Y_metadata: includes censoring information in dictionary key 'censored' + :returns: likelihood evaluated for this point + :rtype: float + """ + assert np.atleast_1d(link_f).shape == np.atleast_1d(y).shape + c = np.zeros_like(y) + if Y_metadata is not None and 'censored' in Y_metadata.keys(): + c = Y_metadata['censored'] + + val = np.log(y) - link_f + val_scaled = val/np.sqrt(self.variance) + val_scaled2 = val/self.variance + uncensored = (1-c)*(val_scaled2) + a = (1- stats.norm.cdf(val_scaled)) + # llg(z) = 1. / (1 - norm_cdf(r / sqrt(s2))). * (1 / sqrt(2 * pi * s2). * exp(-1 / (2. * s2). * r. ^ 2)); + censored = c*( 1./a) * (np.exp(-1.* val**2 /(2*self.variance)) / np.sqrt(2*np.pi*self.variance)) + # censored = c * (1. / (1 - stats.norm.cdf(val_scaled))) * (stats.norm.pdf(val_scaled)) + gradient = uncensored + censored + return gradient + + def d2logpdf_dlink2(self, link_f, y, Y_metadata=None): + """ + Hessian at y, given link(f), w.r.t link(f) + i.e. second derivative logpdf at y given link(f_i) and link(f_j) w.r.t link(f_i) and link(f_j) + The hessian will be 0 unless i == j + + .. math:: + + + :param link_f: latent variables link(f) + :type link_f: Nx1 array + :param y: data + :type y: Nx1 array + :param Y_metadata: includes censoring information in dictionary key 'censored' + :returns: Diagonal of hessian matrix (second derivative of likelihood evaluated at points f) + :rtype: Nx1 array + + .. Note:: + Will return 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)) + """ + # c = Y_metadata['censored'] + # c = np.zeros((y.shape[0],)) + c = np.zeros_like(y) + if Y_metadata is not None and 'censored' in Y_metadata.keys(): + c = Y_metadata['censored'] + + val = np.log(y) - link_f + val_scaled = val/np.sqrt(self.variance) + val_scaled2 = val/self.variance + a = (1 - stats.norm.cdf(val_scaled)) + uncensored = (1-c) *(-1)/self.variance + censored = c*(-np.exp(-val**2/self.variance) / ( 2*np.pi*self.variance*(a**2) ) + + val*np.exp(-(val**2)/(2*self.variance))/( np.sqrt(2*np.pi)*self.variance**(3/2.)*a) ) + hessian = censored + uncensored + return hessian + + def d3logpdf_dlink3(self, link_f, y, Y_metadata=None): + """ + Gradient of the log-likelihood function at y given f, w.r.t shape parameter + + .. math:: + + :param inv_link_f: latent variables link(f) + :type inv_link_f: Nx1 array + :param y: data + :type y: Nx1 array + :param Y_metadata: includes censoring information in dictionary key 'censored' + :returns: derivative of likelihood evaluated at points f w.r.t variance parameter + :rtype: float + """ + c = np.zeros_like(y) + if Y_metadata is not None and 'censored' in Y_metadata.keys(): + c = Y_metadata['censored'] + + val = np.log(y) - link_f + val_scaled = val/np.sqrt(self.variance) + val_scaled2 = val/self.variance + a = (1 - stats.norm.cdf(val_scaled)) + uncensored = 0 + censored = c *( 2*np.exp(-3*(val**2)/(2*self.variance)) / ((a**3)*(2*np.pi*self.variance)**(3/2.)) + - val*np.exp(-(val**2)/self.variance)/ ( (a**2)*np.pi*self.variance**2) + - val*np.exp(-(val**2)/self.variance)/ ( (a**2)*2*np.pi*self.variance**2) + - np.exp(-(val**2)/(2*self.variance))/ ( a*(self.variance**(1.50))*np.sqrt(2*np.pi)) + + (val**2)*np.exp(-(val**2)/(2*self.variance))/ ( a*np.sqrt(2*np.pi*self.variance)*self.variance**2 ) ) + d3pdf_dlink3 = uncensored + censored + return d3pdf_dlink3 + + def dlogpdf_link_dvar(self, link_f, y, Y_metadata=None): + """ + Gradient of the log-likelihood function at y given f, w.r.t variance parameter + + .. math:: + + :param inv_link_f: latent variables link(f) + :type inv_link_f: Nx1 array + :param y: data + :type y: Nx1 array + :param Y_metadata: includes censoring information in dictionary key 'censored' + :returns: derivative of likelihood evaluated at points f w.r.t variance parameter + :rtype: float + """ + + c = np.zeros_like(y) + if Y_metadata is not None and 'censored' in Y_metadata.keys(): + c = Y_metadata['censored'] + + val = np.log(y) - link_f + val_scaled = val/np.sqrt(self.variance) + val_scaled2 = val/self.variance + a = (1 - stats.norm.cdf(val_scaled)) + uncensored = (1-c)*(-0.5/self.variance + (val**2)/(2*(self.variance**2)) ) + censored = c *( val*np.exp(-val**2/ (2*self.variance)) / (a*np.sqrt(2*np.pi)*2*(self.variance**(1.5))) ) + dlogpdf_dvar = uncensored + censored + # dlogpdf_dvar = dlogpdf_dvar*self.variance + return dlogpdf_dvar + + def dlogpdf_dlink_dvar(self, link_f, y, Y_metadata=None): + """ + :param link_f: latent variables link(f) + :type link_f: Nx1 array + :param y: data + :type y: Nx1 array + :param Y_metadata: Y_metadata not used in gaussian + :returns: derivative of log likelihood evaluated at points link(f) w.r.t variance parameter + :rtype: Nx1 array + """ + c = np.zeros_like(y) + if Y_metadata is not None and 'censored' in Y_metadata.keys(): + c = Y_metadata['censored'] + + val = np.log(y) - link_f + val_scaled = val/np.sqrt(self.variance) + val_scaled2 = val/self.variance + a = (1 - stats.norm.cdf(val_scaled)) + uncensored = (1-c)*(-val/(self.variance**2)) + censored = c * (-val*np.exp(-val**2/self.variance)/( 4*np.pi*(self.variance**2)*(a**2)) + + (-1 + (val**2)/self.variance)*np.exp(-val**2/(2*self.variance) ) / + ( a*(np.sqrt(2.*np.pi)*2*self.variance**1.5)) ) + dlik_grad_dsigma = uncensored + censored + # dlik_grad_dsigma = dlik_grad_dsigma*self.variance + return dlik_grad_dsigma + + def d2logpdf_dlink2_dvar(self, link_f, y, Y_metadata=None): + """ + :param link_f: latent variables link(f) + :type link_f: Nx1 array + :param y: data + :type y: Nx1 array + :param Y_metadata: Y_metadata not used in gaussian + :returns: derivative of log likelihood evaluated at points link(f) w.r.t variance parameter + :rtype: Nx1 array + """ + c = np.zeros_like(y) + if Y_metadata is not None and 'censored' in Y_metadata.keys(): + c = Y_metadata['censored'] + val = np.log(y) - link_f + val_scaled = val/np.sqrt(self.variance) + val_scaled2 = val/self.variance + a = (1 - stats.norm.cdf(val_scaled)) + uncensored = (1-c)*( 1./(self.variance**2) ) + censored = c*( val*np.exp(-3*(val**2)/(2*self.variance) )/ ((a**3)*np.sqrt(8*np.pi**3)*self.variance**(5/2.)) + + np.exp(-val**2/self.variance)/((a**2)*4*np.pi*self.variance**2) + - np.exp(-val**2/self.variance)*val**2 / ((a**2)*2*np.pi*self.variance**3) + + np.exp(-val**2/self.variance)/ ( (a**2)*4*np.pi*self.variance**2) + - np.exp(-val**2/ (2*self.variance))*val / ( a*np.sqrt(2*np.pi)*2*self.variance**(5/2.)) + - np.exp(-val**2/self.variance)*(val**2) / ((a**2)*4*np.pi*self.variance**3) + - np.exp(-val**2/ (2*self.variance))*val/ (a*np.sqrt(2*np.pi)*self.variance**(5/2.)) + + np.exp(-val**2/ (2*self.variance))*(val**3) / (a*np.sqrt(2*np.pi)*2*self.variance**(7/2.)) ) + dlik_hess_dsigma = uncensored + censored + return dlik_hess_dsigma + + def dlogpdf_link_dtheta(self, f, y, Y_metadata=None): + """ + :param link_f: latent variables link(f) + :type link_f: Nx1 array + :param y: data + :type y: Nx1 array + :param Y_metadata: Y_metadata not used in gaussian + :returns: derivative of log likelihood evaluated at points link(f) w.r.t variance parameter + :rtype: Nx1 array + """ + dlogpdf_dtheta = np.zeros((self.size, f.shape[0], f.shape[1])) + dlogpdf_dtheta[0,:,:] = self.dlogpdf_link_dvar(f,y,Y_metadata=Y_metadata) + return dlogpdf_dtheta + + def dlogpdf_dlink_dtheta(self, f, y, Y_metadata=None): + """ + :param link_f: latent variables link(f) + :type link_f: Nx1 array + :param y: data + :type y: Nx1 array + :param Y_metadata: Y_metadata not used in gaussian + :returns: derivative of log likelihood evaluated at points link(f) w.r.t variance parameter + :rtype: Nx1 array + """ + dlogpdf_dlink_dtheta = np.zeros((self.size, f.shape[0], f.shape[1])) + dlogpdf_dlink_dtheta[0,:,:] = self.dlogpdf_dlink_dvar(f,y,Y_metadata=Y_metadata) + return dlogpdf_dlink_dtheta + + def d2logpdf_dlink2_dtheta(self, f, y, Y_metadata=None): + """ + :param link_f: latent variables link(f) + :type link_f: Nx1 array + :param y: data + :type y: Nx1 array + :param Y_metadata: Y_metadata not used in gaussian + :returns: derivative of log likelihood evaluated at points link(f) w.r.t variance parameter + :rtype: Nx1 array + """ + d2logpdf_dlink2_dtheta = np.zeros((self.size, f.shape[0], f.shape[1])) + d2logpdf_dlink2_dtheta[0,:,:] = self.d2logpdf_dlink2_dvar(f,y,Y_metadata=Y_metadata) + return d2logpdf_dlink2_dtheta + + def update_gradients(self, grads): + """ + Pull out the gradients, be careful as the order must match the order + in which the parameters are added + """ + self.variance.gradient = grads[0] + + def samples(self, gp, Y_metadata=None): + """ + Returns a set of samples of observations based on a given value of the latent variable. + + :param gp: latent variable + """ + orig_shape = gp.shape + gp = gp.flatten() diff --git a/GPy/likelihoods/loglogistic.py b/GPy/likelihoods/loglogistic.py new file mode 100644 index 00000000..2695bd69 --- /dev/null +++ b/GPy/likelihoods/loglogistic.py @@ -0,0 +1,338 @@ +from __future__ import division +# Copyright (c) 2015 Alan Saul +# Licensed under the BSD 3-clause license (see LICENSE.txt) + +import numpy as np +from scipy import stats,special +import scipy as sp +from GPy.likelihoods.link_functions import Log, Identity +from GPy.likelihoods.likelihood import Likelihood +from GPy.core.parameterization import Param +from GPy.core.parameterization.transformations import Logexp + +class LogLogistic(Likelihood): + """ + .. math:: + $$ p(y_{i}|f_{i}, z_{i}) = \\prod_{i=1}^{n} (\\frac{ry^{r-1}}{\\exp{f(x_{i})}})^{1-z_i} (1 + (\\frac{y}{\\exp(f(x_{i}))})^{r})^{z_i-2} $$ + + .. note: + where z_{i} is the censoring indicator- 0 for non-censored data, and 1 for censored data. + """ + + def __init__(self, gp_link=None, r=1.0): + if gp_link is None: + #Parameterised not as link_f but as f + gp_link = Log() + + super(LogLogistic, self).__init__(gp_link, name='LogLogistic') + self.r = Param('r_log_shape', float(r), Logexp()) + self.link_parameter(self.r) + # self.censored = 'censored' + + + + def pdf_link(self, link_f, y, Y_metadata=None): + """ + Likelihood function given link(f) + + .. math:: + + :param link_f: latent variables link(f) + :type link_f: Nx1 array + :param y: data + :type y: Nx1 array + :param Y_metadata: includes censoring information in dictionary key 'censored' + :returns: likelihood evaluated for this point + :rtype: float + """ + return np.exp(self.logpdf_link(link_f, y, Y_metadata=Y_metadata)) + + + def logpdf_link(self, link_f, y, Y_metadata=None): + """ + Log Likelihood Function given link(f) + + .. math:: + + + :param link_f: latent variables (link(f)) + :type link_f: Nx1 array + :param y: data + :type y: Nx1 array + :param Y_metadata: includes censoring information in dictionary key 'censored' + :returns: likelihood evaluated for this point + :rtype: float + + """ + # c = np.zeros((y.shape[0],)) + c = np.zeros_like(link_f) + if Y_metadata is not None and 'censored' in Y_metadata.keys(): + c = Y_metadata['censored'] + + link_f = np.clip(link_f, 1e-150, 1e100) + # y_link_f = y/link_f + # y_link_f_r = y_link_f**self.r + # y_link_f_r = np.clip(y**self.r, 1e-150, 1e200) / np.clip(link_f**self.r, 1e-150, 1e200) + # y_link_f_r = np.clip((y/link_f)**self.r, 1e-150, 1e200) + y_r = np.clip(y**self.r, 1e-150, 1e200) + link_f_r = np.clip(link_f**self.r, 1e-150, 1e200) + y_link_f_r = np.clip(y_r / link_f_r, 1e-150, 1e200) + #uncensored = (1-c)*(np.log(self.r) + (self.r+1)*np.log(y) - self.r*np.log(link_f) - 2*np.log1p(y_link_f_r)) + #uncensored = (1-c)*(np.log((self.r/link_f)*y_link_f**(self.r-1)) - 2*np.log1p(y_link_f_r)) + + # clever way tp break it into censored and uncensored-parts .. + uncensored = (1-c)*(np.log(self.r) + (self.r-1)*np.log(y) - self.r*np.log(link_f) - 2*np.log1p(y_link_f_r)) + censored = (c)*(-np.log1p(y_link_f_r)) + # + return uncensored + censored + # return uncensored + + def dlogpdf_dlink(self, link_f, y, Y_metadata=None): + """ + Gradient of the log likelihood function at y, given link(f) w.r.t link(f) + + .. math:: + + :param link_f: latent variables (f) + :type link_f: Nx1 array + :param y: data + :type y: Nx1 array + :param Y_metadata: includes censoring information in dictionary key 'censored' + :returns: gradient of likelihood evaluated at points + :rtype: Nx1 array + + """ + # c = Y_metadata['censored'] + # for debugging + # c = np.zeros((y.shape[0],)) + c = np.zeros_like(link_f) + + if Y_metadata is not None and 'censored' in Y_metadata.keys(): + c = Y_metadata['censored'] + + #y_link_f = y/link_f + #y_link_f_r = y_link_f**self.r + y_link_f_r = np.clip(y**self.r, 1e-150, 1e200) / np.clip(link_f**self.r, 1e-150, 1e200) + + #In terms of link_f + # uncensored = (1-c)*( (2*self.r*y**r)/(link_f**self.r + y**self.r) - link_f*self.r) + uncensored = (1-c)*self.r*(y_link_f_r - 1)/(link_f*(1 + y_link_f_r)) + censored = c*(self.r*y_link_f_r/(link_f*y_link_f_r + link_f)) + return uncensored + censored + # return uncensored + + def d2logpdf_dlink2(self, link_f, y, Y_metadata=None): + """ + Hessian at y, given link(f), w.r.t link(f) + i.e. second derivative logpdf at y given link(f_i) and link(f_j) w.r.t link(f_i) and link(f_j) + The hessian will be 0 unless i == j + + .. math:: + + + :param link_f: latent variables link(f) + :type link_f: Nx1 array + :param y: data + :type y: Nx1 array + :param Y_metadata: includes censoring information in dictionary key 'censored' + :returns: Diagonal of hessian matrix (second derivative of likelihood evaluated at points f) + :rtype: Nx1 array + + .. Note:: + Will return 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)) + """ + # c = Y_metadata['censored'] + # c = np.zeros((y.shape[0],)) + c = np.zeros_like(link_f) + + if Y_metadata is not None and 'censored' in Y_metadata.keys(): + c = Y_metadata['censored'] + + y_link_f = y/link_f + y_link_f_r = y_link_f**self.r + + #In terms of link_f + censored = c*(-self.r*y_link_f_r*(y_link_f_r + self.r + 1)/((link_f**2)*(y_link_f_r + 1)**2)) + uncensored = (1-c)*(-self.r*(2*self.r*y_link_f_r + y_link_f**(2*self.r) - 1) / ((link_f**2)*(1+ y_link_f_r)**2)) + hess = censored + uncensored + return hess + + def d3logpdf_dlink3(self, link_f, y, Y_metadata=None): + """ + Third order derivative log-likelihood function at y given link(f) w.r.t link(f) + + .. math:: + + + :param link_f: latent variables link(f) + :type link_f: Nx1 array + :param y: data + :type y: Nx1 array + :param Y_metadata: includes censoring information in dictionary key 'censored' + :returns: third derivative of likelihood evaluated at points f + :rtype: Nx1 array + """ + # c = Y_metadata['censored'] + # for debugging + # c = np.zeros((y.shape[0],)) + c = np.zeros_like(link_f) + + if Y_metadata is not None and 'censored' in Y_metadata.keys(): + c = Y_metadata['censored'] + y_link_f = y/link_f + y_link_f_r = y_link_f**self.r + + #In terms of link_f + censored = c*(self.r*y_link_f_r*(((self.r**2)*(-(y_link_f_r - 1))) + 3*self.r*(y_link_f_r + 1) + 2*(y_link_f_r + 1)**2) + / ((link_f**3)*(y_link_f_r + 1)**3)) + uncensored = (1-c)*(2*self.r*(-(self.r**2)*(y_link_f_r -1)*y_link_f_r + 3*self.r*(y_link_f_r + 1)*y_link_f_r + (y_link_f_r - 1)*(y_link_f_r + 1)**2) + / ((link_f**3)*(y_link_f_r + 1)**3)) + + d3lik_dlink3 = censored + uncensored + return d3lik_dlink3 + + def dlogpdf_link_dr(self, inv_link_f, y, Y_metadata=None): + """ + Gradient of the log-likelihood function at y given f, w.r.t shape parameter + + .. math:: + + :param inv_link_f: latent variables link(f) + :type inv_link_f: Nx1 array + :param y: data + :type y: Nx1 array + :param Y_metadata: includes censoring information in dictionary key 'censored' + :returns: derivative of likelihood evaluated at points f w.r.t variance parameter + :rtype: float + """ + # c = Y_metadata['censored'] + # c = np.zeros((y.shape[0],)) + c = np.zeros_like(y) + if Y_metadata is not None and 'censored' in Y_metadata.keys(): + c = Y_metadata['censored'] + + link_f = inv_link_f #FIXME: Change names consistently... + y_link_f = y/link_f + log_y_link_f = np.log(y) - np.log(link_f) + y_link_f_r = y_link_f**self.r + + #In terms of link_f + censored = c*(-y_link_f_r*log_y_link_f/(1 + y_link_f_r)) + uncensored = (1-c)*(1./self.r + np.log(y) - np.log(link_f) - (2*y_link_f_r*log_y_link_f) / (1 + y_link_f_r)) + + dlogpdf_dr = censored + uncensored + return dlogpdf_dr + + def dlogpdf_dlink_dr(self, inv_link_f, y, Y_metadata=None): + """ + Derivative of the dlogpdf_dlink w.r.t shape parameter + + .. math:: + + :param inv_link_f: latent variables inv_link_f + :type inv_link_f: Nx1 array + :param y: data + :type y: Nx1 array + :param Y_metadata: includes censoring information in dictionary key 'censored' + :returns: derivative of likelihood evaluated at points f w.r.t variance parameter + :rtype: Nx1 array + """ + # c = np.zeros((y.shape[0],)) + c = np.zeros_like(y) + if Y_metadata is not None and 'censored' in Y_metadata.keys(): + c = Y_metadata['censored'] + link_f = inv_link_f + y_link_f = y/link_f + y_link_f_r = y_link_f**self.r + log_y_link_f = np.log(y) - np.log(link_f) + + #In terms of link_f + censored = c*(y_link_f_r*(y_link_f_r + self.r*log_y_link_f + 1)/(link_f*(y_link_f_r + 1)**2)) + uncensored = (1-c)*(y_link_f**(2*self.r) + 2*self.r*y_link_f_r*log_y_link_f - 1) / (link_f*(1 + y_link_f_r)**2) + + # dlogpdf_dlink_dr = uncensored + dlogpdf_dlink_dr = censored + uncensored + return dlogpdf_dlink_dr + + def d2logpdf_dlink2_dr(self, inv_link_f, y, Y_metadata=None): + """ + Gradient of the hessian (d2logpdf_dlink2) w.r.t shape parameter + + .. math:: + + :param inv_link_f: latent variables link(f) + :type inv_link_f: Nx1 array + :param y: data + :type y: Nx1 array + :param Y_metadata: includes censoring information in dictionary key 'censored' + :returns: derivative of hessian evaluated at points f and f_j w.r.t variance parameter + :rtype: Nx1 array + """ + # c = Y_metadata['censored'] + + # c = np.zeros((y.shape[0],)) + c = np.zeros_like(y) + if Y_metadata is not None and 'censored' in Y_metadata.keys(): + c = Y_metadata['censored'] + link_f = inv_link_f + y_link_f = y/link_f + y_link_f_r = y_link_f**self.r + log_y_link_f = np.log(y) - np.log(link_f) + + #In terms of link_f + y_link_f_2r = y_link_f**(2*self.r) + denom2 = (link_f**2)*(1 + y_link_f_r)**2 + denom3 = (link_f**2)*(1 + y_link_f_r)**3 + + censored = c*(-((y_link_f_r + self.r + 1)*y_link_f_r)/denom2 + -(self.r*(y_link_f_r + self.r + 1)*y_link_f_r*log_y_link_f)/denom2 + -(self.r*y_link_f_r*(y_link_f_r*log_y_link_f + 1))/denom2 + +(2*self.r*(y_link_f_r + self.r + 1)*y_link_f_2r*log_y_link_f)/denom3 + ) + + uncensored = (1-c)*(-(2*self.r*y_link_f_r + y_link_f_2r - 1)/denom2 + -(self.r*(2*y_link_f_r + 2*self.r*y_link_f_r*log_y_link_f + 2*y_link_f_2r*log_y_link_f)/denom2) + +(2*self.r*(2*self.r*y_link_f_r + y_link_f_2r - 1)*y_link_f_r*log_y_link_f)/denom3 + ) + d2logpdf_dlink2_dr = censored + uncensored + + return d2logpdf_dlink2_dr + + def dlogpdf_link_dtheta(self, f, y, Y_metadata=None): + dlogpdf_dtheta = np.zeros((self.size, f.shape[0], f.shape[1])) + dlogpdf_dtheta[0, :, :] = self.dlogpdf_link_dr(f, y, Y_metadata=Y_metadata) + return dlogpdf_dtheta + + def dlogpdf_dlink_dtheta(self, f, y, Y_metadata=None): + dlogpdf_dlink_dtheta = np.zeros((self.size, f.shape[0], f.shape[1])) + dlogpdf_dlink_dtheta[0, :, :] = self.dlogpdf_dlink_dr(f, y, Y_metadata=Y_metadata) + return dlogpdf_dlink_dtheta + + def d2logpdf_dlink2_dtheta(self, f, y, Y_metadata=None): + d2logpdf_dlink2_dtheta = np.zeros((self.size, f.shape[0], f.shape[1])) + d2logpdf_dlink2_dtheta[0,:, :] = self.d2logpdf_dlink2_dr(f, y, Y_metadata=Y_metadata) + return d2logpdf_dlink2_dtheta + + def update_gradients(self, grads): + """ + Pull out the gradients, be careful as the order must match the order + in which the parameters are added + """ + self.r.gradient = grads[0] + + def samples(self, gp, Y_metadata=None): + """ + Returns a set of samples of observations based on a given value of the latent variable. + + :param gp: latent variable + """ + orig_shape = gp.shape + gp = gp.flatten() + #rs = np.ones_like(gp)*self.r + #scales = np.ones_like(gp)*np.sqrt(self.sigma2) + #Ysim = sp.stats.fisk.rvs(rs, scale=self.gp_link.transf(gp)) + Ysim = np.array([sp.stats.fisk.rvs(self.r, loc=0, scale=self.gp_link.transf(f)) for f in gp]) + #np.random.fisk(self.gp_link.transf(gp), c=self.r) + return Ysim.reshape(orig_shape) + diff --git a/GPy/likelihoods/weibull.py b/GPy/likelihoods/weibull.py new file mode 100644 index 00000000..8ee9ba7f --- /dev/null +++ b/GPy/likelihoods/weibull.py @@ -0,0 +1,224 @@ +# Copyright (c) 2012 - 2014, GPy authors (see AUTHORS.txt). +# Licensed under the BSD 3-clause license (see LICENSE.txt) + + +import numpy as np +from scipy import stats,special +import scipy as sp +from ..core.parameterization import Param +from ..core.parameterization.transformations import Logexp +from . import link_functions +from .likelihood import Likelihood + +class Weibull(Likelihood): + """ + .. math:: + $$ p(y_{i}|f_{i}, z_{i}) = \\prod_{i=1}^{n} [ r^{1-z_{i}}\\exp(-(1-z_{i})f(x_{i}))y_{i}^{(1-z_{i})(r-1)}\\exp(-y_{i}^{r}/\\exp(f(x_{i})))] $$ + + .. note: + where z_{i} is the censoring indicator- 0 for non-censored data, and 1 for censored data and r is the shape parameter. + """ + def __init__(self,gp_link=None,beta=1.): + if gp_link is None: + gp_link = link_functions.Log() + # gp_link = link_functions.Identity() + super(Weibull, self).__init__(gp_link, name='Weibull') + + self.r = Param('r_weibull_shape', float(beta), Logexp()) + self.link_parameter(self.r) + + # self.r.fix() + + def pdf_link(self, link_f, y, Y_metadata=None): + """ + Likelihood function given link(f) + + .. math:: + p(y_{i}|\\lambda(f_{i})) = \\frac{\\beta^{\\alpha_{i}}}{\\Gamma(\\alpha_{i})}y_{i}^{\\alpha_{i}-1}e^{-\\beta y_{i}}\\\\ + \\alpha_{i} = \\beta y_{i} + + :param link_f: latent variables link(f) + :type link_f: Nx1 array + :param y: data + :type y: Nx1 array + :param Y_metadata: includes censoring information in dictionary key 'censored' + :returns: likelihood evaluated for this point + :rtype: float + """ + assert np.atleast_1d(link_f).shape == np.atleast_1d(y).shape + + c = np.zeros((link_f.shape[0],)) + + log_objective = np.log(self.r) + (self.r - 1) * np.log(y) - link_f - (np.exp(-link_f)*(y ** self.r)) + # log_objective = stats.weibull_min.pdf(y,c=self.beta,loc=link_f,scale=1.) + return np.exp(log_objective) + + def logpdf_link(self, link_f, y, Y_metadata=None): + """ + Log Likelihood Function given link(f) + + .. math:: + \\ln p(y_{i}|\lambda(f_{i})) = \\alpha_{i}\\log \\beta - \\log \\Gamma(\\alpha_{i}) + (\\alpha_{i} - 1)\\log y_{i} - \\beta y_{i}\\\\ + \\alpha_{i} = \\beta y_{i} + + :param link_f: latent variables (link(f)) + :type link_f: Nx1 array + :param y: data + :type y: Nx1 array + :param Y_metadata: includes censoring information in dictionary key 'censored' + :returns: likelihood evaluated for this point + :rtype: float + + """ + #alpha = self.gp_link.transf(gp)*self.beta sum(log(a) + (a-1).*log(y)- f - exp(-f).*y.^a) + #return (1. - alpha)*np.log(obs) + self.beta*obs - alpha * np.log(self.beta) + np.log(special.gamma(alpha)) + assert np.atleast_1d(link_f).shape == np.atleast_1d(y).shape + log_objective = np.log(self.r) + (self.r - 1) * np.log(y) - link_f - (np.exp(-link_f) * (y ** self.r)) + return log_objective + + def dlogpdf_dlink(self, link_f, y, Y_metadata=None): + """ + Gradient of the log likelihood function at y, given link(f) w.r.t link(f) + + .. math:: + \\frac{d \\ln p(y_{i}|\\lambda(f_{i}))}{d\\lambda(f)} = \\beta (\\log \\beta y_{i}) - \\Psi(\\alpha_{i})\\beta\\\\ + \\alpha_{i} = \\beta y_{i} + + :param link_f: latent variables (f) + :type link_f: Nx1 array + :param y: data + :type y: Nx1 array + :param Y_metadata: includes censoring information in dictionary key 'censored' + :returns: likelihood evaluated for this point + :rtype: Nx1 array + + """ + # grad = (1. - self.beta) / (y - link_f) + + grad = -1 + np.exp(-link_f)*(y ** self.r) + return grad + + def d2logpdf_dlink2(self, link_f, y, Y_metadata=None): + """ + Hessian at y, given link(f), w.r.t link(f) + i.e. second derivative logpdf at y given link(f_i) and link(f_j) w.r.t link(f_i) and link(f_j) + The hessian will be 0 unless i == j + + .. math:: + \\frac{d^{2} \\ln p(y_{i}|\lambda(f_{i}))}{d^{2}\\lambda(f)} = -\\beta^{2}\\frac{d\\Psi(\\alpha_{i})}{d\\alpha_{i}}\\\\ + \\alpha_{i} = \\beta y_{i} + + :param link_f: latent variables link(f) + :type link_f: Nx1 array + :param y: data + :type y: Nx1 array + :param Y_metadata: includes censoring information in dictionary key 'censored' + :returns: likelihood evaluated for this point + :rtype: Nx1 array + + .. Note:: + Will return 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)) + """ + # hess = (self.beta - 1.) / (y - link_f)**2 + hess = -(y ** self.r) * np.exp(-link_f) + return hess + + def d3logpdf_dlink3(self, link_f, y, Y_metadata=None): + """ + Third order derivative log-likelihood function at y given link(f) w.r.t link(f) + + .. math:: + \\frac{d^{3} \\ln p(y_{i}|\lambda(f_{i}))}{d^{3}\\lambda(f)} = -\\beta^{3}\\frac{d^{2}\\Psi(\\alpha_{i})}{d\\alpha_{i}}\\\\ + \\alpha_{i} = \\beta y_{i} + + :param link_f: latent variables link(f) + :type link_f: Nx1 array + :param y: data + :type y: Nx1 array + :param Y_metadata: Y_metadata which is not used in gamma distribution + :returns: third derivative of likelihood evaluated at points f + :rtype: Nx1 array + """ + # d3lik_dlink3 = (1. - self.beta) / (y - link_f)**3 + d3lik_dlink3 = (y ** self.r) * np.exp(-link_f) + return d3lik_dlink3 + + def exact_inference_gradients(self, dL_dKdiag,Y_metadata=None): + return np.zeros(self.size) + + def dlogpdf_link_dr(self, inv_link_f, y, Y_metadata=None): + """ + Gradient of the log-likelihood function at y given f, w.r.t shape parameter + + .. math:: + + :param inv_link_f: latent variables link(f) + :type inv_link_f: Nx1 array + :param y: data + :type y: Nx1 array + :param Y_metadata: includes censoring information in dictionary key 'censored' + :returns: derivative of likelihood evaluated at points f w.r.t variance parameter + :rtype: float + """ + dlogpdf_dr = 1./self.r + np.log(y) - np.exp(-inv_link_f)*(y**self.r)*np.log(y) + return dlogpdf_dr + + def dlogpdf_dlink_dr(self, link_f, y, Y_metadata=None): + """ + First order derivative derivative of loglikelihood wrt r:shape parameter + + :param link_f: latent variables link(f) + :type link_f: Nx1 array + :param y: data + :type y: Nx1 array + :param Y_metadata: Y_metadata which is not used in gamma distribution + :returns: third derivative of likelihood evaluated at points f + :rtype: Nx1 array + """ + # dlogpdf_dlink_dr = self.beta * y**(self.beta - 1) * np.exp(-link_f) + dlogpdf_dlink_dr = np.exp(-link_f)* (y ** self.r) * np.log(y) + return dlogpdf_dlink_dr + + def d2logpdf_dlink2_dr(self, link_f, y, Y_metadata=None): + """ + Gradient of the hessian (d2logpdf_dlink2) w.r.t shape parameter + + .. math:: + + :param inv_link_f: latent variables link(f) + :type inv_link_f: Nx1 array + :param y: data + :type y: Nx1 array + :param Y_metadata: includes censoring information in dictionary key 'censored' + :returns: derivative of hessian evaluated at points f and f_j w.r.t variance parameter + :rtype: Nx1 array + """ + d2logpdf_dlink_dr = -np.exp(-link_f)* (y ** self.r) * np.log(y) + return d2logpdf_dlink_dr + + def d3logpdf_dlink3_dr(self, link_f, y, Y_metadata=None): + d3logpdf_dlink_dr = np.exp(-link_f)* (y ** self.r) * np.log(y) + return d3logpdf_dlink_dr + + def dlogpdf_link_dtheta(self, f, y, Y_metadata=None): + dlogpdf_dtheta = np.zeros((self.size, f.shape[0], f.shape[1])) + dlogpdf_dtheta[0, :, :] = self.dlogpdf_link_dr(f, y, Y_metadata=Y_metadata) + return dlogpdf_dtheta + + def dlogpdf_dlink_dtheta(self, f, y, Y_metadata=None): + dlogpdf_dlink_dtheta = np.zeros((self.size, f.shape[0], f.shape[1])) + dlogpdf_dlink_dtheta[0,:,:] = self.dlogpdf_dlink_dr(f,y,Y_metadata) + return dlogpdf_dlink_dtheta + + def d2logpdf_dlink2_dtheta(self, f, y, Y_metadata=None): + d2logpdf_dlink_dtheta2 = np.zeros((self.size, f.shape[0], f.shape[1])) + d2logpdf_dlink_dtheta2[0,:,:] = self.d2logpdf_dlink2_dr(f,y,Y_metadata) + return d2logpdf_dlink_dtheta2 + + def update_gradients(self, grads): + """ + Pull out the gradients, be careful as the order must match the order + in which the parameters are added + """ + self.r.gradient = grads[0] \ No newline at end of file diff --git a/GPy/testing/likelihood_tests.py b/GPy/testing/likelihood_tests.py index 629705ab..3739995a 100644 --- a/GPy/testing/likelihood_tests.py +++ b/GPy/testing/likelihood_tests.py @@ -123,6 +123,11 @@ class TestNoiseModels(object): self.var = 0.2 self.deg_free = 4.0 + censored = np.zeros_like(self.Y) + random_inds = np.random.choice(self.N, int(self.N / 2), replace=True) + censored[random_inds] = 1 + self.Y_metadata = dict() + self.Y_metadata['censored'] = censored #Make a bigger step as lower bound can be quite curved self.step = 1e-4 @@ -274,6 +279,33 @@ class TestNoiseModels(object): "Y_metadata": {'trials': self.ns}, "laplace": True, }, + "loglogistic_censored": { + "model": GPy.likelihoods.LogLogistic(), + "link_f_constraints": [self.constrain_positive], + "Y": self.positive_Y, + "Y_metadata": self.Y_metadata, + "laplace": True + }, + "weibull_censored": { + "model": GPy.likelihoods.Weibull(), + "link_f_constraints": [self.constrain_positive], + "Y": self.positive_Y, + "Y_metadata": self.Y_metadata, + "laplace": True + }, + "loggaussian": { + "model": GPy.likelihoods.LogGaussian(), + "link_f_constraints": [self.constrain_positive], + "Y": self.positive_Y, + "laplace": True + }, + "loggaussian_censored": { + "model": GPy.likelihoods.LogGaussian(), + "link_f_constraints": [self.constrain_positive], + "Y": self.positive_Y, + "Y_metadata": self.Y_metadata, + "laplace": True + } #, #GAMMA needs some work!"Gamma_default": { #"model": GPy.likelihoods.Gamma(), From 84e39667fa1c9fce67f2854c31ceac1efa51e64f Mon Sep 17 00:00:00 2001 From: Akash Kumar Dhaka Date: Thu, 1 Jun 2017 16:07:50 +0300 Subject: [PATCH 2/3] change import statements for calling locally --- GPy/likelihoods/loglogistic.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/GPy/likelihoods/loglogistic.py b/GPy/likelihoods/loglogistic.py index 2695bd69..41e5e960 100644 --- a/GPy/likelihoods/loglogistic.py +++ b/GPy/likelihoods/loglogistic.py @@ -5,10 +5,11 @@ from __future__ import division import numpy as np from scipy import stats,special import scipy as sp -from GPy.likelihoods.link_functions import Log, Identity -from GPy.likelihoods.likelihood import Likelihood -from GPy.core.parameterization import Param -from GPy.core.parameterization.transformations import Logexp +from ..core.parameterization import Param +from ..core.parameterization.transformations import Logexp +from . import link_functions +from .likelihood import Likelihood +from .link_functions import Log class LogLogistic(Likelihood): """ From 2eef5768774b21335f214e40289d2dd20dfb33c9 Mon Sep 17 00:00:00 2001 From: Akash Kumar Dhaka Date: Thu, 27 Jul 2017 21:11:20 +0300 Subject: [PATCH 3/3] correcting weibull- changing parameterisation from f to exp(f) similar to loglogistic --- GPy/likelihoods/__init__.py | 1 - GPy/likelihoods/weibull.py | 190 ++++++++++++++++++++++++-------- GPy/testing/likelihood_tests.py | 13 --- 3 files changed, 144 insertions(+), 60 deletions(-) diff --git a/GPy/likelihoods/__init__.py b/GPy/likelihoods/__init__.py index aa9464b5..83941093 100644 --- a/GPy/likelihoods/__init__.py +++ b/GPy/likelihoods/__init__.py @@ -9,4 +9,3 @@ from .mixed_noise import MixedNoise from .binomial import Binomial from .weibull import Weibull from .loglogistic import LogLogistic -from .loggaussian import LogGaussian diff --git a/GPy/likelihoods/weibull.py b/GPy/likelihoods/weibull.py index 8ee9ba7f..ba9eb540 100644 --- a/GPy/likelihoods/weibull.py +++ b/GPy/likelihoods/weibull.py @@ -3,54 +3,50 @@ import numpy as np -from scipy import stats,special +from scipy import stats, special import scipy as sp from ..core.parameterization import Param from ..core.parameterization.transformations import Logexp from . import link_functions from .likelihood import Likelihood + class Weibull(Likelihood): """ - .. math:: - $$ p(y_{i}|f_{i}, z_{i}) = \\prod_{i=1}^{n} [ r^{1-z_{i}}\\exp(-(1-z_{i})f(x_{i}))y_{i}^{(1-z_{i})(r-1)}\\exp(-y_{i}^{r}/\\exp(f(x_{i})))] $$ + Implementing Weibull likelihood function ... - .. note: - where z_{i} is the censoring indicator- 0 for non-censored data, and 1 for censored data and r is the shape parameter. """ - def __init__(self,gp_link=None,beta=1.): + + def __init__(self, gp_link=None, beta=1.): if gp_link is None: - gp_link = link_functions.Log() + #Parameterised not as link_f but as f # gp_link = link_functions.Identity() + #Parameterised as link_f + gp_link = link_functions.Log() super(Weibull, self).__init__(gp_link, name='Weibull') self.r = Param('r_weibull_shape', float(beta), Logexp()) self.link_parameter(self.r) - # self.r.fix() - def pdf_link(self, link_f, y, Y_metadata=None): """ Likelihood function given link(f) - .. math:: - p(y_{i}|\\lambda(f_{i})) = \\frac{\\beta^{\\alpha_{i}}}{\\Gamma(\\alpha_{i})}y_{i}^{\\alpha_{i}-1}e^{-\\beta y_{i}}\\\\ - \\alpha_{i} = \\beta y_{i} :param link_f: latent variables link(f) :type link_f: Nx1 array :param y: data :type y: Nx1 array - :param Y_metadata: includes censoring information in dictionary key 'censored' + :param Y_metadata: Y_metadata which is not used in weibull distribution :returns: likelihood evaluated for this point :rtype: float """ assert np.atleast_1d(link_f).shape == np.atleast_1d(y).shape - c = np.zeros((link_f.shape[0],)) - log_objective = np.log(self.r) + (self.r - 1) * np.log(y) - link_f - (np.exp(-link_f)*(y ** self.r)) + # log_objective = np.log(self.r) + (self.r - 1) * np.log(y) - link_f - (np.exp(-link_f) * (y ** self.r)) # log_objective = stats.weibull_min.pdf(y,c=self.beta,loc=link_f,scale=1.) + log_objective = self.logpdf_link(link_f, y, Y_metadata) return np.exp(log_objective) def logpdf_link(self, link_f, y, Y_metadata=None): @@ -65,15 +61,24 @@ class Weibull(Likelihood): :type link_f: Nx1 array :param y: data :type y: Nx1 array - :param Y_metadata: includes censoring information in dictionary key 'censored' + :param Y_metadata: Y_metadata which is not used in poisson distribution :returns: likelihood evaluated for this point :rtype: float """ - #alpha = self.gp_link.transf(gp)*self.beta sum(log(a) + (a-1).*log(y)- f - exp(-f).*y.^a) - #return (1. - alpha)*np.log(obs) + self.beta*obs - alpha * np.log(self.beta) + np.log(special.gamma(alpha)) + # alpha = self.gp_link.transf(gp)*self.beta sum(log(a) + (a-1).*log(y)- f - exp(-f).*y.^a) + # return (1. - alpha)*np.log(obs) + self.beta*obs - alpha * np.log(self.beta) + np.log(special.gamma(alpha)) assert np.atleast_1d(link_f).shape == np.atleast_1d(y).shape - log_objective = np.log(self.r) + (self.r - 1) * np.log(y) - link_f - (np.exp(-link_f) * (y ** self.r)) + c = np.zeros_like(y) + if Y_metadata is not None and 'censored' in Y_metadata.keys(): + c = Y_metadata['censored'] + + # uncensored = (1-c)* (np.log(self.r) + (self.r - 1) * np.log(y) - link_f - (np.exp(-link_f) * (y ** self.r))) + # censored = (-c)*np.exp(-link_f)*(y**self.r) + uncensored = (1-c)*( np.log(self.r)-np.log(link_f)+(self.r-1)*np.log(y) - y**self.r/link_f) + censored = -c*y**self.r/link_f + + log_objective = uncensored + censored return log_objective def dlogpdf_dlink(self, link_f, y, Y_metadata=None): @@ -88,14 +93,21 @@ class Weibull(Likelihood): :type link_f: Nx1 array :param y: data :type y: Nx1 array - :param Y_metadata: includes censoring information in dictionary key 'censored' - :returns: likelihood evaluated for this point + :param Y_metadata: Y_metadata which is not used in gamma distribution + :returns: gradient of likelihood evaluated at points :rtype: Nx1 array """ # grad = (1. - self.beta) / (y - link_f) + c = np.zeros_like(y) + if Y_metadata is not None and 'censored' in Y_metadata.keys(): + c = Y_metadata['censored'] - grad = -1 + np.exp(-link_f)*(y ** self.r) + # uncensored = (1-c)* ( -1 + np.exp(-link_f)*(y ** self.r)) + # censored = c*np.exp(-link_f)*(y**self.r) + uncensored = (1-c)*(-1/link_f + y**self.r/link_f**2) + censored = c*y**self.r/link_f**2 + grad = uncensored + censored return grad def d2logpdf_dlink2(self, link_f, y, Y_metadata=None): @@ -112,8 +124,8 @@ class Weibull(Likelihood): :type link_f: Nx1 array :param y: data :type y: Nx1 array - :param Y_metadata: includes censoring information in dictionary key 'censored' - :returns: likelihood evaluated for this point + :param Y_metadata: Y_metadata which is not used in gamma distribution + :returns: Diagonal of hessian matrix (second derivative of likelihood evaluated at points f) :rtype: Nx1 array .. Note:: @@ -121,7 +133,16 @@ class Weibull(Likelihood): (the distribution for y_i depends only on link(f_i) not on link(f_(j!=i)) """ # hess = (self.beta - 1.) / (y - link_f)**2 - hess = -(y ** self.r) * np.exp(-link_f) + c = np.zeros_like(y) + if Y_metadata is not None and 'censored' in Y_metadata.keys(): + c = Y_metadata['censored'] + + # uncensored = (1-c)* (-(y ** self.r) * np.exp(-link_f)) + # censored = -c*np.exp(-link_f)*y**self.r + uncensored = (1-c)*(1/link_f**2 -2*y**self.r/link_f**3) + censored = -c*2*y**self.r/link_f**3 + hess = uncensored + censored + # hess = -(y ** self.r) * np.exp(-link_f) return hess def d3logpdf_dlink3(self, link_f, y, Y_metadata=None): @@ -141,14 +162,25 @@ class Weibull(Likelihood): :rtype: Nx1 array """ # d3lik_dlink3 = (1. - self.beta) / (y - link_f)**3 - d3lik_dlink3 = (y ** self.r) * np.exp(-link_f) + + c = np.zeros_like(y) + if Y_metadata is not None and 'censored' in Y_metadata.keys(): + c = Y_metadata['censored'] + # uncensored = (1-c)* ((y ** self.r) * np.exp(-link_f)) + # censored = c*np.exp(-link_f)*y**self.r + uncensored = (1-c)*(-2/link_f**3+ 6*y**self.r/link_f**4) + censored = c*6*y**self.r/link_f**4 + + d3lik_dlink3 = uncensored + censored + # d3lik_dlink3 = (y ** self.r) * np.exp(-link_f) return d3lik_dlink3 - def exact_inference_gradients(self, dL_dKdiag,Y_metadata=None): + def exact_inference_gradients(self, dL_dKdiag, Y_metadata=None): return np.zeros(self.size) def dlogpdf_link_dr(self, inv_link_f, y, Y_metadata=None): """ + Gradient of the log-likelihood function at y given f, w.r.t shape parameter .. math:: @@ -161,10 +193,16 @@ class Weibull(Likelihood): :returns: derivative of likelihood evaluated at points f w.r.t variance parameter :rtype: float """ - dlogpdf_dr = 1./self.r + np.log(y) - np.exp(-inv_link_f)*(y**self.r)*np.log(y) + c = np.zeros_like(y) + link_f = inv_link_f + if Y_metadata is not None and 'censored' in Y_metadata.keys(): + c = Y_metadata['censored'] + uncensored = (1-c)* (1./self.r + np.log(y) - y**self.r*np.log(y)/link_f) + censored = (-c*y**self.r*np.log(y)/link_f) + dlogpdf_dr = uncensored + censored return dlogpdf_dr - def dlogpdf_dlink_dr(self, link_f, y, Y_metadata=None): + def dlogpdf_dlink_dr(self, inv_link_f, y, Y_metadata=None): """ First order derivative derivative of loglikelihood wrt r:shape parameter @@ -177,43 +215,92 @@ class Weibull(Likelihood): :rtype: Nx1 array """ # dlogpdf_dlink_dr = self.beta * y**(self.beta - 1) * np.exp(-link_f) - dlogpdf_dlink_dr = np.exp(-link_f)* (y ** self.r) * np.log(y) + # dlogpdf_dlink_dr = np.exp(-link_f) * (y ** self.r) * np.log(y) + c = np.zeros_like(y) + if Y_metadata is not None and 'censored' in Y_metadata.keys(): + c = Y_metadata['censored'] + + link_f = inv_link_f + # uncensored = (1-c)*(np.exp(-link_f)* (y ** self.r) * np.log(y)) + # censored = c*np.exp(-link_f)*(y**self.r)*np.log(y) + uncensored = (1-c)*(y**self.r*np.log(y)/link_f**2) + censored = c*(y**self.r*np.log(y)/link_f**2) + dlogpdf_dlink_dr = uncensored + censored return dlogpdf_dlink_dr def d2logpdf_dlink2_dr(self, link_f, y, Y_metadata=None): """ - Gradient of the hessian (d2logpdf_dlink2) w.r.t shape parameter - .. math:: - - :param inv_link_f: latent variables link(f) - :type inv_link_f: Nx1 array - :param y: data - :type y: Nx1 array - :param Y_metadata: includes censoring information in dictionary key 'censored' - :returns: derivative of hessian evaluated at points f and f_j w.r.t variance parameter - :rtype: Nx1 array + Derivative of hessian of loglikelihood wrt r-shape parameter. + :param link_f: + :param y: + :param Y_metadata: + :return: """ - d2logpdf_dlink_dr = -np.exp(-link_f)* (y ** self.r) * np.log(y) + + c = np.zeros_like(y) + if Y_metadata is not None and 'censored' in Y_metadata.keys(): + c = Y_metadata['censored'] + + # uncensored = (1-c)*( -np.exp(-link_f)* (y ** self.r) * np.log(y)) + # censored = -c*np.exp(-link_f)*(y**self.r)*np.log(y) + uncensored = (1-c)*-2*y**self.r*np.log(y)/link_f**3 + censored = c*-2*y**self.r*np.log(y)/link_f**3 + d2logpdf_dlink_dr = uncensored + censored + return d2logpdf_dlink_dr def d3logpdf_dlink3_dr(self, link_f, y, Y_metadata=None): - d3logpdf_dlink_dr = np.exp(-link_f)* (y ** self.r) * np.log(y) - return d3logpdf_dlink_dr + """ + + :param link_f: + :param y: + :param Y_metadata: + :return: + """ + c = np.zeros_like(y) + if Y_metadata is not None and 'censored' in Y_metadata.keys(): + c = Y_metadata['censored'] + + uncensored = (1-c)* ((y**self.r)*np.exp(-link_f)*np.log1p(y)) + censored = c*np.exp(-link_f)*(y**self.r)*np.log(y) + d3logpdf_dlink3_dr = uncensored + censored + return d3logpdf_dlink3_dr def dlogpdf_link_dtheta(self, f, y, Y_metadata=None): + """ + + :param f: + :param y: + :param Y_metadata: + :return: + """ dlogpdf_dtheta = np.zeros((self.size, f.shape[0], f.shape[1])) dlogpdf_dtheta[0, :, :] = self.dlogpdf_link_dr(f, y, Y_metadata=Y_metadata) return dlogpdf_dtheta def dlogpdf_dlink_dtheta(self, f, y, Y_metadata=None): + """ + + :param f: + :param y: + :param Y_metadata: + :return: + """ dlogpdf_dlink_dtheta = np.zeros((self.size, f.shape[0], f.shape[1])) - dlogpdf_dlink_dtheta[0,:,:] = self.dlogpdf_dlink_dr(f,y,Y_metadata) + dlogpdf_dlink_dtheta[0, :, :] = self.dlogpdf_dlink_dr(f, y, Y_metadata) return dlogpdf_dlink_dtheta def d2logpdf_dlink2_dtheta(self, f, y, Y_metadata=None): + """ + + :param f: + :param y: + :param Y_metadata: + :return: + """ d2logpdf_dlink_dtheta2 = np.zeros((self.size, f.shape[0], f.shape[1])) - d2logpdf_dlink_dtheta2[0,:,:] = self.d2logpdf_dlink2_dr(f,y,Y_metadata) + d2logpdf_dlink_dtheta2[0, :, :] = self.d2logpdf_dlink2_dr(f, y, Y_metadata) return d2logpdf_dlink_dtheta2 def update_gradients(self, grads): @@ -221,4 +308,15 @@ class Weibull(Likelihood): Pull out the gradients, be careful as the order must match the order in which the parameters are added """ - self.r.gradient = grads[0] \ No newline at end of file + self.r.gradient = grads[0] + + def samples(self, gp, Y_metadata=None): + """ + Returns a set of samples of observations conditioned on a given value of latent variable f. + + :param gp: latent variable + """ + orig_shape = gp.shape + gp = gp.flatten() + weibull_samples = np.array([sp.stats.weibull_min.rvs(self.r, loc=0, scale=self.gp_link.transf(f)) for f in gp]) + return weibull_samples.reshape(orig_shape) \ No newline at end of file diff --git a/GPy/testing/likelihood_tests.py b/GPy/testing/likelihood_tests.py index 3739995a..ac681ecc 100644 --- a/GPy/testing/likelihood_tests.py +++ b/GPy/testing/likelihood_tests.py @@ -292,19 +292,6 @@ class TestNoiseModels(object): "Y": self.positive_Y, "Y_metadata": self.Y_metadata, "laplace": True - }, - "loggaussian": { - "model": GPy.likelihoods.LogGaussian(), - "link_f_constraints": [self.constrain_positive], - "Y": self.positive_Y, - "laplace": True - }, - "loggaussian_censored": { - "model": GPy.likelihoods.LogGaussian(), - "link_f_constraints": [self.constrain_positive], - "Y": self.positive_Y, - "Y_metadata": self.Y_metadata, - "laplace": True } #, #GAMMA needs some work!"Gamma_default": {