Added several kernel functions and tests

KDEpy/kde.py

@@ -61,7 +61,7 @@ def fit(self, data, boundaries=None):
 
     def _set_weights(self, weights):
         if weights is None:
-            weights = np.ones_like(self._data)
+            weights = np.asfarray(np.ones_like(self._data))
         weights = weights / np.sum(weights)
         return weights
 
@@ -78,7 +78,9 @@ def evaluate_naive(self, grid_points, weights=None):
         weights : np.array of weights for the data points, must sum to unity
         
         """
-        grid_points = grid_points.astype(float)
+        # Return the array converted to a float type
+        grid_points = np.asfarray(grid_points)
+
         # If no weights are passed, weight each data point as unity
         weights = self._set_weights(weights)
 
@@ -103,7 +105,7 @@ def _eval_sorted(self, data_sorted, weights_sorted, grid_point, len_data):
 
         # The relationship between the desired bandwidth and the original
         # bandwidth of the kernel function
-        bw_scale = self.bw / abs(self.kernel.right_bw + self.kernel.left_bw)
+        bw_scale = self.bw / abs(self.kernel.right_bw - self.kernel.left_bw)
 
         # Compute the bandwidth to the left and to the right of the center
         left_bw = self.kernel.left_bw * bw_scale

KDEpy/kernel_funcs.py

@@ -8,6 +8,15 @@
 
 import numpy as np
 import collections.abc
+import numbers
+
+# In R, the following are implemented:
+# "gaussian", "rectangular", "triangular", "epanechnikov", 
+# "biweight", "cosine" or "optcosine"
+
+# Wikipedia
+# uniform, trinagular, epanechnikov, quartic,triweight, tricube,
+# gaussian, cosine, logistic, sigmoid, silverman
 
 
 def epanechnikov(x):
@@ -35,9 +44,45 @@ def tri(x):
     return out
 
 
+def biweight(x):
+    out = np.zeros_like(x)
+    mask = np.logical_and((x < 1), (x > -1))
+    out[mask] = ((15 / 16) * (1 - x**2)**2)[mask]
+    return out
+
+
+def triweight(x):
+    out = np.zeros_like(x)
+    mask = np.logical_and((x < 1), (x > -1))
+    out[mask] = ((35 / 32) * (1 - x**2)**3)[mask]
+    return out
+
+
+def tricube(x):
+    out = np.zeros_like(x)
+    mask = np.logical_and((x < 1), (x > -1))
+    out[mask] = ((70 / 81) * (1 - np.abs(x)**3)**3)[mask]
+    return out
+
+
+def cosine(x):
+    out = np.zeros_like(x)
+    mask = np.logical_and((x < 1), (x > -1))
+    out[mask] = ((np.pi / 4) * np.cos((np.pi * x) / 2))[mask]
+    return out
+
+
+def logistic(x):
+    return 1 / (2 + 2 * np.cosh(x))
+
+
+def sigmoid(x):
+    return (1 / (np.pi * np.cosh(x)))
+
+
 class Kernel(collections.abc.Callable):
 
-    def __init__(self, function, expected_value=0, left_bw=1, right_bw=1):
+    def __init__(self, function, var=1, support=(-3, 3)):
         """
         Initialize a new kernel function.
         
