BUG: fix robust.norm.Hampel

Author: josef-pkt

statsmodels/robust/norms.py

@@ -629,15 +629,21 @@ def rho(self, z):
         """
 
         z = np.abs(z)
-        a = self.a
-        b = self.b
-        c = self.c
+        z_isscalar = np.isscalar(z)
+        z = np.atleast_1d(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))
-        return v
+        t34 = ~(t1 | t2 )
+        v = np.zeros(z.shape, float)
+        v[t1] = z[t1]**2 * 0.5
+        v[t2] = (a * (z[t2] - a) + 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:
+            return v[0]
+        else:
+            return v
 
     def psi(self, z):
         r"""
@@ -662,13 +668,11 @@ def psi(self, z):
             psi(z) = 0                            for \|z\| > c
         """
         z = np.asarray(z)
-        a = self.a
-        b = self.b
-        c = self.c
+        a = self.a; b = self.b; c = self.c
         t1, t2, t3 = self._subset(z)
         s = np.sign(z)
         z = np.abs(z)
-        v = s * (t1 * z +
+        v =  (t1 * z*s +
                  t2 * a*s +
                  t3 * a*s * (c - z) / (c - b))
         return v
@@ -711,13 +715,16 @@ def weights(self, z):
         return v
 
     def psi_deriv(self, z):
+        """Derivative of psi function, second derivative of rho function.
+        """
         t1, _, t3 = self._subset(z)
         a, b, c = self.a, self.b, self.c
+        z = np.atleast_1d(z)
         # default is t1
         d = np.zeros_like(z)
         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))
         return d
 
 

statsmodels/robust/tests/results/results_norms.py

@@ -0,0 +1,14 @@
+
+
+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],
+    )

statsmodels/robust/tests/test_norms.py

@@ -0,0 +1,22 @@
+
+import pytest
+import numpy as np
+from numpy.testing import assert_allclose
+
+from statsmodels.robust import norms
+from .results import results_norms as res_r
+
+cases = [
+    (norms.Hampel, (1.5, 3.5, 8.), res_r.res_hampel)
+    ]
+
+@pytest.mark.parametrize("case", cases)
+def test_norm(case):
+    ncls, args, res = case
+    norm = ncls(*args)
+    x = np.array([-9., -6, -2, -1, 0, 1, 2, 6, 9])
+
+    assert_allclose(norm.weights(x), res.weights, rtol=1e-12, atol=1e-20)
+    assert_allclose(norm.rho(x), res.rho, rtol=1e-12, atol=1e-20)
+    assert_allclose(norm.psi(x), res.psi, rtol=1e-12, atol=1e-20)
+    assert_allclose(norm.psi_deriv(x), res.psi_deriv, rtol=1e-12, atol=1e-20)