function [b,Yest,resid,SS,F,R2,Vb,VYest]=leastsqu(A,B,noyint);
% LEAST SQUARES REGRESSION ANALYSIS, USING MATRIX EQUATIONS
% from Draper & Smith (1981) Applied Regression Analysis, Wiley
% format: [b,Yest,resid,SS,F,R2,Vb,VY]=leastsqu(A,B);
%         or [b,Yest,resid,SS,F,R2,VY]=leastsqu(A,B,noyint);
% input: B=vector of observations of dependent variable
%        A=vector for independent variable
%        noyint=1, the calculation will be done with no y-intercept
%        noyint~=1, or leastsq(A,B) fits a y-intercept model
% output: b =    parameter vector [y-intercept;slope]
%         Yest = fitted estimates of dependent variable, Yest=X*b;
%         resid= residuals=Y-Yest
%         SS =   Sum of squares ANOVA table with DF and Mean squares
%                Partitioned into pure error & lack of fit
%                if repeated observations are present
%          F =   statistics and probabilities for SS table
%                [MSnum MSden dfnum dfden F prob(F)]
%                Test of regression and lack of fit if repeated X.
%          R2 =   R-squared;Vb= Variance-covariance matrix for b.
%          VYest =     Var (Yest)
% E. Gallagher, Environmental Sciences Program, UMASS/Boston
% Eugene.Gallagher@umb.edu, rev. 11/27/10
% Uses stixbox functions.
% uses ties.m, fprob, betai, betacm, gammln
[rA,cA]=size(A);
[rB,cB]=size(B);
% if rA*cA~=rB*cB
  %   error('DATA vectors must be the same length')
  % commented out for multiple regression
% end
% if max([min([rA cA]) min([rB cB])])~=1;
  %   error('DATA must be row or column vectors');
% end
if rA < cA       % change row vectors to col vectors
     A=A';
end
onesc=ones(rA,1);
if nargin>2;
     X=A;         % no y-intercept will be fitted.
else
     X=[onesc A];  % will fit a y-intercept
end
if rB < cB
     Y=B';
else
     Y=B;
end
b=X\Y;
if nargout==1
   return  % stop here if only y-intercept & slope are needed
end;
Yest=X*b;
resid=Y-Yest;
if nargout>3           % only calculated if ANOVA table needed
   n=length(Y);
%   R=X/(X'*X)*X';    % This is equivalent to R=X*inv(X'*X)*X'
                      % Draper & Smith use this R matrix, however
                      % this blew up the AT Matlab with N=211: too large
%   SSb1b0=Y'*(R-onesc*onesc'/n)*Y;
    SSb1b0=b'*X'*Y-Y'*onesc*onesc'*Y/n;     % Alternate to using R
% Here is the full ANOVA table from Draper & Smith (p. 85) with d.f.
   RegSS=b'*X'*Y;
   ResSS=Y'*Y-b'*X'*Y;
   TotSS=Y'*Y;
   % We can further partition the ANOVA table:
   Meanb0=n*(mean(Y)^2);
   % Now for the pure error sums of squares and mean squares, if there
   % are repeated measurements of the independent variable.
   if cA==1
      [T,R,i]=ties(A);             % ties.m is my own little m.file
   else
      T=[0 0]';
   end
   if sum(T)~=0
        % Use equations from page 36 of Draper & Smith
        k=1;
        SSPE=0;
        for j=1:length(T)
              tiedY=Y(i(k:k+T(j)-1));  % a column vector
              meanYj=mean(tiedY);
              SSPE=SSPE+sum((tiedY-ones(size(tiedY))*meanYj).^2);
              k=k+T(j);
        end
        if length(T)>1
             pedf=sum(T-ones(size(T)));
        else
             pedf=T-1;
        end
        MSPE=SSPE/pedf;
        LackfitSS=Y'*Y-b'*X'*Y-SSPE;
   else
        SSPE=0;
   % there were no repeated measurements of the independent variable
   end
   TotcSS=TotSS-Meanb0;
   % The F statistic for the regression:
   RF=SSb1b0/(ResSS/(n-length(b)));
   PRF=fprob(1,n-length(b),RF);  % Uses Gallagher's fprob.m
   if SSPE~=0
        SS=[SSb1b0 1 SSb1b0
           ResSS n-length(b) ResSS/(n-length(b))
           LackfitSS n-length(b)-pedf LackfitSS/(n-length(b)-pedf)
           SSPE pedf SSPE/pedf
           TotcSS n-1 TotcSS/(n-1)];
           % The lack of fit F statistic is LFF:
           LFF=(LackfitSS/(n-length(b)-pedf))/(SSPE/pedf);
           PLFF=fprob(n-length(b)-pedf,pedf,LFF);
           F=[SS(1,3) SS(2,3) SS(1,2) SS(2,2) RF PRF
              SS(3,3) SS(4,3) SS(3,2) SS(4,2) LFF PLFF];
           if PLFF>0.05
              s2=SS(2,3);    % if no lack of fit ResMS estimates s^2
           else
              s2=SS(4,3);   % if significant lack of fit, pure error MS
                            % estimates s^2
              disp('Warning: Significant lack of fit')
           end
   else
        SS=[SSb1b0 1 SSb1b0
           ResSS n-length(b) ResSS/(n-length(b))
           TotcSS n-1 TotcSS/(n-1)];
           F=[SS(1,3) SS(2,3) SS(1,2) SS(2,2) RF PRF];
           s2=SS(2,3);
   end
   % R^2 is easy to calculate:
   R2=(b'*X'*Y-Meanb0)/(Y'*Y-Meanb0);
   Vb=inv(X'*X)*s2;                   % Equation 2.3.2 in Draper & Smith
%  calculation of the variance of Yest (Draper & Smith, p. 28 & 83)
   [c,maxsize]=computer;
   Xo=X';
   [r,c]=size(Xo);
   if c^2<maxsize;          % use matrix approach on page 83
       VYest=diag(Xo'*inv(X'*X)*Xo*s2);
   else
        N=length(A);
        AmnA=A-ones(N,1)*mean(A);
        % See equations on p. 28-29
        VYest=s2*(ones(N,1)/N+AmnA.^2./sum(AmnA.^2));
   end
end