@@ -47,27 +92,54 @@ def __init__(self, function, expected_value=0, left_bw=1, right_bw=1):
         left_bw: support to the right
         """
         self.function = function
-        self.expected_value = expected_value
-        self.left_bw = left_bw
-        self.right_bw = right_bw
+        self.var = var
+        self.support = support
+        self.finite_support = np.all(np.isfinite(np.array(self.support)))
+        assert self.support[0] < self.support[1]
 
     def evaluate(self, x, bw=1):
         """
         Evaluate the kernel.
         """
-        real_bw = (bw / (self.left_bw + self.right_bw))
+
+        # If x is a number, convert it to a length-1 NumPy vector
+        if isinstance(x, numbers.Number):
+            x = np.asarray_chkfinite([x])
+        else:
+            x = np.asarray_chkfinite(x)
+
+        # Scale the function
+        real_bw = bw / np.sqrt(self.var)
         return self.function(x / real_bw) / real_bw
 
     def __call__(self, *args, **kwargs):
         return self.evaluate(*args, **kwargs)
 
 
-gaussian = Kernel(gaussian, 0, 3, 3)
-box = Kernel(box, 0, 1, 1)
-tri = Kernel(tri, 0, 1, 1)
-epa = Kernel(epanechnikov, 0, 1, 1)
+gaussian = Kernel(gaussian, var=1, support=(-np.inf, np.inf))
+box = Kernel(box, var=1 / 3, support=(-1, 1))
+tri = Kernel(tri, var=1 / 6, support=(-1, 1))
+epa = Kernel(epanechnikov, var=1 / 5, support=(-1, 1))
+biweight = Kernel(biweight, var=1 / 7, support=(-1, 1))
+triweight = Kernel(triweight, var=1 / 9, support=(-1, 1))
+tricube = Kernel(tricube, var=35 / 243, support=(-1, 1))
+cosine = Kernel(cosine, var=(1 - (8 / np.pi**2)), support=(-1, 1))
+logistic = Kernel(logistic, var=(np.pi**2 / 3), support=(-np.inf, np.inf))
+sigmoid = Kernel(sigmoid, var=(np.pi**2 / 4), support=(-np.inf, np.inf))
 
 _kernel_functions = {'gaussian': gaussian,
                      'box': box,
                      'tri': tri,
-                     'epa': epa}
+                     'epa': epa,
+                     'biweight': biweight,
+                     'triweight': triweight,
+                     'tricube': tricube,
+                     'cosine': cosine,
+                     'logistic': logistic,
+                     'sigmoid': sigmoid}
+
+
+if __name__ == "__main__":
+    import pytest
+    # --durations=10  <- May be used to show potentially slow tests
+    pytest.main(args=['.', '--doctest-modules', '-v'])

KDEpy/tests/test_kde.py

@@ -0,0 +1,56 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Tests.
+"""
+import numpy as np
+from KDEpy import KDE
+import pytest
+
+
+class TestNaiveKDE():
+
+    @pytest.mark.parametrize("kernel, bw, n, expected_result", 
+                             [('box', 0.1, 5, np.array([2.101278e-19, 
+                                                        3.469447e-18, 
+                                                        1.924501e+00, 
+                                                        0.000000e+00, 
+                                                        9.622504e-01])),
+                              ('box', 0.2, 5, np.array([3.854941e-18, 
+                                                        2.929755e-17,
+                                                        9.622504e-01,
+                                                        0.000000e+00,
+                                                        4.811252e-01])),
+                              ('box', 0.6, 3, np.array([0.1603751,
+                                                        0.4811252,
+                                                        0.4811252])),
+                              ('tri', 0.6, 3, np.array([0.1298519,
+                                                        0.5098009,
+                                                        0.3865535])),
+                              ('epa', 0.1, 6, np.array([0.000000e+00,
+                                                        7.285839e-17,
+                                                        2.251871e-01,
+                                                        1.119926e+00,
+                                                        0.000000e+00,
+                                                        1.118034e+00])),
+                              ('biweight', 2, 5, np.array([0.1524078,
+                                                           0.1655184,
+                                                           0.1729870,
+                                                           0.1743973,
+                                                           0.1696706]))])
+    def test_against_R_density(self, kernel, bw, n, expected_result):
+        """
+        Test against the following function call in R:
+            d <- density(c(0, 0.1, 1), kernel="{kernel}", bw={bw}, 
+            n={n}, from=-1, to=1);
+            d$y
+        """
+        data = np.array([0, 0.1, 1])
+        x = np.linspace(-1, 1, num=n)
+        y = KDE(kernel, bw=bw).fit(data).evaluate_naive(x)
+        assert np.allclose(y, expected_result, atol=10**(-2.7))
+
+
+if __name__ == "__main__":
+    # --durations=10  <- May be used to show potentially slow tests
+    pytest.main(args=['.', '--doctest-modules', '-v'])

KDEpy/tests/test_kernel_funcs.py

