Merge pull request #527 from adhaka/ll-surv

merged in the changes about Ll surv

GPy/likelihoods/__init__.py

@@ -7,4 +7,5 @@ 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

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()

GPy/likelihoods/loglogistic.py

@@ -0,0 +1,339 @@
+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 ..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):
+    """
+    .. 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)
+

GPy/likelihoods/weibull.py

@@ -0,0 +1,322 @@
+# 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):
+    """
+    Implementing Weibull likelihood function ...
+
+    """
+
+    def __init__(self, gp_link=None, beta=1.):
+        if gp_link is None:
+            #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)
+
+    def pdf_link(self, link_f, y, Y_metadata=None):
+        """
+        Likelihood function given link(f)
+
+
+        :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 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 = 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):
+        """
+        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: 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))
+        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)* (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):
+        """
+        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: 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']
+
+        # 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):
+        """
+        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: 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::
+            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
+        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):
+        """
+        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
+
+        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):
+        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
+        """
+        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, inv_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)
+        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):
+        """
+
+        Derivative of hessian of loglikelihood wrt r-shape parameter.
+        :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)*( -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):
+        """
+
+        :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)
+        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)
+        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]
+
+    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)

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,20 @@ 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
+            }
             #,
             #GAMMA needs some work!"Gamma_default": {
             #"model": GPy.likelihoods.Gamma(),