kalman filter - "Machine Learning An Algorithmic Perspective"

 

ch16. Graphical Models

• A simple 1D Kalman Filter

# Code from Chapter 16 of Machine Learning: An Algorithmic Perspective (2nd Edition)
# by Stephen Marsland (http://stephenmonika.net)

# You are free to use, change, or redistribute the code in any way you wish for
# non-commercial purposes, but please maintain the name of the original author.
# This code comes with no warranty of any kind.

# Stephen Marsland, 2008, 2014

# The 1D Kalman filter

import pylab as pl
import numpy as np

def Kalman(obs=None,mu_init=np.array([-0.37727]),cov_init=0.1*np.ones((1)),nsteps=50):

    ndim = np.shape(mu_init)[0]
    
    if obs==None:
        mu_init = np.tile(mu_init,(1,nsteps))
        cov_init = np.tile(cov_init,(1,nsteps))
        obs = np.random.normal(mu_init,cov_init,(ndim,nsteps))
    
    Sigma_x = np.eye(ndim)*1e-5
    A = np.eye(ndim)
    H = np.eye(ndim)
    mu_hat = 0
    cov = np.eye(ndim)
    R = np.eye(ndim)*0.01
    
    m = np.zeros((ndim,nsteps),dtype=float)
    ce = np.zeros((ndim,nsteps),dtype=float)
    
    for t in range(1,nsteps):
        # Make prediction
        mu_hat_est = np.dot(A,mu_hat)
        cov_est = np.dot(A,np.dot(cov,np.transpose(A))) + Sigma_x

        # Update estimate
        error_mu = obs[:,t] - np.dot(H,mu_hat_est)
        error_cov = np.dot(H,np.dot(cov,np.transpose(H))) + R
        K = np.dot(np.dot(cov_est,np.transpose(H)),np.linalg.inv(error_cov))
        mu_hat = mu_hat_est + np.dot(K,error_mu)
        #m[:,:,t] = mu_hat
        m[:,t] = mu_hat
        if ndim>1:
            cov = np.dot((np.eye(ndim) - np.dot(K,H)),cov_est)
        else:
            cov = (1-K)*cov_est 
        ce[:,t] = cov                                
    
    pl.figure()
    pl.plot(obs[0,:],'ko',ms=6)
    pl.plot(m[0,:],'k-',lw=3)
    pl.plot(m[0,:]+20*ce[0,:],'k--',lw=2)
    pl.plot(m[0,:]-20*ce[0,:],'k--',lw=2)
    pl.legend(['Noisy Datapoints','Kalman estimate','Covariance'])
    pl.xlabel('Time')
    
    
    pl.show()
    
Kalman()

 

• A complete Kalman Filter

# Code from Chapter 16 of Machine Learning: An Algorithmic Perspective (2nd Edition)
# by Stephen Marsland (http://stephenmonika.net)

# You are free to use, change, or redistribute the code in any way you wish for
# non-commercial purposes, but please maintain the name of the original author.
# This code comes with no warranty of any kind.

# Stephen Marsland, 2008, 2014

import numpy as np
import pylab as pl

def Kalman_update(A,H,Q,R,y,x,Sig,B=None,u=None):

	if B is None:
		xpred = np.dot(A,x) 
	else:
		xpred = np.dot(A,x) + np.dot(B,u)

	SigPred = np.dot(A,np.dot(Sig,A.T)) + Q

	e = y - np.dot(H,xpred)
	Sinv = np.linalg.inv(np.dot(H,np.dot(SigPred,H.T)) + R)
	K = np.dot(SigPred,np.dot(H.T,Sinv))

	xnew = xpred + np.dot(K,e)
	SigNew = np.dot((np.eye(np.shape(A)[0]) - np.dot(K,H)),SigPred)

	return xnew.T,SigNew
	
def Kalman_smoother_update(A,Q,B,u,xs_t,Sigs_t,xfilt,Sigfilt,Sigfilt_t):

	if B is None:
		xpred = np.dot(A,xfilt)
	else:
		xpred = np.dot(A,xfilt) + np.dot(B,u)

	SigPred = np.dot(A,np.dot(Sigfilt,A.T)) + Q
	J = np.dot(Sigfilt,np.dot(A.T,np.linalg.inv(SigPred)))
	xs = xfilt + np.dot(J,(xs_t - xpred))
	Sigs = Sigfilt + np.dot(J,np.dot((Sigs_t - SigPred),J.T))

	return xs.T, Sigs


