function W = hampel(res, scale, param)
% function W = hampel(res, scale, param)
%
% computes Hampel's weight function for robust regression:
% 
% hampel(x) = 	{ 1,			 |x| < param(1),
%		{ a/|x|,		 a <= |x| < b,
%		{ a/|x|*(c-|x|)/(c-b),	 b <= |x| <  c,
%		{ 0,			 |x| >= c.
%
%	must have 0 <= a <= b <= c
%
% 
% arguments:
%		res:	vector of residuals
%		scale:	robust estimate of scale, such as MAD
%		param:	parameter vector for the Hampel function
%			param = [a, b, c].
%			A standard choice is [2 4 8], which is the
%			default.   The resulting influence
%			function is nearly identical to the biweight 
%			with parameter 8.
%
%
% returns:
%		W:	the vector of Hampel weights
%
%
% P.B. Stark  stark@stat.berkeley.edu
% 11 July 1997.

if ((param(1) < 0) | (param(2) < param(1)) | (param(3) < param(2))),
	error([' illegal choice of parameters in Hampel: ' ...
		num2str(param) ]')
end

if (nargin == 3),
	a = param(1);
	b = param(2);
	c = param(3);
else
	a = 2;
	b = 4;
	c = 8;
end

W = ones(size(res));
W(abs(res) >= a) = a*(abs(res(abs(res) >= a))).^(-1);
W(abs(res) >= b) = W(abs(res) >= b).*(-abs(res(abs(res) >= b)) +c) ...
		/(c-b);
Z = zeros(size(res));
W(abs(res) > c) = Z(abs(res) > c);

return;