@@ -0,0 +1,70 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Tests.
+"""
+import numpy as np
+from scipy.integrate import quad
+from KDEpy import KDE
+import pytest
+
+
+class TestKernelFunctions():
+
+    @pytest.mark.parametrize("fname, function", 
+                             list(KDE._available_kernels.items()))
+    def test_integral_unity(self, fname, function):
+        """
+        Verify that all available kernel functions have an integral evaluating
+        to unity. This is a requirement of the kernel functions.
+        """
+
+        if function.finite_support:
+            # When bw=1, the function is scaled by the standard deviation
+            # so that std(f) = 1. Divide by standard deviation to get
+            # the integration limits.
+            a, b = tuple(s / np.sqrt(function.var) for s in function.support)
+        else:
+            a, b = -30, 30
+        integral, abserr = quad(function, a=a, b=b)
+        assert np.isclose(integral, 1)
+
+    @pytest.mark.parametrize("fname, function", 
+                             list(KDE._available_kernels.items()))
+    def test_monotonic_decreasing(self, fname, function):
+        """
+        Verify that all available kernel functions decrease away from their 
+        mode.
+        """
+
+        if function.finite_support:
+            a, b = tuple(s / np.sqrt(function.var) for s in function.support)
+            x = np.linspace(a, b)
+        else:
+            x = np.linspace(-50, 50)
+        y = function(x)
+        diffs_left = np.diff(y[x <= 0])
+        diffs_right = np.diff(y[x >= 0])
+        assert np.all(diffs_right <= 0)
+        assert np.all(diffs_left >= 0)
+
+    @pytest.mark.parametrize("fname, function", 
+                             list(KDE._available_kernels.items()))
+    def test_non_negative(self, fname, function):
+        """
+        Verify that all available kernel functions are non-negative.
+        """
+
+        if function.finite_support:
+            a, b = tuple(s / np.sqrt(function.var) for s in function.support)
+            x = np.linspace(a, b)
+        else:
+            x = np.linspace(-20, 20)
+        y = function(x)
+        assert np.all(y >= 0)
+
+
+if __name__ == "__main__":
+    import pytest
+    # --durations=10  <- May be used to show potentially slow tests
+    pytest.main(args=['.', '--doctest-modules', '-v'])

docs/source/_static/custom.css

@@ -0,0 +1,3 @@
+.prompt {
+    display: none;
+}

docs/source/conf.py

@@ -47,7 +47,8 @@
     'sphinx.ext.napoleon',
     'sphinx.ext.autosummary',
     'sphinx.ext.extlinks',
-    'numpydoc'
+    'numpydoc',
+    'nbsphinx'
 ]
 
 # Add any paths that contain templates here, relative to this directory.
@@ -72,7 +73,7 @@
 # List of patterns, relative to source directory, that match files and
 # directories to ignore when looking for source files.
 # This pattern also affects html_static_path and html_extra_path .
-exclude_patterns = []
+exclude_patterns = ['_build', '**.ipynb_checkpoints']
 
 # The name of the Pygments (syntax highlighting) style to use.
 pygments_style = 'sphinx'
@@ -83,13 +84,28 @@
 # The theme to use for HTML and HTML Help pages.  See the documentation for
 # a list of builtin themes.
 #
-html_theme = 'classic'
+html_theme = 'alabaster'
 
 # Theme options are theme-specific and customize the look and feel of a theme
 # further.  For a list of options available for each theme, see the
 # documentation.
 #
-# html_theme_options = {}
+html_theme_options = {
+    # 'logo': 'logo.png',
+    'github_user': 'tommyod',
+    'github_repo': 'KDEpy',
+    'github_button': True,
+    'github_banner': True,
+    'travis_button': True,
+    'show_powered_by': False,
+    'font_family': '-apple-system,BlinkMacSystemFont,"Segoe UI",Roboto,\
+        "Helvetica Neue",Arial,sans-serif,"Apple Color Emoji","Segoe UI Emoji"\
+        ,"Segoe UI Symbol"',
+    'font_size': '15px',
+    'head_font_family': '-apple-system,BlinkMacSystemFont,"Segoe UI",Roboto,\
+        "Helvetica Neue",Arial,sans-serif,"Apple Color Emoji","Segoe UI Emoji"\
+        ,"Segoe UI Symbol"'
+}
 
 # Add any paths that contain custom static files (such as style sheets) here,
 # relative to this directory. They are copied after the builtin static files,
@@ -106,6 +122,8 @@
 #
 # html_sidebars = {}
 
+html_show_sphinx = False
+
 
 # -- Options for HTMLHelp output ---------------------------------------------
 

docs/source/index.rst

@@ -52,6 +52,9 @@ Doctest example:
    :maxdepth: 2
    :caption: Contents:
 
+   intro_kde
+   notebook.ipynb
+