def Kalman_filter(y,A,H,Q,R,x0,Sig0,B=None,u=None):

	obs_size,T = np.shape(y)
	state_size = np.shape(A)[0]

	x = np.zeros((state_size,T))
	Sig = np.zeros((state_size,state_size,T))

	[x[:,0],Sig[:,:,0]] = Kalman_update(A,H,Q,R,y[:,0].reshape(len(y),1),x0,Sig0,B,u)
	for t in range(1,T):
		prevx = x[:,t-1].reshape(state_size,1)
		prevSig = Sig[:,:,t-1]
		[x[:,t],Sig[:,:,t]] = Kalman_update(A,H,Q,R,y[:,t].reshape(len(y),1),prevx,prevSig,B,u)

	return x,Sig

def Kalman_smoother(y,A,H,Q,R,x0,Sig0,B=None,u=None):

	obs_size,T = np.shape(y)
	state_size = np.shape(A)[0]

	xs = np.zeros((state_size,T))
	Sigs = np.zeros((state_size,state_size,T))

	[xfilt,Sigfilt] = Kalman_filter(y,A,H,Q,R,x0,Sig0,B,u)

	xs[:,T-1] = xfilt[:,T-1]
	Sigs[:,:,T-1] = Sigfilt[:,:,T-1]

	for t in range(T-2,-1,-1):
		[xs[:,t],Sigs[:,:,t]] = Kalman_smoother_update(A,Q,B,u,xs[:,t+1].reshape(len(xs),1),Sigs[:,:,t+1],xfilt[:,t].reshape(len(xfilt),1),Sigfilt[:,:,t],Sigfilt[:,:,t+1])

	return xs,Sigs

def lds_sample(A,H,Q,R,state0,T):
	# x(t+1) = Ax(t) +  state_noise(t), state_noise ~ N(O,Q), x(0) = state0
	# y(t) = Hx(t) + obs_noise(t), obs_noise~N(O,R)

	state_noise_samples = np.random.multivariate_normal(np.zeros((len(Q))),Q,T).T
	obs_noise_samples = np.random.multivariate_normal(np.zeros((len(R))),R,T).T

	x = np.zeros((np.shape(H)[1],T))
	y = np.zeros((np.shape(H)[0],T))

	x[:,0] = state0.T
	y[:,0] = np.dot(H,x[:,0]) + obs_noise_samples[:,0]

	for t in range(1,T):
		x[:,t] = np.dot(A,x[:,t-1]) + state_noise_samples[:,t]
		y[:,t] = np.dot(H,x[:,t-1]) + obs_noise_samples[:,t]

	return [x,y]

def Kalman_demo():
	state_size = 4
	observation_size = 2
	A = np.array([[1,0,1,0],[0,1,0,1],[0,0,1,0],[0,0,0,1]],dtype=float)
	H = np.array([[1,0,0,0],[0,1,0,0]],dtype=float)

	Q = 0.1*np.eye((state_size))
	R = np.eye(observation_size,dtype=float)

	x0 = np.array([[10],[10],[1],[0]],dtype=float)
	Sig0 = 10. * np.eye(state_size)

	np.random.seed(3)
	T = 15
	
	[x,y] = lds_sample(A,H,Q,R,x0,T)

	[xfilt,Sigfilt] = Kalman_filter(y,A,H,Q,R,x0,Sig0)
	[xsmooth,Sigsmooth] = Kalman_smoother(y,A,H,Q,R,x0,Sig0)

	dfilt = x[[0,1],:] - xfilt[[0,1],:]
	mse_filt = np.sqrt(np.sum(dfilt**2))

	dsmooth = x[[0,1],:] - xsmooth[[0,1],:]
	mse_smooth = np.sqrt(np.sum(dsmooth**2))

	plot_track(x,y,xfilt,Sigfilt)
	plot_track(x,y,xsmooth,Sigsmooth)
	
