function [beta,s,w,tn] = rewls(x,y,beta0,s0)

%function [beta,s,w,tn] = rewls(x,y,beta0,s0)
%
%Computes the one-step REWLSE using hard-rejection weights.
%
%INPUT ARGUMENTS:
%
%	x:      explanatory variables (n by p matrix)
%	y:      response variables (n by 1)
%	beta0:	initial regression estimator (p by 1)
%	s0:     initial scale estimator
%
%OUTPUT ARGUMENTS
%
%	beta:	REWLS regression estimator (p by 1)
%	s:      REWLS scale estimator
%	w:      REWLS weights (n by 1)
%	tn:     REWLS adaptive threshold value
%

n = length(y);
ra = abs(y-x*beta0)/s0;
raord = sort(ra);
dif = (2*normcdf(raord)-1) - (0:(n-1))'/n;

eta = 2.5;
i = find(raord>=eta);
if isempty(i)
   dn = 0;
else
   dn = max(dif(i));
end

if dn>0,
   tn = raord(n-floor(n*dn));
else
   tn = Inf;
end

w = ra<tn;
if sum(w)==0
   beta = beta0;
   s = 0;
else
   beta = inv(x'*diag(w)*x)*x'*(w.*y);
   s = sqrt(w'*((y-x*beta).^2)/(sum(w)-1));
end