Repository: statsmodels/statsmodels

Issue:
BUG/TST robust norms bugs and unit tests
Robust `norms` do not have their own unit tests. They are only indirectly tested through RLM results.

old PR #1341 included some fixes for Hampel bugs and shifted tukey biweight.
see also #8697 and #8696 for Hampel bugs

Methods of norms used in RLM

- fit only uses `weights`
- cov_params, bcov_scaled uses `psi_deriv` and `psi`
- deviance uses `rho` (`__call__`)

unit tests possible against R robustbase
(in my old PR, I did not have reference numbers, I used mostly graphical inspection and checking formulas)

e.g.
```
> Mpsi(c(-9, -6, -2, -1, 0, 1, 2, 6, 9), cc = c(1.5, 3.5, 8), psi = "hampel")
[1]  0.0000000 -0.6666667 -1.5000000 -1.0000000  0.0000000  1.0000000  1.5000000
[8]  0.6666667  0.0000000
> Mpsi(c(-9, -6, -2, -1, 0, 1, 2, 6, 9), cc = c(1.5, 3.5, 8), psi = "hampel", deriv=-1)  # unscaled rho
[1] 7.500000 6.833333 1.875000 0.500000 0.000000 0.500000 1.875000 6.833333
[9] 7.500000
```



Discussion:
josef-pkt: from https://github.com/statsmodels/statsmodels/issues/7603#issuecomment-885884017

at least dtype can be implemented with np.promote_types
```
np.promote_types(np.array([990000000, 2100000000], np.int32).dtype, np.float64)
dtype('float64')

np.promote_types(np.array([990000000, 2100000000], np.complex128).dtype, np.float64)
dtype('complex128')
```

Patch:
diff --git a/statsmodels/robust/norms.py b/statsmodels/robust/norms.py
index ba425e579..3715bd82e 100644
--- a/statsmodels/robust/norms.py
+++ b/statsmodels/robust/norms.py
@@ -3,6 +3,16 @@ import numpy as np
 # TODO: add plots to weighting functions for online docs.
 
 
+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.
+    """
+    sign = (x.real >= 0) * 2 - 1
+    return sign * x
+
+
 class RobustNorm:
     """
     The parent class for the norms used for robust regression.
@@ -245,11 +255,17 @@ class HuberT(RobustNorm):
 
             weights(z) = t/\|z\|      for \|z\| > t
         """
-        z = np.asarray(z)
+        z_isscalar = np.isscalar(z)
+        z = np.atleast_1d(z)
+
         test = self._subset(z)
         absz = np.abs(z)
         absz[test] = 1.0
-        return test + (1 - test) * self.t / absz
+        v = test + (1 - test) * self.t / absz
+
+        if z_isscalar:
+            v = v[0]
+        return v
 
     def psi_deriv(self, z):
         """
@@ -259,7 +275,7 @@ class HuberT(RobustNorm):
         -----
         Used to estimate the robust covariance matrix.
         """
-        return np.less_equal(np.abs(z), self.t)
+        return np.less_equal(np.abs(z), self.t).astype(float)
 
 
 # TODO: untested, but looks right.  RamsayE not available in R or SAS?
@@ -390,16 +406,16 @@ class AndrewWave(RobustNorm):
         Returns
         -------
         rho : ndarray
-            rho(z) = a*(1-cos(z/a))     for \|z\| <= a*pi
+            rho(z) = a**2 *(1-cos(z/a))     for \|z\| <= a*pi
 
-            rho(z) = 2*a                for \|z\| > a*pi
+            rho(z) = 2*a                    for \|z\| > a*pi
         """
 
         a = self.a
         z = np.asarray(z)
         test = self._subset(z)
-        return (test * a * (1 - np.cos(z / a)) +
-                (1 - test) * 2 * a)
+        return (test * a**2 * (1 - np.cos(z / a)) +
+                (1 - test) * a**2 * 2)
 
     def psi(self, z):
         r"""
@@ -415,7 +431,7 @@ class AndrewWave(RobustNorm):
         Returns
         -------
         psi : ndarray
-            psi(z) = sin(z/a)       for \|z\| <= a*pi
+            psi(z) = a * sin(z/a)   for \|z\| <= a*pi
 
             psi(z) = 0              for \|z\| > a*pi
         """