def plot_track(x,y,Kx,Sig):
	fig = pl.figure()
	ax = fig.add_subplot(111, aspect='equal')
	pl.plot(x[0,:],x[1,:],'ks-')
	pl.plot(y[0,:],y[1,:],'k*')
	pl.plot(Kx[0,:],Kx[1,:],'kx:')
	pl.legend(('True','Observed','Filtered'))

	obs_size,T = np.shape(y)

	from matplotlib.patches import Ellipse
	# Axes of ellipse are eigenvectors of covariance matrix, lengths are square roots of eigenvalues
	ellsize = np.zeros((obs_size,T))
	ellangle = np.zeros((T))
	for t in range(T):
		[evals,evecs] = np.linalg.eig(Sig[:2,:2,t])
		ellsize[:,t] = np.sqrt(evals)	
		ellangle[t] = np.angle(evecs[0,0]+0.j*evecs[0,1])
		
	ells = [Ellipse(xy=[Kx[0,t],Kx[1,t]] ,width=ellsize[0,t],height=ellsize[1,t], angle=ellangle[t]) for t in range(T)]
	for e in ells:
		ax.add_artist(e)
		e.set_alpha(0.1)
		e.set_facecolor([0.7,0.7,0.7])
	pl.xlabel('x')
	pl.ylabel('y')


def Kalman_demo1d():

	x0 = np.array([-0.37727])
	Sig0 = 0.1*np.ones((1))
	T = 50	

        y = np.random.normal(x0,Sig0,(1,T))

    	A = np.eye(1)
    	H = np.eye(1)
    	Q = np.eye(1)*1e-5
    	R = np.eye(1)*0.01
    
    	xfilt = np.zeros((1,T),dtype=float)
    	Sigfilt = np.zeros((1,T),dtype=float)

	[xfilt,Sigfilt] = Kalman_filter(y,A,H,Q,R,x0,Sig0)
	xfilt = np.squeeze(xfilt)
	Sigfilt = np.squeeze(Sigfilt)

    	pl.figure()
	time = np.arange(T)
    	pl.plot(time,y[0,:],'ko',ms=6)
    	pl.plot(time,xfilt,'k-',lw=3)
    	pl.plot(time,xfilt+20*Sigfilt,'k--',lw=2)
    	pl.plot(time,xfilt-20*Sigfilt,'k--',lw=2)
    	pl.legend(['Noisy Datapoints','Kalman estimate','20*Covariance'])
    	pl.xlabel('Time')

 

 

• The Extended Kalman Filter

# Code from Chapter 16 of Machine Learning: An Algorithmic Perspective (2nd Edition)
# by Stephen Marsland (http://stephenmonika.net)

# You are free to use, change, or redistribute the code in any way you wish for
# non-commercial purposes, but please maintain the name of the original author.
# This code comes with no warranty of any kind.

# Stephen Marsland, 2008, 2014

import numpy as np
import pylab as pl

def EKF_update(Q,R,y,x,t,Sig,B=None,u=None):

	if B is None:
		xpred = f(x,t).reshape(len(x),1)
	else:
		xpred = f(x,t).reshape(len(x),1) + np.dot(B,u)

	A = np.array(jac_f(x,t))
	#A = jac_f(x,t).reshape(1,len(x))
	SigPred = np.dot(A,np.dot(Sig,A.T)) + Q

	H = jac_h(xpred).reshape(1,len(x))
	e = y - h(xpred)
	Sinv = 1./(np.dot(H,np.dot(SigPred,H.T)) + R)
	#Sinv = np.linalg.inv(np.dot(H,np.dot(SigPred,H.T)) + R)
	K = np.dot(SigPred,np.dot(H.T,Sinv))

	xnew = xpred + np.dot(K,e)
	SigNew = np.dot((1 - np.dot(K,H)),SigPred)

	return xnew.T,SigNew
	
def EKF(y,Q,R,x0,Sig0,B=None,u=None):

	obs_size,T = np.shape(y)
	state_size = np.shape(Q)[0]

	x = np.zeros((state_size,T))
	Sig = np.zeros((state_size,state_size,T))

	[x[:,0],Sig[:,:,0]] = EKF_update(Q,R,y[:,0].reshape(len(y),1),x0,0,Sig0,B,u)
	for t in range(1,T):
		prevx = x[:,t-1].reshape(state_size,1)
		prevSig = Sig[:,:,t-1]
		[x[:,t],Sig[:,:,t]] = EKF_update(Q,R,y[:,t].reshape(len(y),1),prevx,t,prevSig,B,u)

	return x,Sig

def f(x,t):
	return np.array([x[1],x[2],0.5*x[0]*(x[1]+x[2])]).T

def jac_f(x,t):
	return np.array([[0,0,0.5*(x[1]+x[2])],[1,0,0.5*x[0]],[0,1,0.5*x[0]]])
	#return x

