function P3_05A2
clear, clc, format short g, format compact
prob_title = (['Heat Transfer Correlation - Linear Regression 2']);
ind_var_name=['Re '];
dep_var_name=['Nu) '];
xyData=[5.624018	10.79958	0.8329091	-0.0544562
5.852202	11.13605	0.8241754	-0.0470916
6.042633	11.34805	0.8197798	-0.0418642
5.407172	10.43998	0.8415672	-0.058689
5.17615	10.03889	0.8586616	-0.0661398
4.743191	7.186144	5.505332	-0.5242486
4.563306	6.836259	5.509388	-0.5395681
4.22391	6.249975	5.525453	-0.5464528
3.893859	5.846439	5.609472	-1.237874
4.025352	4.811371	7.325149	-1.224176
3.686376	3.988984	7.371489	-1.276543
3.850148	4.437934	7.327123	-1.320507
4.54542	7.130099	4.67656	-0.3229639
4.60417	6.928538	5.225747	-0.491023
4.420045	6.142037	6.025866	-0.6694307
3.580737	4.00369	7.171657	-1.298283];
X=xyData(:,2:end);
y=xyData(:,1);
[m,n]=size(X);
freeparm=input(' Input 1 if there is a free parameter, 0 otherwise >  ');
[Beta, ConfInt,ycal, Var, R2]=MlinReg(X,y,freeparm);
disp([' Results,  ' prob_title]);
Res=[];
if freeparm==0, nparm = n-1; else nparm = n; end 
for i=0:nparm
          if freeparm, ii=i+1; else ii=i; end
          Res=[Res; ii Beta(i+1) ConfInt(i+1)];
end
disp('     Parameter No.   Beta      Conf_int');      
disp(Res);
sigmar=0;
Var=0;
for i=1:m
    Var=Var+ ((exp(ycal(i))-exp(y(i))))^2;
    sigmar=sigmar+((exp(ycal(i))-exp(y(i)))/exp(y(i)))^2;
end
disp(['  Variance ', num2str(Var/(m-(nparm+1)))]);
disp([' Relative variance ', num2str(sigmar/(m-(nparm+1)))]);  
plot(y,y-ycal,'*')
title(['Residual Plot, ' prob_title]) % residual plot
xlabel([dep_var_name '(Measured)'])
ylabel('Residual')
function [Beta, ConfInt, ycalc, Var, R2]=MlinReg(X,y,freeparm)
[m,n]=size(X); % m-number of rows, n-number of columns   
if freeparm				
	X=[ones(m,1) X]; % Add column of ones if there is a free parameter				
	npar=n+1;				
else					
	npar=n;				
end						
Beta=X\y;  % Solve X'Beta = Y using QR decomposition					
ycalc=X*Beta; % Calculated dependent variable values
Var=((y-ycalc)'*(y-ycalc))/(m-npar); % variance
ymean=mean(y);
R2=(ycalc-ymean)'*(ycalc-ymean)/((y-ymean)'*(y-ymean));%linear correlation coefficient
% Calculate the confidence intervals
A=X'*X;	
Ainv=A\eye(size(A)); %Calculate the inverse of the X'X matrix
tdistr95=[12.7062 4.3027 3.1824 2.7764 2.5706 2.4469 2.3646 2.306 2.2622 2.2281...
            2.2010 2.1788 2.1604 2.1448 2.1315 2.1199 2.1098 2.1009 2.093 2.086 2.0796... 
            2.0739 2.0687 2.0639 2.0595 2.0555 2.0518 2.0484 2.0452 2.0423 2.0395 2.0369...
            2.0345 2.0322 2.0301 2.0281 2.0262 2.0244 2.0227 2.0211 2.0195 2.0181 2.0167...
            2.0154 2.0141]; % 95 percent probability t-distr. values
if (m-npar)>45
   t=2.07824-0.0017893*(m-npar)+0.000008089*(m-npar)^2; % t for degrees of freedom > 45
else
   t=tdistr95(m-npar);
end
for i=1:npar					
   	ConfInt(i,1)=t*sqrt(Var*Ainv(i,i)); %confidence intervals					
end