Added discrete_bayes module.

Added module providing discrete Bayes filtering functions.

Includes documentation.

I also added the napolean extension to sphinx; this allows me to
use the NumPy style docstrings more completely.

docs/conf.py

@@ -21,7 +21,8 @@
 MOCK_MODULES = ['numpy', 'scipy', 'matplotlib', 'matplotlib.pyplot',
                 'scipy.linalg', 'numpy.linalg', 'matplotlib.pyplot',
                 'numpy.random', 'scipy.sparse', 'scipy.sparse.linalg',
-                'scipy.stats', 'matplotlib.patches']
+                'scipy.stats', 'matplotlib.patches', 'scipy.ndimage.filters',
+                'scipy.ndimage.interpolation']
 
 for mod_name in MOCK_MODULES:
     sys.modules[mod_name] = mock.Mock()
@@ -51,7 +52,8 @@
     'sphinx.ext.intersphinx',
     'sphinx.ext.mathjax',
     'sphinx.ext.viewcode',
-    'sphinx.ext.autodoc'
+    'sphinx.ext.autodoc',
+    'sphinx.ext.napoleon'
 ]
 
 # Add any paths that contain templates here, relative to this directory.

docs/discrete_bayes/discrete_bayes.rst

@@ -0,0 +1,20 @@
+Discrete Bayes
+==============
+
+words go here
+
+.. automodule:: filterpy.discrete_bayes
+
+-----
+
+.. autofunction:: normalize
+
+-----
+
+.. autofunction:: update
+
+-----
+
+.. autofunction:: predict
+
+-----

docs/index.rst

@@ -248,6 +248,15 @@ particle filtering.
 
    monte_carlo/resampling
 
+filterpy.discrete_bayes
+-----------------------
+
+Routines for performing discrete Bayes filtering.
+
+.. toctree::
+   :maxdepth: 1
+
+   discrete_bayes/discrete_bayes
 
 filterpy.gh
 -----------

filterpy/discrete_bayes/__init__.py

@@ -0,0 +1,20 @@
+# -*- coding: utf-8 -*-
+"""Copyright 2015 Roger R Labbe Jr.
+
+FilterPy library.
+http://github.com/rlabbe/filterpy
+
+Documentation at:
+https://filterpy.readthedocs.org
+
+Supporting book at:
+https://github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python
+
+This is licensed under an MIT license. See the readme.MD file
+for more information.
+"""
+
+from __future__ import (absolute_import, division, print_function,
+                        unicode_literals)
+__all__=["discrete_bayes"]
+from .discrete_bayes import *

filterpy/discrete_bayes/discrete_bayes.py

@@ -0,0 +1,128 @@
+# -*- coding: utf-8 -*-
+"""Copyright 2015 Roger R Labbe Jr.
+
+FilterPy library.
+http://github.com/rlabbe/filterpy
+
+Documentation at:
+https://filterpy.readthedocs.org
+
+Supporting book at:
+https://github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python
+
+This is licensed under an MIT license. See the readme.MD file
+for more information.
+"""
+
+
+from __future__ import (absolute_import, division, print_function,
+                        unicode_literals)
+
+import numpy as np
+from scipy.ndimage.filters import convolve
+from scipy.ndimage.interpolation import shift
+
+
+def normalize(pdf):
+    """Normalize distribution `pdf` in-place so it sums to 1.0.
+
+    Returns pdf for convienence, so you can write things like:
+
+    >>> kernel = normalize(randn(7))
+
+    Parameters
+    ----------
+
+    pdf : ndarray
+        discrete distribution that needs to be converted to a pdf. Converted
+        in-place, i.e., this is modified.
+
+    Returns
+    -------
+
+    pdf : ndarray
+        The converted pdf.
+    """
+
+    pdf /= sum(np.asarray(pdf, dtype=float))
+    return pdf
+
+
+def update(likelihood, prior):
+    """ Computes the posterior of a discrete random variable given a
+    discrete likelihood and prior. In a typical application the likelihood
+    will be the likelihood of a measurement matching your current environment,
+    and the prior comes from discrete_bayes.predict().
+
+    Parameters
+    ----------
+
+    likelihood : ndarray, dtype=flaot
+         array of likelihood values
+
+    prior : ndarray, dtype=flaot
+        prior pdf.
+
+    Returns
+    -------
+
+    posterior : ndarray, dtype=float
+        Returns array representing the posterior.
+
+
+    Example
+    -------
+    .. code::
+
+        # self driving car. Sensor returns values that can be equated to positions
+        # on the road. A real likelihood compuation would be much more complicated
+        # than this example.
+
+        likelihood = np.ones(len(road))
+        likelihood[road==z] *= scale_factor
+
+        prior = predict(posterior, velocity, kernel)
+        posterior = update(likelihood, prior)
+    """
+
+    posterior = prior * likelihood
+    return normalize(posterior)
+
+
+
+def predict(pdf, offset, kernel, mode='wrap', cval=0.):
+    """ Performs the discrete Bayes filter prediction step, generating
+    the prior.
+
+    `pdf` is a discrete probability distribution expressing our initial
+    belief.
+
+    `offset` is an integer specifying how much we want to move to the right
+    (negative values means move to the left)
+
+    We assume there is some noise in that offset, which we express in `kernel`.
+    For example, if offset=3 and kernel=[.1, .7., .2], that means we think
+    there is a 70% chance of moving right by 3, a 10% chance of moving 2
+    spaces, and a 20% chance of moving by 4.
+
+    It returns the resulting distribution.
+
+    If `mode='wrap'`, then the probability distribution is wrapped around
+    the array.
+
+    If `mode='constant'`, or any other value the pdf is shifted, with `cval`
+    used to fill in missing elements.
+
+    Example
+    -------
+    .. code::
+
+        belief = [.05, .05, .05, .05, .55, .05, .05, .05, .05, .05]
+        prior = predict(belief, offset=2, kernel=[.1, .8, .1])
+    """
+
+    if mode == 'wrap':
+        return convolve(np.roll(pdf, offset), kernel, mode='wrap')
+    else:
+        return convolve(shift(pdf, offset, cval=cval), kernel,
+                        cval=cval, mode='constant')

