function [b, SR, t] = mreg(X)
% [b, SR, t] = mreg(X)
% 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.

[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('Sum of squares due to regression on x1 = %8.2f\n',SR1);
fprintf('Total sum of squares                   = %8.2f\n\n',t);
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
   SRA = SRI - SRL;   %sum of squares added by this regression
   fprintf('Sum of squares added by x%g, = %8.2f\n',i, SRA);
   SRL = SRI;
end
b = BI;
SR = SRI;