@@ -423,7 +439,7 @@ class AndrewWave(RobustNorm):
         a = self.a
         z = np.asarray(z)
         test = self._subset(z)
-        return test * np.sin(z / a)
+        return test * a * np.sin(z / a)
 
     def weights(self, z):
         r"""
@@ -439,9 +455,9 @@ class AndrewWave(RobustNorm):
         Returns
         -------
         weights : ndarray
-            weights(z) = sin(z/a)/(z/a)     for \|z\| <= a*pi
+            weights(z) = sin(z/a) / (z/a)     for \|z\| <= a*pi
 
-            weights(z) = 0                  for \|z\| > a*pi
+            weights(z) = 0                    for \|z\| > a*pi
         """
         a = self.a
         z = np.asarray(z)
@@ -467,7 +483,7 @@ class AndrewWave(RobustNorm):
         """
 
         test = self._subset(z)
-        return test*np.cos(z / self.a)/self.a
+        return test * np.cos(z / self.a)
 
 
 # TODO: this is untested
@@ -511,12 +527,12 @@ class TrimmedMean(RobustNorm):
         rho : ndarray
             rho(z) = (1/2.)*z**2    for \|z\| <= c
 
-            rho(z) = 0              for \|z\| > c
+            rho(z) = (1/2.)*c**2              for \|z\| > c
         """
 
         z = np.asarray(z)
         test = self._subset(z)
-        return test * z**2 * 0.5
+        return test * z**2 * 0.5 + (1 - test) * self.c**2 * 0.5
 
     def psi(self, z):
         r"""
@@ -619,24 +635,33 @@ class Hampel(RobustNorm):
         Returns
         -------
         rho : ndarray
-            rho(z) = (1/2.)*z**2                    for \|z\| <= a
+            rho(z) = z**2 / 2                     for \|z\| <= a
 
-            rho(z) = a*\|z\| - 1/2.*a**2              for a < \|z\| <= b
+            rho(z) = a*\|z\| - 1/2.*a**2               for a < \|z\| <= b
 
-            rho(z) = a*(c*\|z\|-(1/2.)*z**2)/(c-b)    for b < \|z\| <= c
+            rho(z) = a*(c - \|z\|)**2 / (c - b) / 2    for b < \|z\| <= c
 
-            rho(z) = a*(b + c - a)                  for \|z\| > c
+            rho(z) = a*(b + c - a) / 2                 for \|z\| > c
         """
+        a, b, c = self.a, self.b, self.c
+
+        z_isscalar = np.isscalar(z)
+        z = np.atleast_1d(z)
 
-        z = np.abs(z)
-        a = self.a
-        b = self.b
-        c = self.c
         t1, t2, t3 = self._subset(z)
-        v = (t1 * z**2 * 0.5 +
-             t2 * (a * z - a**2 * 0.5) +
-             t3 * (a * (c * z - z**2 * 0.5) / (c - b) - 7 * a**2 / 6.) +
-             (1 - t1 + t2 + t3) * a * (b + c - a))
+        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
 
     def psi(self, z):
@@ -661,16 +686,23 @@ class Hampel(RobustNorm):
 
             psi(z) = 0                            for \|z\| > c
         """
-        z = np.asarray(z)
-        a = self.a
-        b = self.b
-        c = self.c
+        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)
-        z = np.abs(z)
-        v = s * (t1 * z +
-                 t2 * a*s +
-                 t3 * a*s * (c - z) / (c - b))
+        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
 
     def weights(self, z):
@@ -687,37 +719,51 @@ class Hampel(RobustNorm):
         Returns
         -------
         weights : ndarray
-            weights(z) = 1                            for \|z\| <= a
+            weights(z) = 1                                for \|z\| <= a
 
-            weights(z) = a/\|z\|                        for a < \|z\| <= b
+            weights(z) = a/\|z\|                          for a < \|z\| <= b
 
             weights(z) = a*(c - \|z\|)/(\|z\|*(c-b))      for b < \|z\| <= c
 
-            weights(z) = 0                            for \|z\| > c
+            weights(z) = 0                                for \|z\| > c
         """
