function [beta,s] = lmsreg(x,y,sub,rndst)

% function [beta,s] = lmsreg(x,y,sub,rndst)
% Least Median of Squares regression estimator, computed by subsampling.
%
%INPUT ARGUMENTS
%   x:      dependent variables (n x p)
%   y:      independent variable (n x 1)
%   sub:    number of subsamples (scalar)
%   rndst:  random generator initial state (scalar)
%
%OUTPUT ARGUMENTS
%   beta:   LMS regression estimate (p x 1)
%   s:      Associated scale estimate (calibrated for Normal consistency)

rand('state',rndst)
[n,p] = size(x);
cp = nchoosek(n,p);
if (3*cp)<sub, sub=3*cp; end

beta = zeros(p,1);
bs = zeros(p,1);
s2 = median(y.^2);
isub = zeros(1,n);
res2 = zeros(size(y));
for k = 1:sub,
    isub = randperm(n);
    xs = x(isub(1:p),:); ys = y(isub(1:p));
    if det(xs)~=0,
	bs = xs\ys;
	res2 = (y-x*bs).^2;
	if mean(res2>=s2)<0.5,
		s2 = median(res2);
		beta = bs;
	end
    end
end
s = 1.4826*(1+5/(n-p))*sqrt(s2);