added kalman

markov_models/Kalman/README.md

+# 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

markov_models/Kalman/kalman_filter.m

+%% 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);
+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!
+
+