function P3_11A2
clear, clc, format short g, format compact
prob_title = (['Regression of Rate Data']);
ind_var_name=['Concentration (g-mol/dm^3  '];
dep_var_name=['Reactoion Rate, g-mol/(dm^3-s) '];
xyData=[1.25E-08	0.0002	0.000798
2.50E-08	0.0004	0.000595
4.05E-08	0.0006	0.0004
7.50E-09	0.00015	0.000849
2.80E-08	0.0004	0.000599
3.57E-08	0.00055	0.00045
2.86E-08	0.00045	0.000547
3.44E-08	0.0005	0.000498
2.44E-08	0.0004	0.000599];
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