% MATLAB script file for fitting the linear model, y = m1*x1 + m2*x2 + ...
% + mn*xn + b to an experimental data set where each data point is of the
% form (x1,x2,...,xn,y_hat).
%
% See Supplemental Unit S3 of "A First Course on Kinetics and Reaction
% Engineering" and/or the MATLAB documentation for additional information
% or if questions arise.
%
% Carl R. F. Lund, 4/22/13
%   last rev. 3/28/14
%
% Required input (must be present in the MATLAB workspace): 
%   x (matrix) - x(i,j) is the value of set variable j for experimental
%       data point i; the last column must contains the value 1 in every
%       row
%   y_hat (column vector) - y_hat(i) is the measured value of the response
%       variable for data point i

% Perform the regression analysis
[c,cint,r,rint,stats] = regress(y_hat,x);

% Calculate the 95% confidence limits from the 95% confidence interval
clim = (cint(:,2) - cint(:,1))/2.0;

% Calculate the responses predicted by the final model
y_resp = x*c;

% Report the results, separating the slopes from the intercept
format shortE;
r_squared = stats(1)
m = c(1:length(c)-1)
m_u = clim(1:length(clim)-1)
b = c(length(c))
b_u = clim(length(clim))

% Generate the parity plot
figure;
plot(y_hat,y_hat,'r',y_hat,y_resp,'o')
xlabel('Measured Response')
ylabel('Predicted Response')
title('Parity Plot')

% Generate the residuals plots
n_set = size(x,2) - 1;
for n = 1:n_set
    figure;
    plot(x(:,n),zeros(length(y_hat),1),'r',x(:,n),r,'o')
    xlabel(['x',num2str(n)])
    ylabel('Residual')
    title(['Residuals Plot vs. x',num2str(n)])
end
    
% End of FitLinSR.m