Sequential state estimation of a simple system x_(n+1) = ax_n

a is assumed known

Contents

First set parameters and generate truth

a = 1.1;
sigo = 2; % std of obs about truth
N = 20; % Number of samples
x0 = 1; % Generates truth

xtrue = x0*a.^(1:N);

% Now make observations and plot on top of truth

xo = xtrue + sigo*randn(1,N);
clf
plot(1:N,xtrue,'k+',1:N,xo,'rx')
xlim([0.5 N+0.5])
xlabel('n')

Sequential state estimation

xhat = nan(1,N);
sighat2 = nan(1,N);
xhat(1) = xo(1);
sighat2(1)= sigo^2;
for n = 2:N
  xp(n) = a*xhat(n-1);
  sigp2(n) = a^2*sighat2(n-1);
  sighat2(n) = 1/(1/(sigo^2) + 1/sigp2(n));
  c = sighat2(n)/sigo^2;
  xhat(n) = xp(n) + c*(xo(n)-xp(n));
end

hold on
plot(1:N,xhat,'bx',1:N,xhat+2*sqrt(sighat2),'b^',1:N,xhat-2*sqrt(sighat2),'bv')
legend('truth','obs','est','+2\sigma','-2\sigma','Location','SouthEast')
hold off

Published with MATLAB® R2012b