%Example2_31.m Sequential MAP estimation Example 2.31
clear all; format short e
opt = 1;   % can be opt = 2 for no prior information (OLS)
if( opt == 1 )
    p=2;  mu=[11 14]';
    P0 = [ 0.04 -0.01; -0.01 0.04];
else
    p=2;  mu=[0  0]'; 
    P0 = 1000*eye(2,2);
end
T = [ 20 20 21 100 150 200 250 297 300 302]';
y = [11.47 10.94 11.15 11.85 12.55 13.18 13.48 13.90 14.54 14.36]';
sigma = 0.2; Ta=20; Tb=300;
n=length(T); t=0:1:n-1; u=[t]';
X=[ (Tb-T)/(Tb-Ta)  (T-Ta)/(Tb-Ta) ];     
sig=sigma*ones(n,1);
b=mu; P=P0; B= [b' P(1,1) P(1,2) P(2,2)]; %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' P(1,1) P(1,2) P(2,2)]; 
end
fprintf('%3d  %6.3f  %6.3f  %9.7f %10.7f  %9.7f\n', [ [0:n]' B]')
% 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(7,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')
% else
%     text(7,2.6,'parameter 3'); text(5,1.5,'parameter 2'); 
%     text(5,0.9,'parameter 1');
%     return
% end
% 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(7.8,42,'parameter 5')
% text(7,10,'4')
