function b = regcoefs(X,y,cutoff)

% REGCOEFS   Linear regression coefficients (after studentizing outliers).
%
%  b = regcoefs(X,y)
%  b = regcoefs(X,y,cutoff)
%
%  X is a matrix of predictors, Y is a column vector.
%  A column of ones is added to the front of X making for an intercept term.
%  CUTOFF must be a positive number (e.g. 3.0) or zero. When CUTOFF is
%  set to zero (the default),  b=regcoefs(X,y)  is simply a shorthand for
%    b = [ones(size(X,1),1),X] \ y
%  When CUTOFF>0, however, two regressions are performed to increase the
%  robustness of the coefficient estimates. The first regression
%  identifies and studentizes outliers and the second one calculates B coefs.
%  An "outlier" is a Y value whose standardized residual after the (first) 
%  regression is greater in absolute value than CUTOFF. All outliers, if
%  any, are studentized and then a second regression is performed to 
%  calculate the final robust regression coefficients B.
%
%  See also \, MLDIVIDE, R2SIGNIF.

% Original coding by Alex Petrov, Ohio State University
% $Revision: 1.1 $  $Date: 2000/10/26 15:20 $
%
% Part of the utils toolbox version 1.1 for MATLAB version 5 and up.
% http://alexpetrov.com/softw/utils/
% Copyright (c) Alexander Petrov 1999-2006, http://alexpetrov.com
% Please read the LICENSE and NO WARRANTY statement in ../utils_license.m

if (nargin < 3)
  cutoff = 0 ;    % default: do not studentize residuals
end

X = [ones(size(X,1),1),X] ;    % add an intercept term

b = X \ y ;      % first regression

if (cutoff > 0)
  yhat = X * b ;
  resid = y - yhat ;
  c = cutoff * std(resid) ;    % cutoff in standard deviation units
  outl = find(abs(resid) > c) ;
  if (~isempty(outl))
    y(outl) = yhat(outl) + c .* sign(resid(outl)) ;    % studentize
    b = X \ y ;  % second and final regression
  %else there are no outliers and hence no need for second regression
  end
%else no studentizatoin is requested -- just return B from the first regr.
end
  
%-- Return B
%%%%% End of REGCOEFS.M
