"""Likelihood objects for inference within the GLM framework."""
from __future__ import division
import numpy as np
from scipy.stats import bernoulli, binom, poisson, norm, gamma
from scipy.special import gammaln, expit
from .btypes import Parameter, Positive
from .mathfun.special import safesoftplus, softplus
#
# Likelihood classes
#
class Bernoulli():
r"""
Bernoulli likelihood class for (binary) classification tasks.
A logistic transformation function is used to map the latent function from
the GLM prior into a probability.
.. math::
p(y_i | f_i) = \sigma(f_i) ^ {y_i} (1 - \sigma(f_i))^{1 - y_i}
where :math:`y_i` is a target, :math:`f_i` the value of the latent function
corresponding to the target, and :math:`\sigma(\cdot)` is the logistic
sigmoid.
"""
_params = Parameter()
@property
def params(self):
"""Get this object's Parameter types."""
return self._params
@params.setter
def params(self, params):
"""Set this object's Parameter types."""
self._params = params
def loglike(self, y, f):
r"""
Bernoulli log likelihood.
Parameters
----------
y: ndarray
array of 0, 1 valued integers of targets
f: ndarray
latent function from the GLM prior (:math:`\mathbf{f} =
\boldsymbol\Phi \mathbf{w}`)
Returns
-------
logp: ndarray
the log likelihood of each y given each f under this
likelihood.
"""
# way faster than calling bernoulli.logpmf
y, f = np.broadcast_arrays(y, f)
ll = y * f - softplus(f)
return ll
def Ey(self, f):
r"""
Expected value of the Bernoulli likelihood.
Parameters
----------
f: ndarray
latent function from the GLM prior (:math:`\mathbf{f} =
\boldsymbol\Phi \mathbf{w}`)
Returns
-------
Ey: ndarray
expected value of y, :math:`\mathbb{E}[\mathbf{y}|\mathbf{f}]`.
"""
return expit(f)
def df(self, y, f):
r"""
Derivative of Bernoulli log likelihood w.r.t.\ f.
Parameters
----------
y: ndarray
array of 0, 1 valued integers of targets
f: ndarray
latent function from the GLM prior (:math:`\mathbf{f} =
\boldsymbol\Phi \mathbf{w}`)
Returns
-------
df: ndarray
the derivative :math:`\partial \log p(y|f) / \partial f`
"""
y, f = np.broadcast_arrays(y, f)
return y - expit(f)
def dp(self, y, f, *args):
r"""
Derivative of Bernoulli log likelihood w.r.t.\ the parameters,
:math:`\theta`.
Parameters
----------
y: ndarray
array of 0, 1 valued integers of targets
f: ndarray
latent function from the GLM prior (:math:`\mathbf{f} =
\boldsymbol\Phi \mathbf{w}`)
Returns
-------
dp: list, float or ndarray
the derivative
:math:`\partial \log p(y|f, \theta)/ \partial \theta` for
each parameter. If there is only one parameter, this is not a
list.
"""
return []
def cdf(self, y, f):
r"""
Cumulative density function of the likelihood.
Parameters
----------
y: ndarray
query quantiles, i.e.\ :math:`P(Y \leq y)`.
f: ndarray
latent function from the GLM prior (:math:`\mathbf{f} =
\boldsymbol\Phi \mathbf{w}`)
Returns
-------
cdf: ndarray
Cumulative density function evaluated at y.
"""
return bernoulli.cdf(y, expit(f))
def __repr__(self):
"""Representation."""
return "{}()".format(self.__class__.__name__)
class Binomial(Bernoulli):
r"""
Binomial likelihood class.
A logistic transformation function is used to map the latent function from
the GLM prior into a probability.
.. math::
p(y_i | f_i) = \genfrac(){0pt}{}{n}{y_i}
\sigma(f_i) ^ {y_i} (1 - \sigma(f_i))^{n - y_i}
where :math:`y_i` is a target, :math:`f_i` the value of the latent function
corresponding to the target, :math:`n` is the total possible count, and
:math:`\sigma(\cdot)` is the logistic sigmoid. :math:`n` can also be
applied per observation.
"""
def loglike(self, y, f, n):
r"""
Binomial log likelihood.
Parameters
----------
y: ndarray
array of 0, 1 valued integers of targets
f: ndarray
latent function from the GLM prior (:math:`\mathbf{f} =
\boldsymbol\Phi \mathbf{w}`)
n: ndarray
the total number of observations
Returns
-------
logp: ndarray
the log likelihood of each y given each f under this
likelihood.
"""
ll = binom.logpmf(y, n=n, p=expit(f))
return ll
def Ey(self, f, n):
r"""
Expected value of the Binomial likelihood.
Parameters
----------
f: ndarray
latent function from the GLM prior (:math:`\mathbf{f} =
\boldsymbol\Phi \mathbf{w}`)
n: ndarray
the total number of observations
Returns
-------
Ey: ndarray
expected value of y, :math:`\mathbb{E}[\mathbf{y}|\mathbf{f}]`.
"""
return expit(f) * n
def df(self, y, f, n):
r"""
Derivative of Binomial log likelihood w.r.t.\ f.
Parameters
----------
y: ndarray
array of 0, 1 valued integers of targets
f: ndarray
latent function from the GLM prior (:math:`\mathbf{f} =
\boldsymbol\Phi \mathbf{w}`)
n: ndarray
the total number of observations
Returns
-------
df: ndarray
the derivative :math:`\partial \log p(y|f) / \partial f`
"""
y, f, n = np.broadcast_arrays(y, f, n)
return y - expit(f) * n
def cdf(self, y, f, n):
r"""
Cumulative density function of the likelihood.
Parameters
----------
y: ndarray
query quantiles, i.e.\ :math:`P(Y \leq y)`.
f: ndarray
latent function from the GLM prior (:math:`\mathbf{f} =
\boldsymbol\Phi \mathbf{w}`)
n: ndarray
the total number of observations
Returns
-------
cdf: ndarray
Cumulative density function evaluated at y.
"""
return binom.cdf(y, n=n, p=expit(f))
def __repr__(self):
"""Representation."""
return "{}()".format(self.__class__.__name__)
class Gaussian(Bernoulli):
r"""
A univariate Gaussian likelihood for general regression tasks.
No transformation function is needed since this is (conditionally)
conjugate to the GLM prior.
.. math::
p(y_i | f_i) = \frac{1}{\sqrt{2 \pi \sigma^2}}
\exp\left(- \frac{(y_i - f_i)^2}{2 \sigma^2} \right)
where :math:`y_i` is a target, :math:`f_i` the value of the latent function
corresponding to the target and :math:`\sigma` is the observation noise
(standard deviation).
Parameters
----------
var: Parameter, optional
A scalar Parameter describing the initial point and bounds for
an optimiser to learn the variance parameter of this object.
"""
def __init__(self, var=Parameter(gamma(1., scale=1), Positive())):
"""See class docstring."""
self.params = var
def _check_param(self, param):
if param is None:
return self.params.value
if not self.params.bounds.check(param):
raise ValueError("Input parameter is out of bounds!")
return param
def loglike(self, y, f, var=None):
r"""
Gaussian log likelihood.
Parameters
----------
y: ndarray
array of 0, 1 valued integers of targets
f: ndarray
latent function from the GLM prior (:math:`\mathbf{f} =
\boldsymbol\Phi \mathbf{w}`)
var: float, ndarray, optional
The variance of the distribution, if not input, the initial value
of variance is used.
Returns
-------
logp: ndarray
the log likelihood of each y given each f under this
likelihood.
"""
# way faster than calling norm.logpdf
var = self._check_param(var)
y, f = np.broadcast_arrays(y, f)
ll = - 0.5 * (np.log(2 * np.pi * var) + (y - f)**2 / var)
return ll
def Ey(self, f, var):
r"""
Expected value of the Gaussian likelihood.
Parameters
----------
f: ndarray
latent function from the GLM prior (:math:`\mathbf{f} =
\boldsymbol\Phi \mathbf{w}`)
var: float, ndarray, optional
The variance of the distribution, if not input, the initial value
of variance is used.
Returns
-------
Ey: ndarray
expected value of y, :math:`\mathbb{E}[\mathbf{y}|\mathbf{f}]`.
"""
var = self._check_param(var)
return f
def df(self, y, f, var):
r"""
Derivative of Gaussian log likelihood w.r.t.\ f.
Parameters
----------
y: ndarray
array of 0, 1 valued integers of targets
f: ndarray
latent function from the GLM prior (:math:`\mathbf{f} =
\boldsymbol\Phi \mathbf{w}`)
var: float, ndarray, optional
The variance of the distribution, if not input, the initial value
of variance is used.
Returns
-------
df: ndarray
the derivative :math:`\partial \log p(y|f) / \partial f`
"""
var = self._check_param(var)
y, f = np.broadcast_arrays(y, f)
return (y - f) / var
def dp(self, y, f, var):
r"""
Derivative of Gaussian log likelihood w.r.t.\ the variance
:math:`\sigma^2`.
Parameters
----------
y: ndarray
array of 0, 1 valued integers of targets
f: ndarray
latent function from the GLM prior (:math:`\mathbf{f} =
\boldsymbol\Phi \mathbf{w}`)
var: float, ndarray, optional
The variance of the distribution, if not input, the initial value
of variance is used.
Returns
-------
dp: float
the derivative
:math:`\partial \log p(y|f, \sigma^2)/ \partial \sigma^2`
where :math:`sigma^2` is the variance.
"""
var = self._check_param(var)
y, f = np.broadcast_arrays(y, f)
ivar = 1. / var
return 0.5 * (((y - f) * ivar)**2 - ivar)
def cdf(self, y, f, var):
r"""
Cumulative density function of the likelihood.
Parameters
----------
y: ndarray
query quantiles, i.e.\ :math:`P(Y \leq y)`.
f: ndarray
latent function from the GLM prior (:math:`\mathbf{f} =
\boldsymbol\Phi \mathbf{w}`)
var: float, ndarray, optional
The variance of the distribution, if not input, the initial value
of variance is used.
Returns
-------
cdf: ndarray
Cumulative density function evaluated at y.
"""
var = self._check_param(var)
return norm.cdf(y, loc=f, scale=np.sqrt(var))
def __repr__(self):
"""Representation."""
return "{}(var={})".format(self.__class__.__name__, self.params)
class Poisson(Bernoulli):
r"""
A Poisson likelihood, useful for various Poisson process tasks.
An exponential transformation function and a softplus transformation
function have been implemented.
.. math::
p(y_i | f_i) = \frac{g(f_i)^{y_i} e^{-g(f_i)}}{y_i!}
where :math:`y_i` is a target, :math:`f_i` the value of the latent function
corresponding to the target, and :math:`g(\cdot)` is the tranformation
function, which can be either an exponential function, or a softplus
function (:math:`\log(1 + \exp(f_i)`).
Parameters
----------
tranfcn: string, optional
this may be 'exp' for an exponential transformation function,
or 'softplus' for a softplut transformation function.
"""
def __init__(self, tranfcn='exp'):
"""See class docstring."""
if tranfcn == 'exp' or tranfcn == 'softplus':
self.tranfcn = tranfcn
else:
raise ValueError('Invalid transformation function specified!')
def loglike(self, y, f):
r"""
Poisson log likelihood.
Parameters
----------
y: ndarray
array of integer targets
f: ndarray
latent function from the GLM prior (:math:`\mathbf{f} =
\boldsymbol\Phi \mathbf{w}`)
Returns
-------
logp: ndarray
the log likelihood of each y given each f under this
likelihood.
"""
y, f = np.broadcast_arrays(y, f)
if self.tranfcn == 'exp':
g = np.exp(f)
logg = f
else:
g = softplus(f)
logg = np.log(g)
return y * logg - g - gammaln(y + 1)
def Ey(self, f):
r"""
Expected value of the Poisson likelihood.
Parameters
----------
f: ndarray
latent function from the GLM prior (:math:`\mathbf{f} =
\boldsymbol\Phi \mathbf{w}`)
Returns
-------
Ey: ndarray
expected value of y, :math:`\mathbb{E}[\mathbf{y}|\mathbf{f}]`.
"""
return np.exp(f) if self.tranfcn == 'exp' else softplus(f)
def df(self, y, f):
r"""
Derivative of Poisson log likelihood w.r.t.\ f.
Parameters
----------
y: ndarray
array of 0, 1 valued integers of targets
f: ndarray
latent function from the GLM prior (:math:`\mathbf{f} =
\boldsymbol\Phi \mathbf{w}`)
Returns
-------
df: ndarray
the derivative :math:`\partial \log p(y|f) / \partial f`
"""
y, f = np.broadcast_arrays(y, f)
if self.tranfcn == 'exp':
return y - np.exp(f)
else:
return expit(f) * (y / safesoftplus(f) - 1)
def cdf(self, y, f):
r"""
Cumulative density function of the likelihood.
Parameters
----------
y: ndarray
query quantiles, i.e.\ :math:`P(Y \leq y)`.
f: ndarray
latent function from the GLM prior (:math:`\mathbf{f} =
\boldsymbol\Phi \mathbf{w}`)
Returns
-------
cdf: ndarray
Cumulative density function evaluated at y.
"""
mu = np.exp(f) if self.tranfcn == 'exp' else softplus(f)
return poisson.cdf(y, mu=mu)
def __repr__(self):
"""Representation."""
return "{}(tranfcn='{}')".format(self.__class__.__name__, self.tranfcn)
