function [varargout]=regress(varargin)
//
// (c) David Cousens 12/1/2006
// Performs a multulinear least squares fit of the dependent variable y a
// (nx1) vector on a matrix x of independent variables (nxm)with each variable
// in the columns of x. If weights w are specified for each value of y a 
// weighted fit is performed. By default the data is fitted with an intercept
// (m+1 fit parameters). An optional flag allows fitting with no intercept 
// (m parameters).
//
// [a[,[resid][,[yp],[cov][,[rr][,[Syx][,[SSt][,[SSe][,[SSr]]]]]]]]]]
//     =regress(y,x [,[w][,[flag]]])
//
// Input Variables
//
//  y -     dependent (real or complex)variable n x 1 vector 
//  x -     independent (real or complex)variables  n x m matrix with each variable
//          in the columns
//  sy -     standard deviation r(eal) of the dependent variable y nx1 vector-un weighted if
//          not specified
//  flag -  boolean -- %t(default): fit with intercept (m +1 fitted parameters)
//                  %f         : fir with no intercept (m fitted parameters)
//
// Output variables
//
//  a -     fitted paramteers (real or complex - (m+1) x 1 vector or m x 1 vector if flag=%f)
//  resid - residuals of the fitted parameters n x 1 vector
//  yp -    predicted values for the input data points
//  cov -   covariance matrix ((m+1) x (m+1) or m x m if flag=%f.) Errors in fitted parameters
//          are the diagonal elements
//  rr -    zero order correlation coefficients - (m+1) x (m+1) upper diagonal matrix
//          rr(1,1)= R coefficient of multiple correlation on y on columns of x
//          rr(1,k), k=2:m+1= correlation coefficient of y on x(:,k-1)
//          rr(j,k), k=2:m+1,j=2:k contains the correlation coefficient of x(:,j-1) on x(:,k-1)
//
//  Syx -   standard error of the estimate of y on x
//  SSt -   total Sum of squares 
//  SSe -   explained Sum of squares
//  SSr -   residual Sum of squares
//
// © David Cousens 12/1/2006
//
  [lhs,rhs]=argn(),
  if ~rhs then 
    // demo of regress multlinear regression
    	title_demo = [
		'';
		'Demo of regress()';
		'========================================';
		''];
		s_mat=[
    "at=2;bt=3;y0=3";
    "x=[];for j=1:5,x=[x;[(1:5)''+0.1*rand(5,1) j*ones(5,1)]];end";
    "y=y0+at*x(:,1)+bt*x(:,2)+0.5*rand(size(x,1),1,''normal'')";
    "scf(1);clf;";
    "plot3d((1:5),(1:5),matrix(y,5,5))";
    "xtitle(''data for regression of y on x1,x2'',''x1'',''x2'')";
    "xname(''data for regression of y on x1,x2'')";
    "sy=ones(y)";
    "flag=%t";
    "[a,resid,yp,cov,rr,Syx,SSt,SSe,SSr]=regress(y,x,sy,flag)";
    "disp('''')";
    "disp(''fitted parameters='')";
    "disp(a.'')";
    "disp(''predicted values'')";
    "disp(yp.'')";
    "disp(''residuals='')";
    "disp(resid.'')";
    "disp(''covariance matrix='')";
    "disp(cov)";
    "disp(''correlation coefficients-zero order='')";
    "disp(rr)";
    "disp(''     Syx          SSt          SSe          SSr'')";
    "disp([Syx,SSt,SSe,SSr])";
    "scf(2)";
    "subplot(1,2,1)";
    "plot(x(:,1),y,''b+'')";
    "plot(x(:,1),y-resid,''r--o'')";
    "xtitle(''regress y-x1'',''y'',''x1'')";
    "subplot(1,2,2)";
    "plot(x(:,2),y,''b+'')";
    "plot(x(:,2),yp,''r--o'')";
    "xtitle(''regress y-x2'',''y'',''x2'')";
    "xname(''individual plots of y on x1 and x2'')";
    "scf(3);clf";
    "plot3d((1:5),(1:5),matrix(yp,5,5))";
    "xtitle(''regression of y on x1,x2'',''x1'',''x2'')";
    "xname(''regression of y on x1,x2'')";
    "scf(4);clf";
    "plot3d((1:5),(1:5),matrix(resid,5,5))";
    "xtitle(''residuals of regression of y on x1,x2'',''x1'',''x2'')";
    "xname(''residuals of regression of y on x1,x2'');";
    ];
    write(%io(2),title_demo);
  	 write(%io(2),s_mat);
  	 write(%io(2),' ');
  	 execstr(s_mat);
  	 outlist=list(0);
    //[a,resid,yp,cov,rr,Syx,SSt,SSe,SSr]=return(a,resid,yp,cov,rr,Syx,SSt,SSe,SSr);
  else
    if rhs<2 then
      error('not enough input variables-must have y,x-running demo'),
    else
      y=varargin(1),
      x=varargin(2),
    end
    if rhs>2 then
      sy=varargin(3),
    else
      sy=ones(y),
    end
    if rhs>3 then, flag=varargin(4),else flag=%t,end
    //
    // Check the input variables and setup dimensions
    //
    if (type(x) <> 1)|(type(y) <> 1)|(type(sy) <> 1)|(type(flag) <> 4) then,
      error('y,x,sy variables must be real or complex and flag boolean'),
    end
    [ny,my]= size(y),
    if my<>1 then, error('y must be a column vector'), end
    [nx,mx]=size(x),
    [nsy,msy]=size(sy),
    if my<>1 then, error('sy must be a column vector'), end
    if (nx<>ny)&(nx<>nsy) then, error('column lengths of y,x,sy must be equal'), end
    n=nx;clear ny;clear ny;clear nsy;
    if flag then, m=mx,M=m+1, else m=mx,M=m,end
    clear mx;clear my;clear msy;
    //
    // Calculate the fitting matrix and vector
    //
    if flag then
      xp=[ones(n,1),x],
    else
      xp=x,
    end
    // apply the wieghtings
    yp=y./sy,
    for j=1:m, xp(:,j)=xp(:,j)./sy, end
    //
    b=xp'*y,
    X=xp'*xp,
    //
    // Invert the matrix and calculate parameters and stats
    //
    if det(X)<>0 then
      cov=inv(X),
      a=cov*b,
      yp=xp*a, // calculate the predicted values
      resid=y-yp, // residuals
      if lhs >=4 then
        SSr=(resid./sy)'*(resid./sy),
        yb=sum(yp)/X(1,1), // weighted mean value of ybar
        SSt=((y-yb)./sy)'*((y-yb)./sy), // variance of the data points about yb
        SSe=((yp-yb)./sy)'*((yp-yb)./sy), // variance of the estimates about yb
        check=((y-yp)./sy)'*((yp-yb)./sy), // should be zero
        check2=SSt-SSe-SSr,  //likewise
        Syx=sqrt(SSr/X(1,1)),  // standard error of estimate of y on x
        Sy=sqrt(SSt/n),
        //
        // calculate the correlation coefficients
        //
        R=sqrt(1-SSr/SSt), // coefficient of multiple correlation  of y on all x
        //
        // zero order correlation coefficients
        xb=sum(xp,'r')/X(1,1),// vector of means of the independent variables
        //
        Xb=[yb*ones(n,1)];for j=2:M,Xb=[Xb,xb(j)*ones(n,1)], end,
        if flag then, XmXb=[yp,xp(:,2:M)]-Xb, else, XmXb=[yp,xp], end,
        SSv=sum(XmXb.^2,'r')/X(1,1),// vector of standard deviations of the variables (row)
        Sxy=[];for j=1:m+1, Sxy=[Sxy,sum(XmXb(:,j).*XmXb(:,1),'r')/X(1.1)], end // vector of the
        rr=zeros(m+1,m+1),
        for k=2:m+1,
          for j=1:k // upper diagonal of matrix only since symmetric
            rr(j,k)=Sxy(k)/sqrt((SSv(j)*SSv(k))),
          end
        end
        rr(1)=R,
      end
      outlist=list(),
      if lhs>=1 then, outlist($+1)=a, end,
      if lhs>=2 then, outlist($+1)=resid,end,
      if lhs>=3 then, outlist($+1)=yp,end,
      if lhs>=4 then, outlist($+1)=cov, end
      if lhs>=5 then, outlist($+1)=rr, end,
      if lhs>=6 then, outlist($+1)=Syx,end,
      if lhs>=7 then, outlist($+1)=SSt, end,
      if lhs>=8 then, outlist($+1)=SSe, end,
      if lhs>=9 then, outlist($+1)=SSr, end,  
    else
      error('fitting matrix is not invertible'),
    end
  end
  varargout=outlist,
endfunction