-        z = np.asarray(z)
-        a = self.a
-        b = self.b
-        c = self.c
+        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)
 
-        v = np.zeros_like(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))
-        v[np.where(np.isnan(v))] = 1.  # TODO: for some reason 0 returns a nan?
+
+        if z_isscalar:
+            v = v[0]
         return v
 
     def psi_deriv(self, z):
-        t1, _, t3 = self._subset(z)
+        """Derivative of psi function, second derivative of rho function.
+        """
         a, b, c = self.a, self.b, self.c
-        # default is t1
-        d = np.zeros_like(z)
+
+        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))
+        d[t3] = -(a * np.sign(zt3) * zt3) / (np.abs(zt3) * (c - b))
+
+        if z_isscalar:
+            d = d[0]
         return d
 
 
@@ -764,7 +810,8 @@ class TukeyBiweight(RobustNorm):
             rho(z) = 0                              for \|z\| > R
         """
         subset = self._subset(z)
-        return -(1 - (z / self.c)**2)**3 * subset * self.c**2 / 6.
+        factor = self.c**2 / 6.
+        return -(1 - (z / self.c)**2)**3 * subset * factor + factor
 
     def psi(self, z):
         r"""


Test patch:
diff --git a/statsmodels/robust/tests/results/results_norms.py b/statsmodels/robust/tests/results/results_norms.py
new file mode 100644
index 000000000..2be610239
--- /dev/null
+++ b/statsmodels/robust/tests/results/results_norms.py
@@ -0,0 +1,33 @@
+
+
+from statsmodels.tools.testing import Holder
+
+res_hampel = Holder(
+    rho=[7.5, 6.833333333333333, 1.875, 0.5, 0.0, 0.5, 1.875,
+         6.833333333333333, 7.5],
+    psi=[0.0, -0.6666666666666666, -1.5, -1.0, 0.0, 1.0, 1.5,
+         0.6666666666666666, 0.0],
+    psi_deriv=[0.0, -0.3333333333333333, 0.0, 1.0, 1.0, 1.0, 0.0,
+               -0.3333333333333333, 0.0],
+    weights=[0.0, 0.1111111111111111, 0.75, 1.0, 1.0, 1.0, 0.75,
+             0.1111111111111111, 0.0],
+    )
+
+res_biweight = Holder(
+    rho=[2.6666666666666665, 2.6666666666666665, 1.5416666666666665,
+         0.46940104166666663, 0.0, 0.46940104166666663, 1.5416666666666665,
+         2.6666666666666665, 2.6666666666666665],
+    psi=[0.0, 0.0, -1.125, -0.87890625, 0.0, 0.87890625, 1.125, 0.0, 0.0],
+    psi_deriv=[0.0, 0.0, -0.1875, 0.64453125, 1.0, 0.64453125, -0.1875, 0.0,
+               0.0],
+    weights=[0.0, 0.0, 0.5625, 0.87890625, 1.0, 0.87890625, 0.5625, 0.0, 0.0],
+    )
+
+res_huber = Holder(
+    rho=[11.200487500000001, 7.165487499999999, 1.7854875000000001, 0.5, 0.0,
+         0.5, 1.7854875000000001, 7.165487499999999, 11.200487500000001],
+    psi=[-1.345, -1.345, -1.345, -1.0, 0.0, 1.0, 1.345, 1.345, 1.345],
+    psi_deriv=[0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0],
+    weights=[0.14944444444444444, 0.22416666666666665, 0.6725, 1.0, 1.0, 1.0,
+             0.6725, 0.22416666666666665, 0.14944444444444444],
+    )
diff --git a/statsmodels/robust/tests/test_norms.py b/statsmodels/robust/tests/test_norms.py
new file mode 100644
index 000000000..9c0ba5f5a
--- /dev/null
+++ b/statsmodels/robust/tests/test_norms.py
@@ -0,0 +1,96 @@
+
+import pytest
+import numpy as np
+from numpy.testing import assert_allclose
+
+from statsmodels.robust import norms
+from statsmodels.tools.numdiff import (
+    _approx_fprime_scalar,
+    # _approx_fprime_cs_scalar,  # not yet
+    )
+from .results import results_norms as res_r
+
+cases = [
+    (norms.Hampel, (1.5, 3.5, 8.), res_r.res_hampel),
+    (norms.TukeyBiweight, (4,), res_r.res_biweight),
+    (norms.HuberT, (1.345,), res_r.res_huber),
+    ]
+
+norms_other = [
+    (norms.LeastSquares, ()),
+    (norms.TrimmedMean, (1.9,)),  # avoid arg at integer used in example
+    (norms.AndrewWave, ()),
+    (norms.RamsayE, ()),
+    # norms.MQuantileNorm,  # requires keywords in init
+    ]
+
+dtypes = ["int", np.float64, np.complex128]
+
+
+@pytest.mark.parametrize("dtype", dtypes)
+@pytest.mark.parametrize("case", cases)
+def test_norm(case, dtype):
+    ncls, args, res = case
+    if ncls in [norms.HuberT] and dtype == np.complex128:
+        # skip for now
+        return
+
+    norm = ncls(*args)
+    x = np.array([-9, -6, -2, -1, 0, 1, 2, 6, 9], dtype=dtype)
+
+    weights = norm.weights(x)
+    rho = norm.rho(x)
+    psi = norm.psi(x)
+    psi_deriv = norm.psi_deriv(x)
+    assert_allclose(weights, res.weights, rtol=1e-12, atol=1e-20)
+    assert_allclose(rho, res.rho, rtol=1e-12, atol=1e-20)
+    assert_allclose(psi, res.psi, rtol=1e-12, atol=1e-20)
+    assert_allclose(psi_deriv, res.psi_deriv, rtol=1e-12, atol=1e-20)
+
+    dtype2 = np.promote_types(dtype, "float")
+    assert weights.dtype == dtype2
+    assert rho.dtype == dtype2
+    assert psi.dtype == dtype2
+    assert psi_deriv.dtype == dtype2
+
+    psid = _approx_fprime_scalar(x, norm.rho)
+    assert_allclose(psid, res.psi, rtol=1e-6, atol=1e-8)
+    psidd = _approx_fprime_scalar(x, norm.psi)
+    assert_allclose(psidd, res.psi_deriv, rtol=1e-6, atol=1e-8)
+
+    # complex step derivatives are not yet supported if method uses np.abs
+    # psid = _approx_fprime_cs_scalar(x, norm.rho)
+    # assert_allclose(psid, res.psi, rtol=1e-12, atol=1e-20)
+    # psidd = _approx_fprime_cs_scalar(x, norm.psi)
+    # assert_allclose(psidd, res.psi_deriv, rtol=1e-12, atol=1e-20)
+
+    # check scalar value
+    methods = ["weights", "rho", "psi", "psi_deriv"]
+    for meth in methods:
+        resm = [getattr(norm, meth)(xi) for xi in x]
+        assert_allclose(resm, getattr(res, meth))
+
+
+@pytest.mark.parametrize("case", norms_other)
+def test_norms_consistent(case):
+    # test that norm methods are consistent with each other
+    ncls, args = case
+    norm = ncls(*args)
+    x = np.array([-9, -6, -2, -1, 0, 1, 2, 6, 9], dtype=float)
+
+    weights = norm.weights(x)
+    rho = norm.rho(x)  # not used
+    psi = norm.psi(x)
+    psi_deriv = norm.psi_deriv(x)
+
+    # check location and u-shape of rho
+    assert rho[4] == 0
+    assert np.all(np.diff(rho[4:]) >= 0)
+    assert np.all(np.diff(rho[:4]) <= 0)
+
+    # avoid zero division nan:
+    assert_allclose(weights, (psi + 1e-50) / (x + 1e-50), rtol=1e-6, atol=1e-8)
+    psid = _approx_fprime_scalar(x, norm.rho)
+    assert_allclose(psi, psid, rtol=1e-6, atol=1e-8)
+    psidd = _approx_fprime_scalar(x, norm.psi)
+    assert_allclose(psi_deriv, psidd, rtol=1e-6, atol=1e-8)