filterpy/discrete_bayes/tests/test_discrete_bayes.py

@@ -0,0 +1,47 @@
+# -*- coding: utf-8 -*-
+"""Copyright 2015 Roger R Labbe Jr.
+
+FilterPy library.
+http://github.com/rlabbe/filterpy
+
+Documentation at:
+https://filterpy.readthedocs.org
+
+Supporting book at:
+https://github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python
+
+This is licensed under an MIT license. See the readme.MD file
+for more information.
+"""
+
+from filterpy.discrete_bayes import predict, update, normalize
+from numpy.random import randn, randint
+import numpy as np
+
+
+def _predict(distribution, offset, kernel):
+    """ explicit convolution with wraparound"""
+
+    N = len(distribution)
+    kN = len(kernel)
+    width = int((kN - 1) / 2)
+
+    prior = np.zeros(N)
+    for i in range(N):
+        for k in range (kN):
+            index = (i + (width-k) - offset) % N
+            prior[i] += distribution[index] * kernel[k]
+    return prior
+
+
+def test_predictions():
+    s = 0.
+
+    for k in range(3, 22, 2):  # different kernel sizes
+        for _ in range(1000):
+            a = randn(100)
+            kernel = normalize(randn(k))
+            move = randint(1, 200)
+            s += sum(predict(a, move, kernel) - _predict(a, move, kernel))
+
+        assert s < 1.e-8, "sum of difference = {}".format(s)

filterpy/kalman/kalman_filter.py

@@ -361,21 +361,22 @@ def test_matrix_dimensions(self, z=None, H=None, R=None, F=None, Q=None):
         Hx = dot(H, x)
 
         if z_shape == (): # scalar or np.array(scalar)
-            assert Hx.ndim == 1 or shape(Hx) == (1,1), "shape of z should be {}, not {}".format(
+            assert Hx.ndim == 1 or shape(Hx) == (1,1), \
+            "shape of z should be {}, not {} for the given H".format(
                    shape(Hx), z_shape)
 
         elif shape(Hx) == (1,):
-            assert z_shape[0] == 1, 'Shape of z must be {}'.format(shape(Hx))
+            assert z_shape[0] == 1, 'Shape of z must be {} for the given H'.format(shape(Hx))
 
         else:
             assert (z_shape == shape(Hx) or
                     (len(z_shape) == 1 and shape(Hx) == (z_shape[0], 1))), \
-                    "shape of z should be {}, not {}".format(
+                    "shape of z should be {}, not {} for the given H".format(
                     shape(Hx), z_shape)
 
         if np.ndim(Hx) > 1 and shape(Hx) != (1,1):
             assert shape(Hx) == z_shape, \
-               'shape of z should be {}, but it is {}'.format(
+               'shape of z should be {} for the given H, but it is {}'.format(
                shape(Hx), z_shape)
 
 

filterpy/kalman/sigma_points.py

@@ -125,6 +125,8 @@ def sigma_points(self, x, P):
 
         if  np.isscalar(P):
             P = np.eye(n)*P
+        else:
+            P = np.asarray(P)
 
         lambda_ = self.alpha**2 * (n + self.kappa) - n
         U = self.sqrt((lambda_ + n)*P)

filterpy/kalman/tests/test_ekf.py

@@ -67,6 +67,10 @@ def hx(x):
 
 
 
+def fx(x, dt):
+    return np.dot(rk.F, x)
+
+
 rk.x = array([-10., 90., 1100.])
 rk.R *= 10
 rk.Q = array([[0, 0, 0],