function P2_05B
clear, clc, format short g, format compact
prob_title = (['VP Benzene with Clapeyrons Equation']);
ind_var_name=['1/T (1/K)'];
dep_var_name=['log (Vapor Pressure) '];
%fname=input('Please enter the data file name  >  ');
%xyData=load(fname);
xyData=[3.936212	0.0034473
4.187182	0.003307
4.351874	0.0032135
4.483787	0.0031379
4.590541	0.003076
4.677342	0.0030254
4.744379	0.0029861
4.804923	0.0029504
4.863234	0.0029159
4.909957	0.0028882
4.96069	0.0028579
5.00877	0.0028291
5.044736	0.0028075
5.079688	0.0027865
5.123133	0.0027602
5.159958	0.002738
5.192568	0.0027181
5.224274	0.0026988
5.255417	0.0026799
5.286748	0.0026607
5.317938	0.0026419
5.349161	0.0026225];
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);
disp(['  Variance ', num2str(Var)]);
disp(['  Correlation Coefficient ', num2str(R2)]);
plot(y,y-ycal,'*')
title(['Residual Plot, ' prob_title]) % residual plot
xlabel([dep_var_name '(Measured)'])
ylabel('Residual')
pause
plot(X(:,1),ycal, 'r-',X(:,1),y,'b+' )
title(['Calculated/Experimental Data ' prob_title])
xlabel([ind_var_name])
ylabel([dep_var_name])
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