function P3_03B
clear, clc,	format short g, format compact
prob_title = (['Thermal Conductivity of Gaseous Propane ']);
ind_var_name=['Temp. (K)'];
dep_var_name=['Thermal Conductivity (W/m*K) '];
xyData=[ 231.07	1.14E-02
240	1.21E-02
260	1.39E-02
280	1.59E-02
300	1.80E-02
320	2.02E-02
340	2.26E-02
360	2.52E-02
380	2.78E-02
400	3.06E-02
420	3.34E-02
440	3.63E-02
460	3.93E-02
480	4.24E-02
500	4.55E-02
520	4.87E-02
540	5.20E-02
560	5.53E-02
580	5.86E-02
600	6.19E-02];
x=xyData(:,1);
y=xyData(:,2);
[m,n]=size(x);
freeparm=input(' Input 1 if there is a free parameter, 0 otherwise   ');
degree=input(' Enter the degree of the polynomial ');
[Beta, ycal,ConfInt, Var, R2, n]=PolyReg(x,y,degree,freeparm);
disp([' Results,' prob_title]);
Res=[];
for i=0:n-1
    if freeparm==0; 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)]);
pause					
plot(y,y-ycal,'*')	% residual plot				
title(['Residual plot, ' prob_title])					
xlabel([dep_var_name '(measured)'])					
ylabel('residual')
pause
%
%Plot of experimental and calculated data
%	
    for i=1:m
        index(i)=i;
    end
	plot(index,ycal, 'r-',index,y,'b+' )						
	title(['Calculated/experimental data ' prob_title])						
	xlabel(['Point No.'])						
	ylabel([dep_var_name])
function [Beta, ycal,ConfInt, Var, R2, n]=PolyReg(x,y,degree,freeparm)
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 freeparm==1
      n=degree+1;
   else
      n=degree;
   end
   m=size(x,1);
   for i=1:m					
      for j=1:n	
         if freeparm==1
      			p=j-1;
   			else
      			p=j;
   			end
			X(i,j)=x(i)^p; %Calculate powers of the independent variable and put them in the matrix X			
		end				
   end	
      	Beta=X\y;  % Solve X'Beta = Y using QR decomposition					
	ycal=X*Beta; % Calculated dependent variable values
	Var=sum((y-ycal)'*(y-ycal))/(m-n); % variance
    if (m-n)>45
        t=2.07824-0.0017893*(m-n)+0.000008089*(m-n)^2;
    else
        t=tdistr95(m-n);
    end
    A=X'*X;
    Ainv=A\eye(size(A)); %Calculate the inverse of the X'X matrix 
    for i=1:n					
        ConfInt(i,1)=t*sqrt(Var*Ainv(i,i)); %confidence intervals					
    end
    ymean=mean(y);
    R2=(ycal-ymean)'*(ycal-ymean)/((y-ymean)'*(y-ymean));%linear correlation coefficient