import numpy as np
from . import _tables
def _cabs(x):
"""
Absolute value function that changes complex sign based on real sign
This could be useful for complex step derivatives of functions that
need abs. Not yet used.
Parameters
----------
x : ndarray
Input array, real or complex.
Returns
-------
ndarray
The absolute value of `x`, with the sign of the real part of `x`
applied for complex input.
"""
sign = (x.real >= 0) * 2 - 1
return sign * x
[docs]
class RobustNorm:
"""
The parent class for the norms used for robust regression
Lays out the methods expected of the robust norms to be used
by statsmodels.RLM.
See Also
--------
statsmodels.robust.robust_linear_model.RLM
Notes
-----
Currently only M-estimators are available.
References
----------
PJ Huber. 'Robust Statistics' John Wiley and Sons, Inc., New York, 1981.
DC Montgomery, EA Peck. 'Introduction to Linear Regression Analysis',
John Wiley and Sons, Inc., New York, 2001.
R Venables, B Ripley. 'Modern Applied Statistics in S'
Springer, New York, 2002.
"""
continuous = 1
def __repr__(self):
return self.__class__.__name__
[docs]
def rho(self, z):
"""
The robust criterion estimator function
Abstract method:
-2 loglike used in M-estimator
Parameters
----------
z : array_like
1d array
Returns
-------
ndarray
The value of the robust criterion function.
"""
raise NotImplementedError
[docs]
def psi(self, z):
"""
Derivative of rho. Sometimes referred to as the influence function
Abstract method:
psi = rho'
Parameters
----------
z : array_like
1d array
Returns
-------
ndarray
The value of the psi function.
"""
raise NotImplementedError
[docs]
def weights(self, z):
"""
Return the value of psi(z) / z
Abstract method:
psi(z) / z
Parameters
----------
z : array_like
1d array
Returns
-------
ndarray
The value of the weighting function.
"""
raise NotImplementedError
[docs]
def psi_deriv(self, z):
"""
Derivative of psi. Used to obtain robust covariance matrix
See statsmodels.robust.robust_linear_model.RLM for more information.
Abstract method:
psi_derive = psi'
Parameters
----------
z : array_like
1d array
Returns
-------
ndarray
The value of the derivative of the psi function.
"""
raise NotImplementedError
def __call__(self, z):
"""
Return the value of estimator rho applied to an input
Parameters
----------
z : array_like
1d array
Returns
-------
ndarray
The value of the robust criterion function.
"""
return self.rho(z)
[docs]
class LeastSquares(RobustNorm):
"""
Least squares rho for M-estimation and its derived functions
See Also
--------
statsmodels.robust.norms.RobustNorm
"""
continuous = 2
redescending = "not"
[docs]
def max_rho(self):
return np.inf
[docs]
def rho(self, z):
r"""
The robust criterion function for the least squares estimator
.. math::
\rho(z) = \frac{1}{2} z^2
Parameters
----------
z : array_like
1d array
Returns
-------
rho : ndarray
The value of the robust criterion function.
"""
return z**2 * 0.5
[docs]
def psi(self, z):
r"""
The psi function for the least squares estimator
The analytic derivative of rho.
.. math::
\psi(z) = z
Parameters
----------
z : array_like
1d array
Returns
-------
psi : ndarray
The value of the psi function.
"""
return np.asarray(z)
[docs]
def weights(self, z):
r"""
The least squares estimator weighting function for the IRLS algorithm
The psi function scaled by the input z.
.. math::
w(z) = 1
Parameters
----------
z : array_like
1d array
Returns
-------
weights : ndarray
The value of the weighting function.
"""
z = np.asarray(z)
return np.ones(z.shape, np.float64)
[docs]
def psi_deriv(self, z):
r"""
The derivative of the least squares psi function
.. math::
\psi'(z) = 1
Parameters
----------
z : array_like
1d array
Returns
-------
psi_deriv : ndarray
The value of the derivative of the psi function.
Notes
-----
Used to estimate the robust covariance matrix.
"""
z = np.asarray(z)
return np.ones(z.shape, np.float64)
[docs]
class HuberT(RobustNorm):
"""
Huber's T for M estimation
Parameters
----------
t : float, optional
The tuning constant for Huber's t function. The default value is
1.345.
See Also
--------
statsmodels.robust.norms.RobustNorm
"""
continuous = 1
redescending = "not"
def __init__(self, t=1.345):
self.t = t
def _set_tuning_param(self, c, inplace=False):
"""
Set and change the tuning parameter of the Norm
Parameters
----------
c : float
The new value for the tuning parameter `t`.
inplace : bool, optional
If True, modify this instance in place and return it. If
False, return a new instance with the new tuning parameter.
The default is False.
Returns
-------
RobustNorm
The norm instance with the new tuning parameter.
Warnings
--------
This needs to wipe cached attributes that depend on the parameter.
"""
if inplace:
self.t = c
return self
else:
return self.__class__(t=c)
[docs]
def max_rho(self):
return np.inf
def _subset(self, z):
"""
Huber's T is defined piecewise over the range of z
Parameters
----------
z : array_like
1d array
Returns
-------
ndarray
Boolean array indicating which elements of `z` fall in the
quadratic (rather than linear) region of rho.
"""
z = np.asarray(z)
return np.less_equal(np.abs(z), self.t)
[docs]
def rho(self, z):
r"""
The robust criterion function for Huber's t estimator
.. math::
\rho(z) = \begin{cases}
\frac{1}{2} z^2 & \text{for } |z| \le t \\
t|z| - \frac{1}{2} t^2 & \text{for } |z| > t
\end{cases}
Parameters
----------
z : array_like
1d array
Returns
-------
rho : ndarray
The value of the robust criterion function.
"""
z = np.asarray(z)
test = self._subset(z)
return (test * 0.5 * z**2 +
(1 - test) * (np.abs(z) * self.t - 0.5 * self.t**2))
[docs]
def psi(self, z):
r"""
The psi function for Huber's t estimator
The analytic derivative of rho.
.. math::
\psi(z) = \begin{cases}
z & \text{for } |z| \le t \\
\text{sign}(z) \cdot t & \text{for } |z| > t
\end{cases}
Parameters
----------
z : array_like
1d array
Returns
-------
psi : ndarray
The value of the psi function.
"""
z = np.asarray(z)
test = self._subset(z)
return test * z + (1 - test) * self.t * np.sign(z)
[docs]
def weights(self, z):
r"""
Huber's t weighting function for the IRLS algorithm
The psi function scaled by z.
.. math::
w(z) = \begin{cases}
1 & \text{for } |z| \le t \\
t/|z| & \text{for } |z| > t
\end{cases}
Parameters
----------
z : array_like
1d array
Returns
-------
weights : ndarray
The value of the weighting function.
"""
z_isscalar = np.isscalar(z)
z = np.atleast_1d(z)
test = self._subset(z)
absz = np.abs(z)
absz[test] = 1.0
v = test + (1 - test) * self.t / absz
if z_isscalar:
v = v[0]
return v
[docs]
def psi_deriv(self, z):
"""
Derivative of the Huber T psi function.
Parameters
----------
z : array_like
Input residual values.
Returns
-------
ndarray
The value of the derivative of the psi function. Indicator values
equal to 1 when '|z| <= t' and 0 otherwise.
Notes
-----
Used when computing the robust covariance matrix in robustlinear models.
"""
return np.less_equal(np.abs(z), self.t).astype(float)
[docs]
class RamsayE(RobustNorm):
"""
Ramsay's Ea for M estimation
Parameters
----------
a : float, optional
The tuning constant for Ramsay's Ea function. The default value is
0.3.
See Also
--------
statsmodels.robust.norms.RobustNorm
"""
continuous = 2
redescending = "soft"
def __init__(self, a=.3):
self.a = a
def _set_tuning_param(self, c, inplace=False):
"""
Set and change the tuning parameter of the Norm
Parameters
----------
c : float
The new value for the tuning parameter `a`.
inplace : bool, optional
If True, modify this instance in place and return it. If
False, return a new instance with the new tuning parameter.
The default is False.
Returns
-------
RobustNorm
The norm instance with the new tuning parameter.
Warnings
--------
This needs to wipe cached attributes that depend on the parameter.
"""
# todo : change default to inplace=False, when tools are fixed
if inplace:
self.a = c
return self
else:
return self.__class__(a=c)
[docs]
def max_rho(self):
return np.inf
[docs]
def rho(self, z):
r"""
The robust criterion function for Ramsay's Ea
.. math::
\rho(z) = a^{-2} (1 - \exp(-a|z|)(1 + a|z|))
Parameters
----------
z : array_like
1d array
Returns
-------
rho : ndarray
The value of the robust criterion function.
"""
z = np.asarray(z)
return (1 - np.exp(-self.a * np.abs(z)) *
(1 + self.a * np.abs(z))) / self.a**2
[docs]
def psi(self, z):
r"""
The psi function for Ramsay's Ea estimator
The analytic derivative of rho.
.. math::
\psi(z) = z \exp(-a|z|)
Parameters
----------
z : array_like
1d array
Returns
-------
psi : ndarray
The value of the psi function.
"""
z = np.asarray(z)
return z * np.exp(-self.a * np.abs(z))
[docs]
def weights(self, z):
r"""
Ramsay's Ea weighting function for the IRLS algorithm
The psi function scaled by z.
.. math::
w(z) = \exp(-a|z|)
Parameters
----------
z : array_like
1d array
Returns
-------
weights : ndarray
The value of the weighting function.
"""
z = np.asarray(z)
return np.exp(-self.a * np.abs(z))
[docs]
def psi_deriv(self, z):
"""
The derivative of Ramsay's Ea psi function
Parameters
----------
z : array_like
1d array
Returns
-------
ndarray
The value of the derivative of the psi function.
Notes
-----
Used to estimate the robust covariance matrix.
"""
a = self.a
x = np.exp(-a * np.abs(z))
dx = -a * x * np.sign(z)
y = z
dy = 1
return x * dy + y * dx
[docs]
class AndrewWave(RobustNorm):
"""
Andrew's wave for M estimation
Parameters
----------
a : float, optional
The tuning constant for Andrew's Wave function. The default value is
1.339.
See Also
--------
statsmodels.robust.norms.RobustNorm
"""
continuous = 1
redescending = "hard"
def __init__(self, a=1.339):
self.a = a
def _set_tuning_param(self, c, inplace=False):
"""
Set and change the tuning parameter of the Norm
Parameters
----------
c : float
The new value for the tuning parameter `a`.
inplace : bool, optional
If True, modify this instance in place and return it. If
False, return a new instance with the new tuning parameter.
The default is False.
Returns
-------
RobustNorm
The norm instance with the new tuning parameter.
Warnings
--------
This needs to wipe cached attributes that depend on the parameter.
"""
if inplace:
self.a = c
return self
else:
return self.__class__(a=c)
[docs]
def max_rho(self):
return 2 * self.a**2
def _subset(self, z):
"""
Andrew's wave is defined piecewise over the range of z
Parameters
----------
z : array_like
1d array
Returns
-------
ndarray
Boolean array indicating which elements of `z` fall in the
nonzero region of psi.
"""
z = np.asarray(z)
return np.less_equal(np.abs(z), self.a * np.pi)
[docs]
def rho(self, z):
r"""
The robust criterion function for Andrew's wave
.. math::
\rho(z) = \begin{cases}
a^2 (1 - \cos(z/a)) & \text{for } |z| \le a\pi \\
2a^2 & \text{for } |z| > a\pi
\end{cases}
Parameters
----------
z : array_like
1d array
Returns
-------
rho : ndarray
The value of the robust criterion function.
"""
a = self.a
z = np.asarray(z)
test = self._subset(z)
return (test * a**2 * (1 - np.cos(z / a)) +
(1 - test) * a**2 * 2)
[docs]
def psi(self, z):
r"""
The psi function for Andrew's wave
The analytic derivative of rho.
.. math::
\psi(z) = \begin{cases}
a \sin(z/a) & \text{for } |z| \le a\pi \\
0 & \text{for } |z| > a\pi
\end{cases}
Parameters
----------
z : array_like
1d array
Returns
-------
psi : ndarray
The value of the psi function.
"""
a = self.a
z = np.asarray(z)
test = self._subset(z)
return test * a * np.sin(z / a)
[docs]
def weights(self, z):
r"""
Andrew's wave weighting function for the IRLS algorithm
The psi function scaled by z.
.. math::
w(z) = \begin{cases}
\sin(z/a) / (z/a) & \text{for } |z| \le a\pi \\
0 & \text{for } |z| > a\pi
\end{cases}
Parameters
----------
z : array_like
1d array
Returns
-------
weights : ndarray
The value of the weighting function.
"""
a = self.a
z = np.asarray(z)
test = self._subset(z)
ratio = z / a
small = np.abs(ratio) < np.finfo(np.double).eps
if np.any(small):
weights = np.ones_like(ratio)
large = ~small
ratio = ratio[large]
weights[large] = test[large] * np.sin(ratio) / ratio
else:
weights = test * np.sin(ratio) / ratio
return weights
[docs]
def psi_deriv(self, z):
"""
The derivative of Andrew's wave psi function
Parameters
----------
z : array_like
1d array
Returns
-------
ndarray
The value of the derivative of the psi function.
Notes
-----
Used to estimate the robust covariance matrix.
"""
test = self._subset(z)
return test * np.cos(z / self.a)
[docs]
class TrimmedMean(RobustNorm):
"""
Trimmed mean function for M-estimation
Parameters
----------
c : float, optional
The tuning constant for the trimmed mean function. The default
value is 2.0.
See Also
--------
statsmodels.robust.norms.RobustNorm
"""
continuous = 0
redescending = "hard"
def __init__(self, c=2.):
self.c = c
def _set_tuning_param(self, c, inplace=False):
"""
Set and change the tuning parameter of the Norm
Parameters
----------
c : float
The new value for the tuning parameter `c`.
inplace : bool, optional
If True, modify this instance in place and return it. If
False, return a new instance with the new tuning parameter.
The default is False.
Returns
-------
RobustNorm
The norm instance with the new tuning parameter.
Warnings
--------
This needs to wipe cached attributes that depend on the parameter.
"""
if inplace:
self.c = c
return self
else:
return self.__class__(c=c)
[docs]
def max_rho(self):
return self.rho(self.c)
def _subset(self, z):
"""
Least trimmed mean is defined piecewise over the range of z
Parameters
----------
z : array_like
1d array
Returns
-------
ndarray
Boolean array indicating which elements of `z` fall in the
quadratic (rather than constant) region of rho.
"""
z = np.asarray(z)
return np.less_equal(np.abs(z), self.c)
[docs]
def rho(self, z):
r"""
The robust criterion function for least trimmed mean
.. math::
\rho(z) = \begin{cases}
\frac{1}{2} z^2 & \text{for } |z| \le c \\
\frac{1}{2} c^2 & \text{for } |z| > c
\end{cases}
Parameters
----------
z : array_like
1d array
Returns
-------
rho : ndarray
The value of the robust criterion function.
"""
z = np.asarray(z)
test = self._subset(z)
return test * z**2 * 0.5 + (1 - test) * self.c**2 * 0.5
[docs]
def psi(self, z):
r"""
The psi function for least trimmed mean
The analytic derivative of rho.
.. math::
\psi(z) = \begin{cases}
z & \text{for } |z| \le c \\
0 & \text{for } |z| > c
\end{cases}
Parameters
----------
z : array_like
1d array
Returns
-------
psi : ndarray
The value of the psi function.
"""
z = np.asarray(z)
test = self._subset(z)
return test * z
[docs]
def weights(self, z):
r"""
Least trimmed mean weighting function for the IRLS algorithm
The psi function scaled by z.
.. math::
w(z) = \begin{cases}
1 & \text{for } |z| \le c \\
0 & \text{for } |z| > c
\end{cases}
Parameters
----------
z : array_like
1d array
Returns
-------
weights : ndarray
The value of the weighting function.
"""
z = np.asarray(z)
test = self._subset(z)
return test
[docs]
def psi_deriv(self, z):
"""
The derivative of least trimmed mean psi function
Parameters
----------
z : array_like
1d array
Returns
-------
ndarray
The value of the derivative of the psi function.
Notes
-----
Used to estimate the robust covariance matrix.
"""
test = self._subset(z)
return test
[docs]
class Hampel(RobustNorm):
"""
Hampel function for M-estimation
Parameters
----------
a, b, c : float, optional
The tuning constants for Hampel's function. The default values are
a, b, c = 2, 4, 8.
See Also
--------
statsmodels.robust.norms.RobustNorm
"""
continuous = 1
redescending = "hard"
def __init__(self, a=2., b=4., c=8.):
self.a = a
self.b = b
self.c = c
def _set_tuning_param(self, c, inplace=False):
"""
Set and change the tuning parameter of the Norm
The tuning constants `a` and `b` are re-derived from `c` as
``a = c / 4`` and ``b = c / 2``.
Parameters
----------
c : float
The new value for the tuning parameter `c`.
inplace : bool, optional
If True, modify this instance in place and return it. If
False, return a new instance with the new tuning parameters.
The default is False.
Returns
-------
RobustNorm
The norm instance with the new tuning parameters.
Warnings
--------
This needs to wipe cached attributes that depend on the parameter.
"""
a = c / 4
b = c / 2
if inplace:
self.c = c
self.a = a
self.b = b
return self
else:
return self.__class__(a=a, b=b, c=c)
[docs]
def max_rho(self):
return self.rho(self.c)
def _subset(self, z):
"""
Hampel's function is defined piecewise over the range of z
Parameters
----------
z : array_like
1d array
Returns
-------
tuple of ndarray
Three boolean arrays, ``t1``, ``t2`` and ``t3``, indicating
which elements of `z` fall in each of the piecewise regions of
rho.
"""
z = np.abs(np.asarray(z))
t1 = np.less_equal(z, self.a)
t2 = np.less_equal(z, self.b) * np.greater(z, self.a)
t3 = np.less_equal(z, self.c) * np.greater(z, self.b)
return t1, t2, t3
[docs]
def rho(self, z):
r"""
The robust criterion function for Hampel's estimator
.. math::
\rho(z) = \begin{cases}
\frac{1}{2} z^2 & \text{for } |z| \le a \\
a|z| - \frac{1}{2} a^2 & \text{for } a < |z| \le b \\
\frac{a(c - |z|)^2}{2(c - b)} & \text{for } b < |z| \le c \\
\frac{a(b + c - a)}{2} & \text{for } |z| > c
\end{cases}
Parameters
----------
z : array_like
1d array
Returns
-------
rho : ndarray
The value of the robust criterion function.
"""
a, b, c = self.a, self.b, self.c
z_isscalar = np.isscalar(z)
z = np.atleast_1d(z)
t1, t2, t3 = self._subset(z)
t34 = ~(t1 | t2)
dt = np.promote_types(z.dtype, "float")
v = np.zeros(z.shape, dtype=dt)
z = np.abs(z)
v[t1] = z[t1]**2 * 0.5
# v[t2] = (a * (z[t2] - a) + a**2 * 0.5)
v[t2] = (a * z[t2] - a**2 * 0.5)
v[t3] = a * (c - z[t3])**2 / (c - b) * (-0.5)
v[t34] += a * (b + c - a) * 0.5
if z_isscalar:
v = v[0]
return v
[docs]
def psi(self, z):
r"""
The psi function for Hampel's estimator
The analytic derivative of rho.
.. math::
\psi(z) = \begin{cases}
z & \text{for } |z| \le a \\
a \cdot \text{sgn}(z) & \text{for } a < |z| \le b \\
a \cdot \text{sgn}(z) \frac{c - |z|}{c - b} & \text{for } b < |z| \le c \\
0 & \text{for } |z| > c
\end{cases}
Parameters
----------
z : array_like
1d array
Returns
-------
psi : ndarray
The value of the psi function.
"""
a, b, c = self.a, self.b, self.c
z_isscalar = np.isscalar(z)
z = np.atleast_1d(z)
t1, t2, t3 = self._subset(z)
dt = np.promote_types(z.dtype, "float")
v = np.zeros(z.shape, dtype=dt)
s = np.sign(z)
za = np.abs(z)
v[t1] = z[t1]
v[t2] = a * s[t2]
v[t3] = a * s[t3] * (c - za[t3]) / (c - b)
if z_isscalar:
v = v[0]
return v
[docs]
def weights(self, z):
r"""
Hampel weighting function for the IRLS algorithm
The psi function scaled by z.
.. math::
w(z) = \begin{cases}
1 & \text{for } |z| \le a \\
a/|z| & \text{for } a < |z| \le b \\
\frac{a(c - |z|)}{|z|(c - b)} & \text{for } b < |z| \le c \\
0 & \text{for } |z| > c
\end{cases}
Parameters
----------
z : array_like
1d array
Returns
-------
weights : ndarray
The value of the weighting function.
"""
a, b, c = self.a, self.b, self.c
z_isscalar = np.isscalar(z)
z = np.atleast_1d(z)
t1, t2, t3 = self._subset(z)
dt = np.promote_types(z.dtype, "float")
v = np.zeros(z.shape, dtype=dt)
v[t1] = 1.0
abs_z = np.abs(z)
v[t2] = a / abs_z[t2]
abs_zt3 = abs_z[t3]
v[t3] = a * (c - abs_zt3) / (abs_zt3 * (c - b))
if z_isscalar:
v = v[0]
return v
[docs]
def psi_deriv(self, z):
"""
Derivative of psi function, second derivative of rho function
Parameters
----------
z : array_like
1d array
Returns
-------
ndarray
The value of the derivative of the psi function.
Notes
-----
Used to estimate the robust covariance matrix.
"""
a, b, c = self.a, self.b, self.c
z_isscalar = np.isscalar(z)
z = np.atleast_1d(z)
t1, _, t3 = self._subset(z)
dt = np.promote_types(z.dtype, "float")
d = np.zeros(z.shape, dtype=dt)
d[t1] = 1.0
zt3 = z[t3]
d[t3] = -(a * np.sign(zt3) * zt3) / (np.abs(zt3) * (c - b))
if z_isscalar:
d = d[0]
return d
[docs]
class TukeyBiweight(RobustNorm):
"""
Tukey's biweight function for M-estimation
Parameters
----------
c : float, optional
The tuning constant for Tukey's Biweight. The default value is
c = 4.685.
Notes
-----
Tukey's biweight is sometimes called bisquare.
"""
continuous = 2
redescending = "hard"
def __init__(self, c=4.685):
self.c = c
def __repr__(self):
return f"{self.__class__.__name__}(c={self.c})"
[docs]
@classmethod
def get_tuning(cls, bp=None, eff=None):
"""
Tuning parameter for given breakdown point or efficiency
This currently only returns values from a table.
Parameters
----------
bp : float in [0.05, 0.5] or None
Required breakdown point
Either bp or eff has to be specified, but not both.
eff : float or None
Required asymptotic efficiency.
Either bp or eff has to be specified, but not both.
Returns
-------
float
The tuning parameter.
"""
if ((bp is None and eff is None) or
(bp is not None and eff is not None)):
raise ValueError("exactly one of bp and eff needs to be provided")
if bp is not None:
return _tables.tukeybiweight_bp[bp]
elif eff is not None:
return _tables.tukeybiweight_eff[eff]
def _set_tuning_param(self, c, inplace=False):
"""
Set and change the tuning parameter of the Norm
Parameters
----------
c : float
The new value for the tuning parameter `c`.
inplace : bool, optional
If True, modify this instance in place and return it. If
False, return a new instance with the new tuning parameter.
The default is False.
Returns
-------
RobustNorm
The norm instance with the new tuning parameter.
Warnings
--------
This needs to wipe cached attributes that depend on the parameter.
"""
# todo : change default to inplace=False, when tools are fixed
if inplace:
self.c = c
return self
else:
return self.__class__(c=c)
[docs]
def max_rho(self):
return self.rho(self.c)
def _subset(self, z):
"""
Tukey's biweight is defined piecewise over the range of z
Parameters
----------
z : array_like
1d array
Returns
-------
ndarray
Boolean array indicating which elements of `z` fall in the
nonzero region of psi.
"""
z = np.abs(np.asarray(z))
return np.less_equal(z, self.c)
[docs]
def rho(self, z):
r"""
The robust criterion function for Tukey's biweight estimator
Parameters
----------
z : array_like
1d array
Returns
-------
rho : ndarray
rho(z) = (c**2 / 6) * (1 - (1 - (z/c)**2)**3) for \|z\| <= c
rho(z) = c**2 / 6 for \|z\| > c
"""
subset = self._subset(z)
factor = self.c**2 / 6.
return -(1 - (z / self.c)**2)**3 * subset * factor + factor
[docs]
def psi(self, z):
r"""
The psi function for Tukey's biweight estimator
The analytic derivative of rho
Parameters
----------
z : array_like
1d array
Returns
-------
psi : ndarray
psi(z) = z*(1 - (z/c)**2)**2 for \|z\| <= c
psi(z) = 0 for \|z\| > c
"""
z = np.asarray(z)
subset = self._subset(z)
return z * (1 - (z / self.c)**2)**2 * subset
[docs]
def weights(self, z):
r"""
Tukey's biweight weighting function for the IRLS algorithm
The psi function scaled by z
Parameters
----------
z : array_like
1d array
Returns
-------
weights : ndarray
weights(z) = (1 - (z/c)**2)**2 for \|z\| <= c
weights(z) = 0 for \|z\| > c
"""
z = np.asarray(z)
subset = self._subset(z)
return (1 - (z / self.c)**2)**2 * subset
[docs]
def psi_deriv(self, z):
"""
The derivative of Tukey's biweight psi function
Parameters
----------
z : array_like
1d array
Returns
-------
ndarray
The value of the derivative of the psi function.
Notes
-----
Used to estimate the robust covariance matrix.
"""
subset = self._subset(z)
return subset * ((1 - (z/self.c)**2)**2
- (4*z**2/self.c**2) * (1-(z/self.c)**2))
class TukeyQuartic(RobustNorm):
"""
Variant of Tukey's biweight function with power 4 for M-estimation
Parameters
----------
c : float, optional
The tuning constant for Tukey's Biweight. The default value is
c = 3.61752.
k : int, optional
The power used in the inner term in place of the power 2 used by
Tukey's biweight. The default value is 4.
Notes
-----
This is a variation of Tukey's biweight (bisquare) function where
the weight function has power 4 instead of power 2 in the inner term.
"""
continuous = 2
redescending = "hard"
def __init__(self, c=3.61752, k=4):
# TODO: c needs to be changed if k != 4
# also, I think implementation assumes k is even integer
self.c = c
self.k = k
def _set_tuning_param(self, c, inplace=False):
"""
Set and change the tuning parameter of the Norm
Parameters
----------
c : float
The new value for the tuning parameter `c`.
inplace : bool, optional
If True, modify this instance in place and return it. If
False, return a new instance with the new tuning parameter.
The default is False.
Returns
-------
RobustNorm
The norm instance with the new tuning parameter.
Warnings
--------
This needs to wipe cached attributes that depend on the parameter.
"""
if inplace:
self.c = c
return self
else:
return self.__class__(c=c, k=self.k)
def max_rho(self):
return self.rho(self.c)
def _subset(self, z):
"""
TukeyQuartic is defined piecewise over the range of z
Parameters
----------
z : array_like
1d array
Returns
-------
ndarray
Boolean array indicating which elements of `z` fall in the
nonzero region of psi.
"""
z = np.abs(np.asarray(z))
return np.less_equal(z, self.c)
def rho(self, z):
r"""
The robust criterion function for TukeyQuartic norm
.. math::
\rho(z) = \begin{cases}
\frac{1}{2} z^2 \left(1 - \frac{4}{k + 2} x^k + \frac{1}{k + 1} x^{2k}\right) & \text{for } |z| \le c \\
\rho(c) & \text{for } | z | > c
\end{cases}
where :math:`x = z / c`.
Parameters
----------
z : array_like
1d array
Returns
-------
rho : ndarray
The value of the robust criterion function.
"""
c = self.c
k = self.k
subset = self._subset(z)
x = z / c
rhoc = 1 / 2 * c**2 * (1 - 4 / (k + 2) + 1 / (k + 1))
# integral x (1 - x^k)^2 dx =
# 1/2 x^2 (x^(2 k)/(k + 1) - (4 x^k)/(k + 2) + 1) + constant
# integral x (1 - (x/c)^k)^2 dx =
# 1/2 x^2 (-(4 (x/c)^k)/(k + 2) + (x/c)^(2 k)/(k + 1) + 1) +
# constant
rh = (
subset
* 1
/ 2
* z**2
* (1 - 4 / (k + 2) * x**k + 1 / (k + 1) * x ** (2 * k))
+ (1 - subset) * rhoc
)
return rh
def psi(self, z):
r"""
The psi function of TukeyQuartic norm
The analytic derivative of rho.
Parameters
----------
z : array_like
1d array
Returns
-------
psi : ndarray
psi(z) = z*(1 - (z/c)**k)**2 for \|z\| <= c
psi(z) = 0 for \|z\| > c
"""
k = self.k
z = np.asarray(z)
subset = self._subset(z)
return z * (1 - (z / self.c)**k)**2 * subset
def weights(self, z):
r"""
TukeyQuartic weighting function for the IRLS algorithm
The psi function scaled by z.
Parameters
----------
z : array_like
1d array
Returns
-------
weights : ndarray
weights(z) = (1 - (z/c)**k)**2 for \|z\| <= c
weights(z) = 0 for \|z\| > c
"""
k = self.k
z = np.asarray(z)
subset = self._subset(z)
return (1 - (z / self.c)**k)**2 * subset
def psi_deriv(self, z):
"""
The derivative of the TukeyQuartic psi function
Parameters
----------
z : array_like
1d array
Returns
-------
ndarray
The value of the derivative of the psi function.
Notes
-----
Used to estimate the robust covariance matrix.
"""
c = self.c
k = self.k
subset = self._subset(z)
x = z / c
# d/dx(x (1 - (x/c)^k)^2) = -(1 - (x/c)^k) (2 k (x/c)^k + (x/c)^k - 1)
return subset * (1 - x**k) * (1 - (2 * k + 1) * x**k)
class StudentT(RobustNorm):
"""
Robust norm based on t distribution
Rho is a rescaled version of the t-loglikelihood function after dropping
constant terms.
The norms are rescaled so that the largest weights are 1 and
the second derivative of the rho function at zero is equal to 1.
The maximum likelihood estimator based on the loglikelihood
function of the t-distribution is available in
``statsmodels.miscmodels``, which can be used to also
estimate scale and degrees of freedom by MLE.
Parameters
----------
c : float, optional
The tuning constant for the StudentT norm. The default value is
2.3849.
df : float, optional
The degrees of freedom of the t distribution used to derive the
norm. The default value is 4.
"""
continuous = 2
redescending = "soft"
def __init__(self, c=2.3849, df=4):
self.c = c
self.df = df
def _set_tuning_param(self, c, inplace=False):
"""
Set and change the tuning parameter of the Norm
Parameters
----------
c : float
The new value for the tuning parameter `c`.
inplace : bool, optional
If True, modify this instance in place and return it. If
False, return a new instance with the new tuning parameter.
The default is False.
Returns
-------
RobustNorm
The norm instance with the new tuning parameter.
Warnings
--------
This needs to wipe cached attributes that depend on the parameter.
"""
if inplace:
self.c = c
return self
else:
return self.__class__(c=c, df=self.df)
def max_rho(self):
return np.inf
def rho(self, z):
"""
The rho function of the StudentT norm
Parameters
----------
z : array_like
1d array
Returns
-------
rho : ndarray
rho(z) = (c**2 * df / 2.) * log(df + (z / c)**2) - const
The ``const`` shifts the rho function so that rho(0) = 0.
"""
c = self.c
df = self.df
z = np.asarray(z)
const = (c**2 * df / 2.) * np.log(df) if df != 0 else 0
return (c**2 * df / 2.) * np.log(df + (z / c)**2) - const
def psi(self, z):
"""
The psi function of the StudentT norm
The analytic derivative of rho.
Parameters
----------
z : array_like
1d array
Returns
-------
psi : ndarray
psi(z) = z * df / (df + (z / c)**2)
"""
c = self.c
df = self.df
z = np.asarray(z)
return z * df / (df + (z / c)**2)
def weights(self, z):
"""
The weighting function for the IRLS algorithm of the StudentT norm
The psi function scaled by the input z
Parameters
----------
z : array_like
1d array
Returns
-------
weights : ndarray
weights(z) = df / (df + (z / c)**2)
"""
c = self.c
df = self.df
z = np.asarray(z)
return df / (df + (z / c)**2)
def psi_deriv(self, z):
"""
The derivative of the psi function of the StudentT norm
Parameters
----------
z : array_like
1d array
Returns
-------
psi_deriv : ndarray
psi_deriv(z) = -2 * df * x**2 / (df + x**2)**2 + df / (df + x**2)
where ``x = z / c``.
Notes
-----
Used to estimate the robust covariance matrix.
"""
c = self.c
df = self.df
x = np.asarray(z) / c
return - 2 * df * x**2 / (df + x**2)**2 + df / (df + x**2)
[docs]
class MQuantileNorm(RobustNorm):
"""
M-quantiles objective function based on a base norm
This norm has the same asymmetric structure as the objective function
in QuantileRegression but replaces the L1 absolute value by a chosen
base norm.
rho_q(u) = abs(q - I(q < 0)) * rho_base(u)
or, equivalently,
rho_q(u) = q * rho_base(u) if u >= 0
rho_q(u) = (1 - q) * rho_base(u) if u < 0
Parameters
----------
q : float
M-quantile, must be between 0 and 1
base_norm : RobustNorm instance
Basic norm that is transformed into an asymmetric M-quantile norm
Notes
-----
This is mainly for base norms that are not redescending, like HuberT or
LeastSquares. (See Jones for the relationship of M-quantiles to quantiles
in the case of non-redescending Norms.)
Expectiles are M-quantiles with the LeastSquares as base norm.
References
----------
.. [*] Bianchi, Annamaria, and Nicola Salvati. 2015. “Asymptotic Properties
and Variance Estimators of the M-Quantile Regression Coefficients
Estimators.” Communications in Statistics - Theory and Methods 44 (11):
2416-29. doi:10.1080/03610926.2013.791375.
.. [*] Breckling, Jens, and Ray Chambers. 1988. “M-Quantiles.”
Biometrika 75 (4): 761-71. doi:10.2307/2336317.
.. [*] Jones, M. C. 1994. “Expectiles and M-Quantiles Are Quantiles.”
Statistics & Probability Letters 20 (2): 149-53.
doi:10.1016/0167-7152(94)90031-0.
.. [*] Newey, Whitney K., and James L. Powell. 1987. “Asymmetric Least
Squares Estimation and Testing.” Econometrica 55 (4): 819-47.
doi:10.2307/1911031.
"""
continuous = 1
def __init__(self, q, base_norm):
self.q = q
self.base_norm = base_norm
def _get_q(self, z):
nobs = len(z)
mask_neg = (z < 0) # if self.q < 0.5 else (z <= 0) # maybe symmetric
qq = np.empty(nobs)
qq[mask_neg] = 1 - self.q
qq[~mask_neg] = self.q
return qq
[docs]
def rho(self, z):
"""
The robust criterion function for MQuantileNorm
Parameters
----------
z : array_like
1d array
Returns
-------
rho : ndarray
The value of the robust criterion function.
"""
qq = self._get_q(z)
return qq * self.base_norm.rho(z)
[docs]
def psi(self, z):
"""
The psi function for MQuantileNorm estimator
The analytic derivative of rho
Parameters
----------
z : array_like
1d array
Returns
-------
psi : ndarray
The value of the psi function.
"""
qq = self._get_q(z)
return qq * self.base_norm.psi(z)
[docs]
def weights(self, z):
"""
MQuantileNorm weighting function for the IRLS algorithm
The psi function scaled by z, psi(z) / z
Parameters
----------
z : array_like
1d array
Returns
-------
weights : ndarray
The value of the weighting function.
"""
qq = self._get_q(z)
return qq * self.base_norm.weights(z)
[docs]
def psi_deriv(self, z):
"""
The derivative of MQuantileNorm function
Parameters
----------
z : array_like
1d array
Returns
-------
psi_deriv : ndarray
The value of the derivative of the psi function.
Notes
-----
Used to estimate the robust covariance matrix.
"""
qq = self._get_q(z)
return qq * self.base_norm.psi_deriv(z)
def __call__(self, z):
"""
Return the value of estimator rho applied to an input
Parameters
----------
z : array_like
1d array
Returns
-------
ndarray
The value of the robust criterion function.
"""
return self.rho(z)
[docs]
def estimate_location(a, scale, norm=None, axis=0, initial=None,
maxiter=30, tol=1.0e-06):
"""
Estimate a robust location parameter using an M-estimator.
This function iteratively computes the location parameterthat satisfies:
sum(psi((a-mu)/scale))=0
where 'psi' is the influence function defined by the selected robust norm.
Parameters
----------
a : ndarray
Array over which the location parameter is to be estimated
scale : ndarray
Scale parameter to be used in M-estimator
norm : RobustNorm, optional
Robust norm used in the M-estimator. The default is HuberT().
axis : int, optional
Axis along which to estimate the location parameter. The default is 0.
initial : ndarray, optional
Initial condition for the location parameter. Default is None, which
uses the median of a.
maxiter : int, optional
Maximum number of iterations. The default is 30.
tol : float, optional
Toleration for convergence. The default is 1e-06.
Returns
-------
mu : ndarray
Estimate of location
"""
if norm is None:
norm = HuberT()
if initial is None:
mu = np.median(a, axis)
else:
mu = initial
for _ in range(maxiter):
W = norm.weights((a-mu)/scale)
denom = np.sum(W, axis)
if np.any(denom <= 0):
return mu
nmu = np.sum(W * a, axis) / denom
if np.all(np.less(np.abs(mu - nmu), scale * tol)):
return nmu
else:
mu = nmu
raise ValueError(f"location estimator failed to converge in {maxiter:d} iterations")
