% Example2_28b.m
% Sequential analysis of example 2.13, line source problem.
clear all; format short g; 
% data below from Example 2.13
t=[0.7916; 1.2916; 2.0833; 3.4582; 5.7081; 9.2913; 15.1660; 24.7490; 40.3320 ];
y =[65.2010; 66.2620; 67.3580; 68.5790; 69.8210; 71.0220; 72.0390; 73.2130; 74.3840 ];
C =[0.0200; 0.0200; 0.0200; 0.0200; 0.0200; 0.0200; 0.0800; 0.0800; 0.0800 ];
n= 9; p=2;

% use the line below for logarithmic time variation
X=[ones(n,1) (0.25*36.3/pi)*log(t)]; 

% alternatively, use the line below for exponential integral
% X=[ones(n,1) (0.25*36.3/pi)*expint((0.25/3.2)./t)];

Cr=C.*C; 
mu=zeros(p,1); %mu is the prior parameter estimate vector
sig=C; Pz=1.0e+5*sig(1)^2; b=mu;
P0=Pz*[ 4000 0; 0 1]; P=P0; 
B=[mu' P(1,1) P(1,2) P(2,2)]; %P0 is V sub(beta)
for ii=1:n			%SEQUENTIAL
% these lines from Example 2.9.1 code
    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
% end lines from Example 2.9.1
   B=[B; b' P(1,1) P(1,2) P(2,2)]; 
end
B(:,3)=B(:,3)*1000; B(:,4)=B(:,4)*1000; B(:,5)=B(:,5)*1000;
fprintf(' i     T0         K=1/k     P(1,1)     P(1,2)       P(2,2)\n')
fprintf(' %1d   %6.2f    %8.4f %10.4g  %10.4g  %10.4g\n',[ [0:n]' B ]')
fprintf('Final values:  T0=%6.3f F,           s.e.T0=%6.4f F\n',b(1), sqrt(P(1,1)) )
fprintf('               K =%7.5f hr-ft-F/Btu, s.e.k=%9.6f hr-ft-F/Btu\n', ...
             b(2), sqrt(P(2,2))  )
fprintf('               Standard error in temperature: %8.5f\n', ...
             sqrt( [y-X*b]'*[y-X*b]/length(y) ) );
figure(1); clf;
plot( [1:n], B(2:end,1), 'ro', [1:n], 50*B(2:end,2), 'bx', 'linewidth', 1 )
legend('T0', '50 times K', 'location', 'southeast')
ylabel('parameters T0 and 50*K')
xlabel('sequential step')