%Example2_21.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;
% loop to evaluat the different models
for ii=1:6
   if ii == 1; X=ones(10,1); p=1;
       elseif ii==2; X=[ones(10,1) T]; p=2; R1=R;
       elseif ii==3; X=[ones(10,1) P]; p=2; R2=R; R=R1;
       elseif ii==4; X=[ones(10,1) T P]; p=3; R=R2;
       elseif ii==5; X=[ones(10,1) T T.^2]; p=3; R5=R; R=R2;
       else; X=[ones(10,1) T T.^2 T.^3]; p=4;
   end
Ftest = finv( 0.95, 1, n-p );
b=lscov(X,y,V); df=n-p; DR=R-(y-X*b)'*VI*(y-X*b); R=(y-X*b)'*VI*(y-X*b);
s2=R/df; F=DR/s2; Conx=cond(X); B(ii,:)=[ii p df R s2 DR F Ftest Conx];

end
fprintf('Model p   D.o.F.   R Est.   var.          DR     F=DR/s^2      F_test    Cond(X)\n');
fprintf('%3d  %3d   %3d  %7.4f   %7.4f   %9.4f   %9.2f   %9.2f  %8g\n', B')
%
% use the Model 5
X=[ones(10,1) T T.^2]; p=3; R=R2;
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)
figure(1); clf;
plot(T(1:2:end), y(1:2:end), 'bx', T(2:2:end), y(2:2:end), 'go', 'linewidth', 1.0 ); figure(1)
hold on
Tvals = [ 60:5:375 ];
plot( Tvals, b(1) +b(2)*Tvals + b(3).*Tvals.^2, 'k-', 'linewidth', 1.0)
legend('About 275W', 'About 550W', 'Model 5')
xlabel('Temperature, F')
ylabel('Thermal Conductivity, Btu/hr-ft-F')
hold off
% print the 600 dpi figure to tiff format
fig=figure(1);
print(fig,fullfile(pwd, 'Figure2_7.tif'),'-dtiff','-r600')