# Source code for GPy.likelihoods.loggaussian

# Copyright (c) 2012 - 2014, GPy authors (see AUTHORS.txt).

import numpy as np
from scipy import stats, special
from ..core.parameterization import Param
from ..core.parameterization.transformations import Logexp
from .likelihood import Likelihood

[docs]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.

"""

self.sigma = Param('sigma', sigma, Logexp())
self.variance = Param('variance', sigma**2, Logexp())

"""
: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
"""

"""
: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_like(y)

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

"""
derivative of logpdf wrt link_f param
.. math::

: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_like(y)

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

"""
The hessian will be 0 unless i == j

.. math::

: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
"""
# c = np.zeros((y.shape[0],))
c = np.zeros_like(y)

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

"""
Gradient of the log-likelihood function at y given f, w.r.t shape parameter

.. math::

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

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

"""
Gradient of the log-likelihood function at y given f, w.r.t variance parameter

.. math::

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

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

"""
:param y: data
:type y: Nx1 array
:returns: derivative of log likelihood evaluated at points link(f) w.r.t variance parameter
:rtype: Nx1 array
"""
c = np.zeros_like(y)

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

"""
:param y: data
:type y: Nx1 array
:returns: derivative of log likelihood evaluated at points link(f) w.r.t variance parameter
:rtype: Nx1 array
"""
c = np.zeros_like(y)
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

"""
:param y: data
:type y: Nx1 array
: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]))
return dlogpdf_dtheta

"""
:param y: data
:type y: Nx1 array
:returns: derivative of log likelihood evaluated at points link(f) w.r.t variance parameter
:rtype: Nx1 array
"""

"""
:param y: data
:type y: Nx1 array
:returns: derivative of log likelihood evaluated at points link(f) w.r.t variance parameter
:rtype: Nx1 array
"""

"""
Pull out the gradients, be careful as the order must match the order
in which the parameters are added
"""