Added section on fat tails.

exp/distributions.py

@@ -0,0 +1,42 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Sat Jan 17 20:17:14 2015
+
+@author: rlabbe
+"""
+
+from scipy.stats import t, norm
+import matplotlib.pyplot as plt
+import numpy as np
+import math
+import random
+
+
+df =4
+mu = 10
+std = 2
+
+x = np.linspace(-5, 20, 100)
+               
+plt.plot(x, t.pdf(x, df=df, loc=mu, scale=std), 'r-', lw=5, label='t pdf')
+
+
+x2 = np.linspace(mu-10, mu+10, 100)
+plt.plot(x, norm.pdf(x, mu, std), 'b-', lw=5,  label='gaussian pdf')
+plt.legend()
+plt.figure()
+
+def student_t(df, mu, std): # nu equals number of degrees of freedom
+    x = random.gauss(0, std)
+    y = 2.0*random.gammavariate(0.5*df, 2.0)
+    return x / (math.sqrt(y/df)) + mu
+    
+    
+N = 100000
+ys = [student_t(2.7, 100, 2) for i in range(N)]
+plt.hist(ys, 10000, histtype='step')
+
+ys = [random.gauss(100,2) for i in range(N)]
+plt.hist(ys, 10000, histtype='step',color='r')
+
+plt.show()

exp/doppler.py

@@ -0,0 +1,140 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Sat Jan 17 15:04:08 2015
+
+@author: rlabbe
+"""
+
+
+
+from numpy.random import randn
+import numpy as np
+from numpy import array
+import matplotlib.pyplot as plt
+from filterpy.common import Q_discrete_white_noise
+from filterpy.kalman import KalmanFilter
+from math import sqrt
+import math
+import random
+
+pos_var = 1.
+dop_var = 2.
+dt = 1/20
+
+
+def rand_student_t(df, mu=0, std=1):
+    """return random number distributed by student's t distribution with
+    `df` degrees of freedom with the specified mean and standard deviation.
+    """
+    x = random.gauss(0, std)
+    y = 2.0*random.gammavariate(0.5*df, 2.0)
+    return x / (math.sqrt(y/df)) + mu
+    
+
+
+np.random.seed(124)
+class ConstantVelocityObject(object):
+    def __init__(self, x0, vel, noise_scale):
+        self.x = np.array([x0, vel])
+        self.noise_scale = noise_scale
+        self.vel = vel
+
+
+    def update(self, dt):
+        pnoise = abs(randn()*self.noise_scale)
+        if self.x[1] > self.vel:
+            pnoise = -pnoise
+            
+        self.x[1] += pnoise
+        self.x[0] += self.x[1]*dt
+
+        return self.x.copy()
+        
+        
+    def sense_pos(self):
+        return self.x[0] + [randn()*self.noise_scale[0]]
+
+        
+    def vel(self):
+        return self.x[1] + [randn()*self.noise_scale[1]]
+
+
+def sense(x, noise_scale=(1., 0.01)):
+    return x + [randn()*noise_scale[0], randn()*noise_scale[1]]
+
+def sense_t(x, noise_scale=(1., 0.01)):
+    return x + [rand_student_t(2)*noise_scale[0], 
+                rand_student_t(2)*noise_scale[1]]
+
+
+
+
+
+R2 = np.zeros((2, 2))
+R1 = np.zeros((1, 1))
+
+R2[0, 0] = pos_var
+R2[1, 1] = dop_var
+
+R1[0,0] = dop_var
+
+H2 = np.array([[1., 0.], [0., 1.]])
+H1 = np.array([[0., 1.]])
+
+
+
+kf = KalmanFilter(2, 1)
+kf.F = array([[1, dt], 
+              [0, 1]])
+kf.H = array([[1, 0], 
+              [0, 1]])
+kf.Q = Q_discrete_white_noise(2, dt, .02)
+kf.x = array([0., 0.])
+kf.R = R1
+kf.H = H1
+kf.P *= 10
+
+
+random.seed(1234)
+track = []
+zs = []
+ests = []
+sensor_noise = (sqrt(pos_var), sqrt(dop_var))
+obj = ConstantVelocityObject(0., 1., noise_scale=.18)
+
+for i in range(30000):
+    x = obj.update(dt/10)
+    z = sense_t(x, sensor_noise)
+    if i % 10 != 0:
+        continue
+
+    if i % 60 == 87687658760:
+        kf.predict()
+        kf.update(z, H=H2, R=R2)        
+    else:
+        kf.predict()
+        kf.update(z[1], H=H1, R=R1) 
+        
+    ests.append(kf.x)
+        
+    track.append(x)
+    zs.append(z)
+
+
+track = np.asarray(track)
+zs = np.asarray(zs)
+ests = np.asarray(ests)
+
+plt.figure()
+plt.subplot(211)
+plt.plot(track[:,0], label='track')
+plt.plot(zs[:,0], label='measurements')
+plt.plot(ests[:,0], label='filter')
+plt.legend(loc='best')
+
+plt.subplot(212)
+plt.plot(track[:,1], label='track')
+plt.plot(zs[:,1], label='measurements')
+plt.plot(ests[:,1], label='filter')
+plt.show()
+