from __future__ import absolute_import
# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
# vi: set ft=python sts=4 ts=4 sw=4 et:
import numpy as np
import scipy.stats
[docs]class Link(object):
"""
A generic link function for one-parameter exponential
family, with call, inverse and deriv methods.
"""
[docs] def initialize(self, Y):
return np.asarray(Y).mean() * np.ones(Y.shape)
def __call__(self, p):
return NotImplementedError
[docs] def inverse(self, z):
return NotImplementedError
[docs] def deriv(self, p):
return NotImplementedError
[docs]class Logit(Link):
"""
The logit transform as a link function:
g'(x) = 1 / (x * (1 - x))
g^(-1)(x) = exp(x)/(1 + exp(x))
"""
tol = 1.0e-10
[docs] def clean(self, p):
"""
Clip logistic values to range (tol, 1-tol)
INPUTS:
p -- probabilities
OUTPUTS: pclip
pclip -- clipped probabilities
"""
return np.clip(p, Logit.tol, 1. - Logit.tol)
def __call__(self, p):
"""
Logit transform
g(p) = log(p / (1 - p))
INPUTS:
p -- probabilities
OUTPUTS: z
z -- logit transform of p
"""
p = self.clean(p)
return np.log(p / (1. - p))
[docs] def inverse(self, z):
"""
Inverse logit transform
h(z) = exp(z)/(1+exp(z))
INPUTS:
z -- logit transform of p
OUTPUTS: p
p -- probabilities
"""
t = np.exp(z)
return t / (1. + t)
[docs] def deriv(self, p):
"""
Derivative of logit transform
g(p) = 1 / (p * (1 - p))
INPUTS:
p -- probabilities
OUTPUTS: y
y -- derivative of logit transform of p
"""
p = self.clean(p)
return 1. / (p * (1 - p))
logit = Logit()
[docs]class Power(Link):
"""
The power transform as a link function:
g(x) = x**power
"""
[docs] def __init__(self, power=1.):
self.power = power
def __call__(self, x):
"""
Power transform
g(x) = x**self.power
INPUTS:
x -- mean parameters
OUTPUTS: z
z -- power transform of x
"""
return np.power(x, self.power)
[docs] def inverse(self, z):
"""
Inverse of power transform
g(x) = x**(1/self.power)
INPUTS:
z -- linear predictors in GLM
OUTPUTS: x
x -- mean parameters
"""
return np.power(z, 1. / self.power)
[docs] def deriv(self, x):
"""
Derivative of power transform
g(x) = self.power * x**(self.power - 1)
INPUTS:
x -- mean parameters
OUTPUTS: z
z -- derivative of power transform of x
"""
return self.power * np.power(x, self.power - 1)
inverse = Power(power=-1.)
inverse.__doc__ = """
The inverse transform as a link function:
g(x) = 1 / x
"""
sqrt = Power(power=0.5)
sqrt.__doc__ = """
The square-root transform as a link function:
g(x) = sqrt(x)
"""
inverse_squared = Power(power=-2.)
inverse_squared.__doc__ = """
The inverse squared transform as a link function:
g(x) = 1 / x**2
"""
identity = Power(power=1.)
identity.__doc__ = """
The identity transform as a link function:
g(x) = x
"""
[docs]class Log(Link):
"""
The log transform as a link function:
g(x) = log(x)
"""
tol = 1.0e-10
[docs] def clean(self, x):
return np.clip(x, Logit.tol, np.inf)
def __call__(self, x, **extra):
"""
Log transform
g(x) = log(x)
INPUTS:
x -- mean parameters
OUTPUTS: z
z -- log(x)
"""
x = self.clean(x)
return np.log(x)
[docs] def inverse(self, z):
"""
Inverse of log transform
g(x) = exp(x)
INPUTS:
z -- linear predictors in GLM
OUTPUTS: x
x -- exp(z)
"""
return np.exp(z)
[docs] def deriv(self, x):
"""
Derivative of log transform
g(x) = 1/x
INPUTS:
x -- mean parameters
OUTPUTS: z
z -- derivative of log transform of x
"""
x = self.clean(x)
return 1. / x
log = Log()
[docs]class CDFLink(Logit):
"""
The use the CDF of a scipy.stats distribution as a link function:
g(x) = dbn.ppf(x)
"""
[docs] def __init__(self, dbn=scipy.stats.norm):
self.dbn = dbn
def __call__(self, p):
"""
CDF link
g(p) = self.dbn.pdf(p)
INPUTS:
p -- mean parameters
OUTPUTS: z
z -- derivative of CDF transform of p
"""
p = self.clean(p)
return self.dbn.ppf(p)
[docs] def inverse(self, z):
"""
Derivative of CDF link
g(z) = self.dbn.cdf(z)
INPUTS:
z -- linear predictors in GLM
OUTPUTS: p
p -- inverse of CDF link of z
"""
return self.dbn.cdf(z)
[docs] def deriv(self, p):
"""
Derivative of CDF link
g(p) = 1/self.dbn.pdf(self.dbn.ppf(p))
INPUTS:
x -- mean parameters
OUTPUTS: z
z -- derivative of CDF transform of x
"""
p = self.clean(p)
return 1. / self.dbn.pdf(self(p))
probit = CDFLink()
probit.__doc__ = """
The probit (standard normal CDF) transform as a link function:
g(x) = scipy.stats.norm.ppf(x)
"""
cauchy = CDFLink(dbn=scipy.stats.cauchy)
cauchy.__doc__ = """
The Cauchy (standard Cauchy CDF) transform as a link function:
g(x) = scipy.stats.cauchy.ppf(x)
"""
[docs]class CLogLog(Logit):
"""
The complementary log-log transform as a link function:
g(x) = log(-log(x))
"""
def __call__(self, p):
"""
C-Log-Log transform
g(p) = log(-log(p))
INPUTS:
p -- mean parameters
OUTPUTS: z
z -- log(-log(p))
"""
p = self.clean(p)
return np.log(-np.log(p))
[docs] def inverse(self, z):
"""
Inverse of C-Log-Log transform
g(z) = exp(-exp(z))
INPUTS:
z -- linear predictor scale
OUTPUTS: p
p -- mean parameters
"""
return np.exp(-np.exp(z))
[docs] def deriv(self, p):
"""
Derivatve of C-Log-Log transform
g(p) = - 1 / (log(p) * p)
INPUTS:
p -- mean parameters
OUTPUTS: z
z -- - 1 / (log(p) * p)
"""
p = self.clean(p)
return -1. / (np.log(p) * p)
cloglog = CLogLog()