Repository: statsmodels/statsmodels (Python)

Issue:
BUG: incorrect definition of piecewise boolean indicator array in `Hampel.rho` (*robust* module)
#### Describe the bug

There is a bug in `class Hampel` within [statsmodels/robust/norms.py](https://github.com/statsmodels/statsmodels/blob/main/statsmodels/robust/norms.py) in defining piecewise function over `z` in `Hampel.rho` method.

It results in incorrect values for the `Hampel.rho` method outcomes.

The `rho` method uses `_subset` method to calculate boolean indicator arrays for `t1`, `t2`, and `t3` pieces, with the last piece `t4` defining the whole range needed for the Hampel's $\rho$ function being calculated in the `rho` method. This last piece is incorrect in the calculation of the boolean indicator array, I show in the code below.

Note that I can't reproduce the issue perfectly, since `Hampel.rho` does not rely on `Hampel._subset` entirely in defining piecewise boolean indicator arrays, as opposed to `Hampel.psi` or `Hampel.weights`.

#### Code Sample, a copy-pastable example if possible


```python
import statsmodels.api as sm
# create hampel norm with default parameters a, b, c
hampel = sm.robust.norms.Hampel(a=2., b=4., c=8.)
# create specially crafted t1, t2, t3
t1, t2, t3 = hampel._subset([1, 3, 5, 9])
# calculate t4 as in `Hampel.rho` (line 639 calculates that implicitly, without assigning variable `t4`)
t4 = (1 - t1 + t2 + t3)
# t4 gives `array([0, 2, 2, 1])`
```
<details>

**Note**: As you can see, there are many issues on our GitHub tracker, so it is very possible that your issue has been posted before. Please check first before submitting so that we do not have to handle and close duplicates.

**Note**: Please be sure you are using the latest released version of `statsmodels`, or a recent build of `main`. If your problem has been fixed in an unreleased version, you might be able to use `main` until a new release occurs. 

**Note**: If you are using a released version, have you verified that the bug exists in the main branch of this repository? It helps the limited resources if we know problems exist in the current main branch so that they do not need to check whether the code sample produces a bug in the next release.

</details>


If the issue has not been resolved, please file it in the issue tracker.

#### Expected Output

The implicit `t4` that is calculated "on the fly" in `Hampel.rho` should result in an array `array([0, 0, 0, 1])` to pick only the last piecewise boolean indicator array filter.

In fact, it is obvious the "on the fly" calculation in `Hampel.rho` line 639 should be

```python
(1 - (t1 + t2 + t3)) * a * (b + c - a))
```

rather than the existing

```python
(1 - t1 + t2 + t3) * a * (b + c - a))
```

#### Output of ``import statsmodels.api as sm; sm.show_versions()``

<details>

INSTALLED VERSIONS
------------------
Python: 3.8.10.final.0
OS: Linux 5.10.102.1-microsoft-standard-WSL2 #1 SMP Wed Mar 2 00:30:59 UTC 2022 x86_64
byteorder: little
LC_ALL: None
LANG: C.UTF-8

statsmodels
===========

Installed: 0.13.5 (/home/kamil/.cache/pypoetry/virtualenvs/smps-aJ6sAQBH-py3.8/lib/python3.8/site-packages/statsmodels)

Required Dependencies
=====================

cython: Not installed
numpy: 1.22.3 (/home/kamil/.cache/pypoetry/virtualenvs/smps-aJ6sAQBH-py3.8/lib/python3.8/site-packages/numpy)
scipy: 1.8.0 (/home/kamil/.cache/pypoetry/virtualenvs/smps-aJ6sAQBH-py3.8/lib/python3.8/site-packages/scipy)
pandas: 1.4.1 (/home/kamil/.cache/pypoetry/virtualenvs/smps-aJ6sAQBH-py3.8/lib/python3.8/site-packages/pandas)
    dateutil: 2.8.2 (/home/kamil/.cache/pypoetry/virtualenvs/smps-aJ6sAQBH-py3.8/lib/python3.8/site-packages/dateutil)
patsy: 0.5.3 (/home/kamil/.cache/pypoetry/virtualenvs/smps-aJ6sAQBH-py3.8/lib/python3.8/site-packages/patsy)

Optional Dependencies
=====================

matplotlib: 3.5.1 (/home/kamil/.cache/pypoetry/virtualenvs/smps-aJ6sAQBH-py3.8/lib/python3.8/site-packages/matplotlib)
    backend: module://matplotlib_inline.backend_inline 
cvxopt: Not installed
joblib: 1.1.0 (/home/kamil/.cache/pypoetry/virtualenvs/smps-aJ6sAQBH-py3.8/lib/python3.8/site-packages/joblib)

Developer Tools
================

IPython: 8.1.1 (/home/kamil/.cache/pypoetry/virtualenvs/smps-aJ6sAQBH-py3.8/lib/python3.8/site-packages/IPython)
    jinja2: 3.0.3 (/home/kamil/.cache/pypoetry/virtualenvs/smps-aJ6sAQBH-py3.8/lib/python3.8/site-packages/jinja2)
sphinx: Not installed
    pygments: 2.11.2 (/home/kamil/.cache/pypoetry/virtualenvs/smps-aJ6sAQBH-py3.8/lib/python3.8/site-packages/pygments)
pytest: Not installed
virtualenv: Not installed

</details>



Hints:
Note that this affects the same line of code as in #8697 but a different problem is raised in this BUG issue.

Patch:
diff --git a/statsmodels/robust/norms.py b/statsmodels/robust/norms.py
--- a/statsmodels/robust/norms.py
+++ b/statsmodels/robust/norms.py
@@ -3,6 +3,16 @@
 # 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 @@ def weights(self, z):
 
             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 @@ def psi_deriv(self, z):
         -----
         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 @@ def rho(self, z):
         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 @@ def psi(self, z):
         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 @@ def psi(self, z):
         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 @@ def weights(self, z):
         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 @@ def psi_deriv(self, z):
         """
 
         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 @@ def rho(self, z):
         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 @@ def rho(self, z):
         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 @@ def psi(self, z):
 
             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 @@ def weights(self, z):
         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 @@ def rho(self, z):
             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
--- /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
--- /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)