Commit 05629b13 authored by Saad Jbabdi's avatar Saad Jbabdi
Browse files

Merge branch 'add_kalman' into 'master'

Add kalman

See merge request !8
parents ef5e8dd8 45e23f58
# Kalman Filter practical
This is a very short Matlab practical that implements a 1D, linear Kalman Filter.
Please refer to the presentation slides for the derivation of the update equations.
\ No newline at end of file
%% This is a very short example of a very simple implementation of the Kalman filter,
% which I presented in the Graphical Models talk.
% Please refer to the slides to understand the derivation of the update equations
% The Graphical model: (xk is the state parameters and yk are the oservations)
%
% x0->x1 -> ... -> xk-1 -> xk -> ...
% | | |
% V V V
% y1 yk-1 yk
%
%
% We will use the Oxford Weather data as a dataset to play with.
%
data = load('weather_data.txt');
year = data(:,1);
month = data(:,2);
maxTemp = data(:,3);
minTemp = data(:,4);
meanTemp = (maxTemp+minTemp)/2;
hoursFrost = data(:,5);
rain = data(:,6);
hoursSun = data(:,7);
% If you would like to visualise these data, you can use the code below:
figure
subplot(2,3,1),boxplot(rain,round(year/10)*10); xlabel('decade'); ylabel('rain');
subplot(2,3,2),boxplot(rain,month); xlabel('month'); ylabel('rain');
subplot(2,3,3),plot(minTemp,maxTemp,'.k'); xlabel('minTemp'); ylabel('maxTemp');
subplot(2,3,4),plot(meanTemp,hoursFrost,'.k'); xlabel('meanTemp'); ylabel('hoursFrost');
subplot(2,3,5),boxplot(rain,round(meanTemp/3)*3); xlabel('meanTemp'); ylabel('rain');
subplot(2,3,6),boxplot(hoursSun,month); xlabel('month'); ylabel('housrSun');
% choose your favourite data among the above
y = maxTemp(1:200);
n = length(y);
% set up noise precisions
% here we pretend that we know the values for these parameters
% however, there are variations on the algorithm where we can estimate these
% Feel free to change the values of these two and observe the effect this has on the
% model fit
b=1; % observation noise
b0=1; % state noise
% first state is first data point
x0=y(1);
% The forward model is very simple:
% yk ~ N(xk,1/b)
% xk ~ N(xk-1,1/b0)
% or equivalently
% yk = xk + noise(0,1/b) --> observation equation
% xk = xk-1 + noise(0,1/b0) --> state dynamics equation
% In general, both x and y can be multi-dimensional vectors, and the
% state/observation equations can involve matrix multiplications.
% To keep the equations simple, we are here using 1-dimensional
% vectors and the equations do not involve any extra constants.
% initialise the vector of states
xk=zeros(n,1);
bk=zeros(n,1);
% forward updates
xk(1) = (b*y(1)+b0*x0)/(b+b0);
bk(1) = b+b0;
for k=2:n
% update for the precision
bk(k) = b+b0*bk(k-1)/(b0+bk(k-1));
% update for the state parameter
% this is the key equation. it says that the new
% state param xk is equal to the old xk-1 plus
% a correction term that involves the new observation yk,
% weighted by the current estimate for the precision
xk(k) = xk(k-1) + b/bk(k)*(y(k)-xk(k-1));
end
% Now visualise the data, the state parameters (which in this model are simply tring to track the data).
% We also have an estimate of the state precision (bk), which we use to show +/- one standard deviation
figure
plot(y,'.-');
hold on
plot(xk,'r-');
plot(xk+sqrt(1./bk),'r--');
plot(xk-sqrt(1./bk),'r--');
% That's All for now. I told you this one was short!
This diff is collapsed.
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment