% Example2_24.m  Sequential experiments
% Each experiment temperatures in columns of t
t = [ [1:0.1:2]' [5:0.1:6]' [10:0.1:11]' ];
s = [ 0.01 0.02 0.03 ];   % standard deviations for the three experiments
beta = [ 1; 0.1; -0.005 ];  % exact model parameters

close(1)
figure(1); clf; 
hold on

m = 3;             % number of experiments
n = length( t );   % number of data points
p = length(beta);  % number of parameters

bavg = zeros(3,1);

rng(129);   % always use the same random numbers
% Analysis 1
fprintf('Analysis 1\n');
for im = 1:m
    omega = diag( ones(1,n)/s(im)^2 );  
    X = [ ones(n,1) t(:,im) t(:,im).^2 ];
    y = X*beta + randn(n,1)*s(im);
    ydata(:,im) = y;  % save the data 
    plot( t(:,im), y, 'ro', 'linewidth', 1.0 )
    b = X \ y;
%    lscov( X, y, omega )    % this could also be used
    fprintf('Experiment %2d:\n',im)
    fprintf('b = %10.4g  %10.4g  %10.4g\n',b)
    fprintf('cond(X^T*omega*X) = %10.6g\n', cond(X'*omega*X) )
    bavg = bavg + b/3;
end
fprintf('bavg = %10.4g  %10.4g  %10.4g\n',bavg)
hold off
xlabel('t_i')
ylabel('y_i')
text(2,1.12,'Experiment 1')
text(5.5,1.35,'Experiment 2')
text(9.2,1.42,'Experiment 3')
box on

% Analysis 2
fprintf('Analysis 2\n');
yavg = mean(ydata,1);
tmean = mean(t, 1);
X =  [ ones(1,p); tmean; tmean.^2 ]';
b = X \ yavg';
fprintf('b = %10.4g  %10.4g  %10.4g\n',b)
fprintf('cond(X^T*X) = %10.6g\n', cond(X'*X) )
% Analysis 3
fprintf('Analysis 3\n');
VI = zeros(p,p);
b = zeros(p,1);
for im = 1:m
    omega = diag( ones(1,n)/s(im)^2 );  
    X = [ ones(n,1) t(:,im) t(:,im).^2 ];
    y = ydata(:,im);
    % use Eq. 2.7.21a
    ya = [ b; y];
    Xa = [ eye(p,p); X ];
    Phi_inv = [ blkdiag( VI, omega ) ];
    b = (Xa'*Phi_inv*Xa) \ Xa'*Phi_inv*ya;
%     if( im > 1 )     % can use this alternatively for im > 1
%         lscov( Xa, ya, inv(Phi_inv) ) 
%     end
    fprintf('Experiment %2d:\n',im)
    fprintf('b = %10.4g  %10.4g  %10.4g\n',b)
    fprintf('cond(Xa^T*Phi_inv*Xa)  = %10.6g\n', cond(Xa'*Phi_inv*Xa) )
    fprintf('cond(X^T*omega*X + VI) = %10.6g\n', cond(X'*omega*X + VI) )
    VI = VI + X'*omega*X;  % by Eq. 2.1.17 for the next iteration
end
% print the 600 dpi figure to tiff format
fig=figure(1);
print(fig,fullfile(pwd, 'Figure2_9.tif'),'-dtiff','-r600')



