BUG/REF/TST: add consistency checks, fix smaller bugs

Author: josef-pkt

statsmodels/robust/norms.py

@@ -255,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):
         """
@@ -269,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?
@@ -408,7 +414,7 @@ def rho(self, z):
         a = self.a
         z = np.asarray(z)
         test = self._subset(z)
-        return (test * a * (1 - np.cos(z / a)) +
+        return (test * a**2 * (0 - np.cos(z / a)) +
                 (1 - test) * 2 * a)
 
     def psi(self, z):
@@ -433,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"""
@@ -477,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

statsmodels/robust/tests/results/results_norms.py

@@ -22,3 +22,12 @@
                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],
+    )

statsmodels/robust/tests/test_norms.py

@@ -13,6 +13,15 @@
 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]
@@ -22,6 +31,10 @@
 @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)
 
@@ -56,3 +69,23 @@ def test_norm(case, dtype):
     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)
+
+    # 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)