def jac_h(x):
	#return np.array([x[2]*np.sin(x[0]), 0, np.cos(x[0])])
	return np.array([1.0,1.0,0.0]) #np.array([0.4*x]).reshape(1,1)
	#return np.array([1.0,0.0,0.0]) #np.array([0.4*x]).reshape(1,1)

def h(x):
	return x[0]+x[1]

def EKF_demo():

	state_size = 3
	observation_size = 1

	Q = 0.1*np.eye(state_size)
	R = 0.1*np.eye(observation_size)
	x0 = np.array([0,0,1])
	Sig0 = np.eye(state_size)
	
	T = 20
	state = np.zeros((state_size,T))
	y = np.zeros((observation_size,T))
	state[:,0] = x0.T

	for t in range(1,T):
		state[:,t] = f(state[:,t-1],t) + np.random.multivariate_normal(np.zeros((len(Q))),Q)
		y[:,t] = h(state[:,t]) + np.sqrt(R)*np.random.randn()
		#state[:,t] = np.dot(A,state[:,t-1]) + np.random.multivariate_normal(np.zeros((len(Q))),Q)
		#y[:,t] = h(state[:,t]) + np.random.randn()

	[xfilt,Sigfilt] = EKF(y,Q,R,x0,Sig0)

	dfilt = state[0,:] - xfilt[0,:]
	#dfilt = state[[0,1],:] - xfilt[[0,1],:]
	mse_filt = np.sqrt(np.sum(dfilt**2))
	print mse_filt

	ypred = np.zeros((T))
	for t in range(T):
		ypred[t] = h(xfilt[:,t])

	print Sigfilt
	pl.figure()
	pl.plot(np.arange(T),np.squeeze(y),'*')
	pl.plot(np.arange(T),ypred,'k-')
	

https://homepages.ecs.vuw.ac.nz/~marslast/MLbook.html

 

 

 

Stephen Marsland

 

 

 

This webpage contains the code and other supporting material for the textbook "Machine Learning: An Algorithmic Perspective" by Stephen Marsland, published by CRC Press, part of the Taylor and Francis group. The first edition was published in 2009, and a revised and updated second edition is due out towards the end of 2014.

The book is aimed at computer science and engineering undergraduates studing machine learning and artificial intelligence.

The table of contents for the second edition can be found here.

There are lots of Python/NumPy code examples in the book, and the code is available here. Datasets (either the actual data, or links to the appropriate resources) are given at the bottom of the page.

Note that the chapter headings and order below refer to the second edition. However, the titles of the chapters should enable users of the first edition to find the relevant sections. In addition, a zip file of the code for the 1st edition is available here.

All of the code is freely available to use (with appropriate attribution), but comes with no warranty of any kind.

Option 1: Zip file of all code, arranged into chapters
Option 2: Choose what you want from here:

• Chapter 2 (Preliminaries):
  • Plots a 1D Gaussian (Fig 2.14)
  • Plot some 2D Gaussians
• Chapter 3 (Neurons, Neural Networks, and Linear Discriminants):
  • The Perceptron
  • The Linear Regressor
  • Another Perceptron (for use with logic.py)
  • Demonstration of Perceptron with logic functions
  • Demonstration of Linear Regressor with logic functions
  • Demonstration of Perceptron with Pima Indian dataset
  • Demonstration of Linear Regressor with auto-mpg dataset
  • Demonstration of Perceptron with the MNIST dataset
• Chapter 4 (The Multi-Layer Perceptron):
  • The Multi-Layer Perceptron
  • Demonstration of the MLP on logic functions
  • Demonstration of the MLP for classification on the Iris dataset
  • Demonstration of the MLP for regression on data from a sine wave
  • Demonstration of the MLP for time series on the Palmerston North Ozone dataset
  • Demonstration of MLP with the MNIST dataset
• Chapter 5 (Radial Basis Functions and Splines):
  • The Radial Basis Function
  • Linear Least Squares Fitting
  • Demonstration of the RBF on the Iris dataset
• Chapter 6 (Dimensionality Reduction):
  • Linear Discriminant Analysis
  • Principal Components Analysis
  • Factor Analysis
  • Kernel Principal Components Analysis
  • Locally Linear Embedding
  • Isomap
  • Demonstration of PCA
  • Demonstration of kernel PCA on a set of circular data
  • Demonstration of the algorithms on the Iris dataset
• Chapter 7 (Probabilistic Learning):
  • The Gaussian Mixture Model
  • The k-Nearest Neighbour Algorithm
  • The k-Nearest Neighbour Smoother
  • The kd-Tree Algorithm
