% Example2_30.m Sequential solution of Example 2.30
% Estimate polynomial for 
clear all; format short e
opt = 1;   % can be opt = 2 for 3 parameter model
if( opt == 1 )
    p=5;  mu=[300 10 1 .02 .0004]; 
else
    p=3;  mu=[300 10 1 ]; 
end
n=30; t=0:1:n-1; u=[t]';
y=[279.59 264.87 251.53 239.3 228.18 217.24 207.86 199.36 191.65 184.44 ...
177.64 171.41 165.04 159.89 155.19 150.78 146.68 ...
142.93 139.34 136.04 132.94 130.07 127.39 124.88 122.46 120.18 118.04 ...
115.97 114.13 112.35]'; 
if( opt == 1 )
    X=[ones(n,1),u,u.^2,u.^3,u.^4];     
else
    X=[ones(n,1),u,u.^2 ];     
end
sigma=1; sig=sigma*ones(n,1);
b=zeros(p,1); P0=1e+5*diag(mu.*mu); P=P0; B=b'; %P0 is V sub(beta)
for ii=1:n			%SEQUENTIAL
    A=P*X(ii,:)';	%eq. (2.9.8a)	%SEQUENTIAL
    Delta=sig(ii)^2+X(ii,:)*A;	%eq. (2.9.8)	%SEQUENTIAL
    K=A/Delta;	%eq. (2.9.8c)	%SEQUENTIAL
    e=y(ii)-X(ii,:)*b;	%eq. (2.9.8d)	%SEQUENTIAL
    b=b+K*e;	%eq. (2.9.8e)	%SEQUENTIAL
    P=P-K*A';	%eq. (2.9.8f)	%SEQUENTIAL
   B=[B; b']; 
end
B
figure(1); clf;
nos = [0:30]';
plot( nos, B(:,1)/B(end,1),nos, B(:,2)/B(end,2),nos, B(:,3)/B(end,3),...
      'linewidth', 1 )
xlabel('Measurement number')
ylabel('Estimated parameter/(final value)')
if( opt == 1)
    text(8.75,1.5,'parameter 3'); text(7,0.9,'1'); text(7,1.1,'2');
    text(0.5,1.05,'1'); text(0.8,0.5,'2'); text(1.6,0,'3')
    fig=figure(1);
    print(fig,fullfile(pwd, 'Figure2_13a.tif'),'-dtiff','-r600')
    figure(2); clf;
    plot( nos, B(:,4)/B(end,4),nos, B(:,5)/B(end,5), ...
          'linewidth', 1 )
    xlabel('Measurement number')
    ylabel('Estimated parameter/(final value)')
    text(9,42,'parameter 5')
    text(8,10.3,'4')
    fig=figure(2);
    print(fig,fullfile(pwd, 'Figure2_13b.tif'),'-dtiff','-r600')
else
    text(7,2.6,'parameter 3'); text(5,1.5,'parameter 2'); 
    text(5,0.9,'parameter 1');
    if( opt == 2 )
        fig=figure(1);
        print(fig,fullfile(pwd, 'Figure2_14.tif'),'-dtiff','-r600')
    end
    return
end

