function [a,b,R,yhat,resid] = linreg(x,y)
%  LINREG
%      [a,b,R,yhat,resid] = linreg(x,y) performs a linear regression to fit
%            y = a*x + b in the least squares sense.  Observations
%      should appear in y, the independent variable in x.  Returned value 
%      R is the correlation coefficient (not r^2).
%      If x and y are both NxM matrices, then M linear regressions are
%      performed fitting columns of y to corresponding columns of x.
%      a,b, and R are then rowvectors of length M.
%
%      The optional return arguments yhat and resid are the fits (predicted
%      y values) and the residuals (y-yhat).  Both yhat and resid are NxM
%      matrices.
%
%


%   Roger McNichols
%   Texas A&M University
%   11/10/94

% Check to make sure matrices are of same size.  If not an error will occur.
z = x+y;


% calculate regression coefficients

[n,m] = size(x);

xbar = ones(n,1)*mean(x);
ybar = ones(n,1)*mean(y);

Sxx = sum(x.^2) - 1/n*(sum(x)).^2;

Syy = sum(y.^2) - 1/n*(sum(y)).^2;

Sxy = sum(x.*y) - 1/n*(sum(x).*sum(y));

a = Sxy./Sxx;

b = mean(y) - a.*mean(x);

R = Sxy./((Sxx.^.5).*(Syy.^.5));


%  Calculate residuals, fits, etc.

yhat = ones(n,1)*b + (ones(n,1)*a).*x;

resid = y-yhat;