%Example2_22.m Thermal conductivity example
clear all; format short g; 
T = [ 100	90	161	149	227	206	270	247	362	352 ]';
P = [ 545	276	602	275	538	274	550	274	522	272 ]';
y = [ 41.6	42.345	37.7875	39.5375	36.4975	37.3525	35.785	36.36	34.53	33.915 ]';
n = length( y ) ;

Va=[1 4 1 4 1 4 1 4 1 4];  
V=diag(Va); VI=inv(V); R=y'*VI*y;
%
% use the Model 5
X=[ones(10,1) T T.^2]; p=3; 
b=lscov(X,y,V);  R=(y-X*b)'*VI*(y-X*b);
fprintf('\nEst. Cov(b_ML)\n')
fprintf(' %10.4g %10.4g %10.4g \n', inv(X'*VI*X))
s2_est = (y - X*b)'*VI*(y - X*b)/(n-p)

Tvals = [ 60:5:375 ];
P_ML = inv(X'*VI*X);   % Eq. 2.6.5
nv = length( Tvals);
X1 = [ ones( nv, 1) Tvals' (Tvals').^2 ];
% X1 = X;
cov_y = X1*P_ML*X1';
cov_y_ii = diag( cov_y );
se_y = sqrt( cov_y_ii*s2_est )
CI = se_y*tinv(0.975,n-p)


residuals = y - X*b;

figure(1); clf;
plot(T(1:2:end), residuals(1:2:end), 'bx', T(2:2:end), residuals(2:2:end), 'go', 'linewidth', 1 ); figure(1)
axis( [50 400 -1.5 1.5 ] )
hold on
plot( Tvals, CI, 'k-', 'linewidth', 1 )
plot( Tvals, -CI, 'k-', 'linewidth', 1 , 'HandleVisibility', 'off')
legend('About 275W', 'About 550W', 'Conf Int.')
xlabel('Temperature, F')
ylabel('Residuals, Btu/hr-ft-F')
hold off
% print the 600 dpi figure to tiff format
fig=figure(1);
print(fig,fullfile(pwd, 'Figure2_8.tif'),'-dtiff','-r600')
