Skip to content

Commit

Permalink
Docstring, variable name, and bug fixes
Browse files Browse the repository at this point in the history
Not a 'clean' commit because I forgot that I had put imm.py into
.gitignore.

I added a lot of docstrings to the file.

I renamed variables to be closer to the (inconsistent) names used
in the literature.
  • Loading branch information
rlabbe committed Aug 13, 2015
1 parent f2322ff commit b63f69b
Show file tree
Hide file tree
Showing 3 changed files with 55 additions and 50 deletions.
1 change: 0 additions & 1 deletion .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,5 @@ dist/
_build/
_static
_templates
imm.py
quadrature.py
.settings/
90 changes: 50 additions & 40 deletions filterpy/kalman/imm.py
Original file line number Diff line number Diff line change
Expand Up @@ -23,13 +23,44 @@

class IMMEstimator(object):
""" Implements an Interacting Multiple-Model (IMM) estimator.
**References**
Bar-Shalom, Y., Li, X-R., and Kirubarajan, T. "Estimation with
Application to Tracking and Navigation". Wiley-Interscience, 2001.
Crassidis, J and Junkins, J. "Optimal Estimation of
Dynamic Systems". CRC Press, second edition. 2012.
Labbe, R. "Kalman and Bayesian Filters in Python".
https://github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python
"""

def __init__(self, x_shape, filters, p, trans):
def __init__(self, filters, mu, M):
""""
**Parameters**
filters : (N,) array_like of KalmanFilter objects
List of N filters. filters[i] is the ith Kalman filter in the
IMM estimator.
mu : (N,N) ndarray of float
mode probability: mu[i] is the probability that
filter i is the correct one.
M : (N,N) ndarray of float
Markov chain transition matrix. M[i,j] is the probability of
switching from filter i to filter j.
"""

assert len(filters) > 1

self.filters = filters
self.w = p
self.trans = trans
self.mu = mu
self.M = M

x_shape = filters[0].x.shape
try:
self.N = x_shape[0]
except:
Expand All @@ -38,24 +69,11 @@ def __init__(self, x_shape, filters, p, trans):
self.x = np.zeros(x_shape)
self.P = np.zeros((self.N, self.N))

self.cbar = dot(self.trans.T, self.w)
self.cbar = dot(self.M.T, self.mu)
self.n = len(filters)



'''@property
def x(self):
""" The estimated state of the bank of filters."""
return self._x
@property
def P(self):
""" Estimated covariance of the bank of filters."""
return self._P'''



def update(self, z, R=None, H=None):
def update(self, z):
"""
Add a new measurement (z) to the Kalman filter. If z is None, nothing
is changed.
Expand All @@ -64,59 +82,51 @@ def update(self, z, R=None, H=None):
z : np.array
measurement for this update.
R : np.array, scalar, or None
Optionally provide R to override the measurement noise for this
one call, otherwise self.R will be used.
H : np.array, or None
Optionally provide H to override the measurement function for this
one call, otherwise self.H will be used.
"""

L = zeros(len(self.filters))
for i, f in enumerate(self.filters):
f.update(z, R, H)
f.update(z)
L[i] = f.likelihood # prior * likelihood

# compute weights
self.w = self.cbar * L
self.w /= sum(self.w) # normalize
# compute mode probabilities for this step
self.mu = self.cbar * L
self.mu /= sum(self.mu) # normalize


# compute IMM state and covariance
# compute mixed IMM state and covariance
self.x.fill(0.)
self.P.fill(0.)

for f, w in zip(self.filters, self.w):
for f, w in zip(self.filters, self.mu):
self.x += f.x*w

for f, w in zip(self.filters, self.w):
for f, w in zip(self.filters, self.mu):
y = f.x - self.x
self.P += w*(np.outer(y, y) + f.P)


# initial condition IMM state, covariance
xs, Ps = [], []
self.cbar = dot(self.trans.T, self.w)
self.cbar = dot(self.M.T, self.mu)

omega = np.zeros((self.n, self.n))
for i in range(self.n):
omega[i, 0] = self.trans[i, 0]*self.w[i] / self.cbar[0]
omega[i, 1] = self.trans[i, 1]*self.w[i] / self.cbar[1]
omega[i, 0] = self.M[i, 0]*self.mu[i] / self.cbar[0]
omega[i, 1] = self.M[i, 1]*self.mu[i] / self.cbar[1]


# compute mixed states
# compute initial states
for i, (f, w) in enumerate(zip(self.filters, omega.T)):
x = np.zeros(self.x.shape)
P = np.zeros(self.P.shape)
for kf, wj in zip(self.filters, w):
x += kf.x * wj
xs.append(x)

P = np.zeros(self.P.shape)
for kf, wj in zip(self.filters, w):
y = kf.x - x
P += wj*(np.outer(y, y) + kf.P)
xs.append(x)
Ps.append(P)

for i in range(len(xs)):
Expand Down
14 changes: 5 additions & 9 deletions filterpy/kalman/tests/test_imm.py
Original file line number Diff line number Diff line change
Expand Up @@ -155,9 +155,7 @@ def test_imm():
Second edition.
"""



r = 100.
r = 1.
dt = 1.
phi_sim = np.array(
[[1, dt, 0, 0],
Expand All @@ -171,17 +169,16 @@ def test_imm():
[0, dt]])

x = np.array([[2000, 0, 10000, -15.]]).T
#x = np.genfromtxt('c:/users/rlabbe/dropbox/Crassidis/x.csv', delimiter=',')

simxs = []
N = 600
for i in range(N):

x = np.dot(phi_sim, x)
if i >= 400:
x += np.dot(gam, np.array([[.075, .075]]).T)
simxs.append(x)
simxs = np.array(simxs)
#x = np.genfromtxt('c:/users/rlabbe/dropbox/Crassidis/mycode/x.csv', delimiter=',')

zs = np.zeros((N, 2))
for i in range(len(zs)):
Expand All @@ -190,7 +187,7 @@ def test_imm():

try:
#data to test against crassidis' IMM matlab code
zs_tmp = np.genfromtxt('c:/users/rlabbe/dropbox/Crassidis/x.csv', delimiter=',')[:-1]
zs_tmp = np.genfromtxt('c:/users/rlabbe/dropbox/Crassidis/mycode/xx.csv', delimiter=',')[:-1]
zs = zs_tmp
except:
pass
Expand Down Expand Up @@ -232,8 +229,7 @@ def test_imm():
trans = np.array([[0.97, 0.03],
[0.03, 0.97]])


bank = IMMEstimator((6, 1), filters, (0.5, 0.5), trans)
bank = IMMEstimator(filters, (0.5, 0.5), trans)

xs, probs = [], []
cvxs, caxs = [], []
Expand All @@ -248,7 +244,7 @@ def test_imm():
#print(i, ca.likelihood, cano.likelihood, bank.w)

#print('p', bank.p)
probs.append(bank.w.copy())
probs.append(bank.mu.copy())

if DO_PLOT:
xs = np.array(xs)
Expand Down

0 comments on commit b63f69b

Please sign in to comment.