Small updates.

wafo/kdetools/kernels.py

@@ -9,8 +9,9 @@ import warnings
 import numpy as np
 from numpy import pi, sqrt, exp, percentile
 from numpy.fft import fft, ifft
-from scipy import optimize, linalg
+from scipy import optimize
 from scipy.special import gamma
+from scipy.linalg import sqrtm
 from wafo.misc import tranproc  # , trangood
 from wafo.kdetools.gridding import gridcount
 from wafo.dctpack import dct
@@ -289,8 +290,8 @@ class _Kernel(object):
     def kernel(self, x):
         return self._kernel(np.atleast_2d(x))
 
-    def deriv4_6_8_10(self, t, numout=4):
-        raise NotImplementedError('Method not implemented for this kernel!')
+#     def deriv4_6_8_10(self, t, numout=4):
+#         raise NotImplementedError('Method not implemented for this kernel!')
 
     def get_ste_constant(self, n):
         mu2, R = self.stats[:2]
@@ -542,7 +543,7 @@ class Kernel(object):
 
     """
 
-    def __init__(self, name, fun='hste'):  # 'hns'):
+    def __init__(self, name, fun='hste'):
         self.kernel = _MKERNEL_DICT[name[:4]]
         self.get_smoothing = getattr(self, fun)
 
@@ -681,7 +682,7 @@ class Kernel(object):
         """
         return self.hns(data) / 0.93
 
-    def _hmns_scale(self, d):
+    def _hmns_scale(self, n, d):
         name = self.name
         short_name = name[:4].lower()
         scale_dict = dict(epan=(8.0 * (d + 4.0) * (2 * sqrt(pi)) ** d /
@@ -694,7 +695,7 @@ class Kernel(object):
         if d > 2 and short_name in ['biwe', 'triw']:
             raise NotImplementedError('Not implemented for d>2 and '
                                       'kernel {}'.format(name))
-        return scale_dict[short_name]
+        return scale_dict[short_name] * n ** (-1. / (d + 4))
 
     def hmns(self, data):
         """Returns Multivariate Normal Scale Estimate of Smoothing Parameter.
@@ -742,9 +743,7 @@ class Kernel(object):
         d, n = a.shape
         if d == 1:
             return self.hns(data)
-        scale = self._hmns_scale(d)
-        cov_a = np.cov(a)
-        return scale * linalg.sqrtm(cov_a).real * n ** (-1. / (d + 4))
+        return self._hmns_scale(n, d) * np.real(sqrtm(np.cov(a)))
 
     @staticmethod
     def _get_g(k_order_2, mu2, psi_order, n, order):

wafo/kdetools/tests/test_gridding.py

@@ -19,6 +19,7 @@ class TestKdeTools(unittest.TestCase):
 
         dx = x[1] - x[0]
         c = wkg.gridcount(data, x)
+        assert_allclose(c.sum(), len(data))
         assert_allclose(c,
                         [0.1430937435034, 5.864465648665, 9.418694957317207,
                          2.9154367000439, 0.6583089504704, 0.0,
@@ -34,6 +35,7 @@ class TestKdeTools(unittest.TestCase):
         dx = x[1] - x[0]
         X = np.vstack((x, x))
         c = wkg.gridcount(data, X)
+        assert_allclose(c.sum(), N)
         assert_allclose(c,
                         [[0.38922806, 0.8987982,  0.34676493, 0.21042807,  0.],
                          [1.15012203, 5.16513541, 3.19250588, 0.55420752,  0.],
@@ -52,6 +54,7 @@ class TestKdeTools(unittest.TestCase):
         dx = x[1] - x[0]
         X = np.vstack((x, x, x))
         c = wkg.gridcount(data, X)
+        assert_allclose(c.sum(), N)
         assert_allclose(c,
                         [[[8.74229894e-01, 1.27910940e+00, 1.42033973e-01],
                           [1.94778915e+00, 2.59536282e+00, 3.28213680e-01],
@@ -69,14 +72,13 @@ class TestKdeTools(unittest.TestCase):
     @staticmethod
     def test_gridcount_4d():
 
-        N = 20
-        data = np.reshape(DATA2D, (4, -1))
+        N = 10
+        data = np.reshape(DATA2D, (4, N))
         x = np.linspace(0, max(np.ravel(data)) + 1, 3)
         dx = x[1] - x[0]
         X = np.vstack((x, x, x, x))
         c = wkg.gridcount(data, X)
-        assert_allclose(c,
-                        [[[[1.77163904e-01, 1.87720108e-01, 0.0],
+        truth = [[[[1.77163904e-01, 1.87720108e-01, 0.0],
                    [5.72573585e-01, 6.09557834e-01, 0.0],
                    [3.48549923e-03, 4.05931870e-02, 0.0]],
                   [[1.83770124e-01, 2.56357594e-01, 0.0],
@@ -102,10 +104,12 @@ class TestKdeTools(unittest.TestCase):
                      [6.61257395e-02, 2.47717418e-02, 0.0]],
                   [[6.38695629e-04, 5.69610302e-03, 0.0],
                      [1.00358265e-02, 2.44053065e-02, 0.0],
-                                [5.67244468e-03, 2.12498697e-03, 0.0]]]])
+                     [5.67244468e-03, 2.12498697e-03, 0.0]]]]
+        assert_allclose(c.sum(), N)
+        assert_allclose(c, truth)
 
         t = np.trapz(np.trapz(np.trapz(np.trapz(c / dx**4 / N, x), x), x), x)
-        assert_allclose(t, 0.21183518274521254)
+        assert_allclose(t, 0.4236703654904251)
 
 
 if __name__ == "__main__":