import numpy as np
from robpy.univariate.base import RobustScale
from scipy.stats import median_abs_deviation, norm
[docs]
class Tau(RobustScale):
def __init__(
self,
c1: float = 4.5,
c2: float = 3.0,
consistency_correction: bool = True,
):
"""
Implementation of tau estimator of scale of Yohai, V.J. & Zamar, R.H. (1988).
Args:
c1 (float, optional):
Constant for the weight function, defaults to 4.5.
c2 (float, optional):
Constant for the rho function, defaults to 3.0.
consistency_correction (bool, optional):
boolean indicating if consistency for normality should be applied.
Defaults to True.
References:
- Yohai, V.J. & Zamar, R.H. (1988). High breakdown estimates of regression by means of
the minimization of an efficient scale. Journal of the American Statistical
Association, 83(402), 406-413.
"""
super().__init__()
self.c1 = c1
self.c2 = c2
self.consistency_correction = consistency_correction
def _calculate(self, X: np.ndarray):
"""
Args:
X (np.ndarray):
Univariate data.
"""
n = len(X)
sigma0 = median_abs_deviation(X)
weights = self._weight_function((X - np.median(X)) / sigma0)
self.location_ = np.sum(X * weights) / np.sum(weights)
self.scale_ = sigma0 * np.sqrt(
1 / n * np.sum(self._rho_function((X - self.location_) / sigma0))
)
if self.consistency_correction:
# expectation of rho(X/qnorm(3/4)) for X standard normal
b = self.c2 * norm.ppf(3 / 4)
corr = 2 * ((1 - b**2) * norm.cdf(b) - b * norm.pdf(b) + b**2) - 1
self.scale_ = self.scale_ / np.sqrt(corr)
def _weight_function(self, X):
return np.where(
np.abs(X) <= self.c1,
(1 - (X / self.c1) ** 2) ** 2,
0.0,
)
def _rho_function(self, X):
return np.where(X**2 <= self.c2**2, X**2, self.c2**2)