Daniel@0: % Make a point move in the 2D plane
Daniel@0: % State = (x y xdot ydot). We only observe (x y).
Daniel@0: 
Daniel@0: % This code was used to generate Figure 15.9 of "Artificial Intelligence: a Modern Approach",
Daniel@0: % Russell and Norvig, 2nd edition, Prentice Hall, 2003.
Daniel@0: 
Daniel@0: % X(t+1) = F X(t) + noise(Q)
Daniel@0: % Y(t) = H X(t) + noise(R)
Daniel@0: 
Daniel@0: ss = 4; % state size
Daniel@0: os = 2; % observation size
Daniel@0: F = [1 0 1 0; 0 1 0 1; 0 0 1 0; 0 0 0 1]; 
Daniel@0: H = [1 0 0 0; 0 1 0 0];
Daniel@0: Q = 0.1*eye(ss);
Daniel@0: R = 1*eye(os);
Daniel@0: initx = [10 10 1 0]';
Daniel@0: initV = 10*eye(ss);
Daniel@0: 
Daniel@0: seed = 9;
Daniel@0: rand('state', seed);
Daniel@0: randn('state', seed);
Daniel@0: T = 15;
Daniel@0: [x,y] = sample_lds(F, H, Q, R, initx, T);
Daniel@0: 
Daniel@0: [xfilt, Vfilt, VVfilt, loglik] = kalman_filter(y, F, H, Q, R, initx, initV);
Daniel@0: [xsmooth, Vsmooth] = kalman_smoother(y, F, H, Q, R, initx, initV);
Daniel@0: 
Daniel@0: dfilt = x([1 2],:) - xfilt([1 2],:);
Daniel@0: mse_filt = sqrt(sum(sum(dfilt.^2)))
Daniel@0: 
Daniel@0: dsmooth = x([1 2],:) - xsmooth([1 2],:);
Daniel@0: mse_smooth = sqrt(sum(sum(dsmooth.^2)))
Daniel@0: 
Daniel@0: 
Daniel@0: figure(1)
Daniel@0: clf
Daniel@0: %subplot(2,1,1)
Daniel@0: hold on
Daniel@0: plot(x(1,:), x(2,:), 'ks-');
Daniel@0: plot(y(1,:), y(2,:), 'g*');
Daniel@0: plot(xfilt(1,:), xfilt(2,:), 'rx:');
Daniel@0: for t=1:T, plotgauss2d(xfilt(1:2,t), Vfilt(1:2, 1:2, t)); end
Daniel@0: hold off
Daniel@0: legend('true', 'observed', 'filtered', 3)
Daniel@0: xlabel('x')
Daniel@0: ylabel('y')
Daniel@0: 
Daniel@0: 
Daniel@0: 
Daniel@0: % 3x3 inches
Daniel@0: set(gcf,'units','inches');
Daniel@0: set(gcf,'PaperPosition',[0 0 3 3])  
Daniel@0: %print(gcf,'-depsc','/home/eecs/murphyk/public_html/Bayes/Figures/aima_filtered.eps');
Daniel@0: %print(gcf,'-djpeg','-r100', '/home/eecs/murphyk/public_html/Bayes/Figures/aima_filtered.jpg');
Daniel@0: 
Daniel@0: 
Daniel@0: figure(2)
Daniel@0: %subplot(2,1,2)
Daniel@0: hold on
Daniel@0: plot(x(1,:), x(2,:), 'ks-');
Daniel@0: plot(y(1,:), y(2,:), 'g*');
Daniel@0: plot(xsmooth(1,:), xsmooth(2,:), 'rx:');
Daniel@0: for t=1:T, plotgauss2d(xsmooth(1:2,t), Vsmooth(1:2, 1:2, t)); end
Daniel@0: hold off
Daniel@0: legend('true', 'observed', 'smoothed', 3)
Daniel@0: xlabel('x')
Daniel@0: ylabel('y')
Daniel@0: 
Daniel@0: 
Daniel@0: % 3x3 inches
Daniel@0: set(gcf,'units','inches');
Daniel@0: set(gcf,'PaperPosition',[0 0 3 3])  
Daniel@0: %print(gcf,'-djpeg','-r100', '/home/eecs/murphyk/public_html/Bayes/Figures/aima_smoothed.jpg');
Daniel@0: %print(gcf,'-depsc','/home/eecs/murphyk/public_html/Bayes/Figures/aima_smoothed.eps');