Kalman-and-Bayesian-Filters-in-Python/gauss.py

# -*- coding: utf-8 -*-
"""
Created on Tue Apr 22 11:38:49 2014

@author: rlabbe
"""
from __future__ import division, print_function
import math
import matplotlib.pyplot as plt
import numpy.random as random


class gaussian(object):
    def __init__ (self,m,s):
        self.mu    = float(m)
        self.sigma = float(s)

    def __add__ (a,b):
       return gaussian (a.mu + b.mu, a.sigma + b.sigma)

    def __mul__ (a,b):
        m = (a.sigma*b.mu + b.sigma*a.mu) / (a.sigma + b.sigma)
        s = 1. / (1./a.sigma + 1./b.sigma)
        return gaussian (m,s)

    def __call__(self,x):
        return math.exp (-0.5 * (x-self.mu)**2 / self.sigma) / \
        math.sqrt(2.*math.pi*self.sigma)

    def __str__(self):
        return "(%f, %f)" %(self.mu, self.sigma)

    def __getitem__ (self,index):
        """ maybe silly, allows you to access obect as a tuple:
        a = gaussian(3,4)
        print (tuple(a))
        """
        if index == 0:   return self.mu
        elif index == 1: return self.sigma
        else:            raise StopIteration

class KF1D(object):
    def __init__ (self, pos, sigma):
        self.estimate = gaussian(pos,sigma)

    def update(self, Z,var):
        self.estimate = self.estimate * gaussian (Z,var)

    def predict(self, U, var):
        self.estimate = self.estimate + gaussian (U,var)


measurements = [x+5 for x in range(100)]



def fixed_error_kf(measurement_error, noise_factor = 1.0):
    motion_sig = 2.
    mu = 0
    sig = 1000

    f = KF1D (mu,sig)

    ys = []
    errs = []
    xs = []

    for i in range(len(measurements)):
        r = random.randn() * noise_factor
        m = measurements[i] + r
        f.update (m, measurement_error)

        xs.append(m)
        ys.append(f.estimate.mu)
        errs.append (f.estimate.sigma)

        f.predict (1.0, motion_sig)

    plt.clf()

    plt.plot (measurements, 'r')
    plt.plot (xs,'g')
    plt.errorbar (x=range(len(ys)), color='b', y=ys, yerr=errs)
    plt.show()

def varying_error_kf(noise_factor=1.0):
    motion_sig = 2.
    mu = 0
    sig = 1000


    f = KF1D (mu,sig)
    ys = []
    us = []
    errs = []
    xs = []

    for i in range(len(measurements)):
        r = random.randn() * noise_factor
        m = measurements[i] + r
        print (r)
        f.update (m, abs(r*10))
        xs.append(m)
        #print ("measure:" + str(f.estimate))
        ys.append(f.estimate.mu)
        errs.append (f.estimate.sigma)

        f.predict (1.0, motion_sig)
        #print ("predict:" + str(f.estimate))

    plt.clf()
    plt.plot (measurements, 'r')
    plt.plot (xs,'g')
    plt.errorbar (x=range(len(ys)), color='b', y=ys, yerr=errs)
    plt.show()

varying_error_kf( noise_factor=100)