• Chapter 8 (Support Vector Machines):
  • The Support Vector Machine (needs cvxopt)
  • Demonstration of the SVM for classification on the Iris dataset
  • Demonstration of the SVM for the variant of XOR in Figs 8.7 and 8.8
• Chapter 9 (Optimisation and Search):
  • Steepest Descent
  • Newton's method
  • Levenberg-Marquarft
  • Conjugate Gradients
  • The version of the MLP algorithm trained using conjugate gradients
  • Demonstration of the MLP algorithm trained using conjugate gradients on the Iris dataset
  • Demonstration of Levenberg-Marquardt on a least-squares fitting problem
  • Demonstration of four solution methods for the Travelling Salesman Problem
• Chapter 10 (Evolutionary Learning):
  • The Genetic Algorithm
  • A Runner for the GA
  • Population-Based Incremental Learning
  • A knapsack problem fitness function
  • The four peaks fitness function
  • The onemax fitness function
  • Exhaustive search algorithm to solve the knapsack problem
  • A greedy algorithm to solve the knapsack problem
• Chapter 11 (Reinforcement Learning):
  • The SARSA algorithm
  • The TD(0) algorithm
  • Demonstration of the SARSA algorithm on the Cliff problem
  • Demonstration of the TD(0) algorithm on the Cliff problem
• Chapter 12 (Learning with Trees):
  • The decision tree
  • Demonstration of the decision tree on the Party dataset
  • The Party dataset
• Chapter 13 (Decision by Committee: Ensemble Learning):
  • The Boosting algorithm
  • The Bagging algorithm
  • A decision tree with weights
  • The random forest algorithm
  • Demonstration of bagging, stumping and random forests on the Party dataset
  • Demonstration of bagging on the Car Safety dataset
  • Demonstration of bagging and random forests on the mushroom dataset
• Chapter 14 (Unsupervised Learning):
  • The k-Means Algorithm
  • The k-Means Neural Network
  • The Self-Organising Map
  • A Simple 2D Demonstration of the SOM
  • Demonstration of k-Means and the SOM on the Iris dataset
  • Demonstration of SOM and Perceptron together on the Iris dataset
  • More demonstrations of the SOM
• Chapter 15 (Markov Chain Monte Carlo Methods):
  • The Linear Congruential Pseudo-Random Number Generator
  • The Box-Muller method of constructing Gaussian-distributed pseudo-random numbers
  • The Rejection Sampling Algorithm
  • The Importance Sampling Algorithm
  • The Sampling-Importance-Resampling Algorithm
  • The Metropolis-Hastings Algorithm
  • The Gibbs Sampler
• Chapter 16 (Graphical Models):
  • The Gibbs Sampler for the Exam Panic dataset
  • The Hidden Markov Model
  • A simple 1D Kalman Filter
  • A complete Kalman Filter
  • The Extended Kalman Filter
  • The Basic Particle Filter
  • A Tracking Particle Filter
  • The Markov Random Field for Image Denoising
  • A demonstration of finding paths in graphs
  • An image for denoising
• Chapter 17 (Symmetric Weights and Deep Belief Networks):
  • A Hopfield network
  • The Restricted Boltzmann Machine
  • The Deep Belief Network Algorithm
• Chapter 18 (Gaussian Processes):
  • The Gaussian Process for Regression Algorithm
  • The Gaussian Process for Classification Algorithm
  • Demo of the Gaussian Process for Classification
  • Plots of the Weibull and Gaussian distributions (Fig 18.1)
  • Plots of GP samples (Fig 18.2)
  • Very simple data for the Gaussian Process Regression demo

Datasets

Many of the datasets used in the book are available from the UCI Machine Learning Repository. In particular, look for the Iris data, the Pima Indian data, the car safety data, the auto-mpg data, the wine data, and the mushroom data.

Two of the most popular machine learning demonstration datasets are the MNIST set of zip code digits, which is available here, and the binary alpha digits dataset, which can be downloaded here. Finally, there are a couple of smaller datasets that are not available elsewhere, at least in their current form, and so should be downloaded from this website:

• The Palmerston North Ozone dataset
• Training data for the prostate dataset
• Test data for the prostate dataset (variables are log cancer volume, log prostate weight, age, lbph, svi, lcp, Gleason score, pgg45 and the last one is response lpsa)
• The Ruapehu dataset (thanks to Mark Bebbington)
• Short version of the e. coli dataset