function [b, SR, t] = mreg(X)
% [b, SR, t] = mreg(X)  from Middleton (2000)
% mreg(X) performs a multiple linear regression on the data
% in the array X: this array has N rows (samples) and p columns:
% the first p-1 columns are measurements on independent
% variables (or squares, products, etc. calculated from
% such data).  The output consists of b, the vector of
% regression coefficients; SR, the sum of squares due to
% regression; and t, the total sum of squares.
% modified by L Braile 02/28/05
[N,p] = size(X);
m = mean(X);
D = X - ones(N,1)*m;  % convert X to deviations from mean
SS = D'*D;                   
                      % partition the matrix SS
S = SS(1:p-1, 1:p-1);
P = SS(1:p-1,p);
t = SS(p,p);
               % solve for regression on the first variable
b1 = P(1)/S(1,1);
b0 = m(p) - b1*m(1);
r  = P(1)/sqrt(t*S(1,1));
SR1 = r*r*t;
fprintf('\nAnalysis of Variance Table.\n\n');
fprintf('Source            d.f.   Sums of Squares  Mean Squares\n');
fprintf('------------------------------------------------------\n');
fprintf('Regression on x1   1     %8.4f        %8.4f\n',SR1,SR1);
SRL = SR1;	% last sum of squares due to regression
for i = 2:p-1 
   SI = S(1:i,1:i);   % redefine S,B,SR
   PI = P(1:i);
   BI = SI\PI;
   SRI = BI'*PI;      % sum of squares due to this regression
   MSRI = SRI/i;
   fprintf('This regression   %2.0f     %8.4f        %8.4f\n',i,SRI,MSRI);
   SRA = SRI - SRL;   %sum of squares added by this regression
   fprintf('Added by x%g,       1     %8.4f        %8.4f\n',i,SRA,SRA);
   SRL = SRI;
end
fprintf('Residuals         %2.0f     %8.4f        %8.4f\n',N-p,t-SRI,(t-SRI)/(N-p));
fprintf('----------------------------------------\n');
fprintf('Total             %2.0f     %8.4f\n',N-1,t);
b = BI;
SR = SRI;