function [tmcd,smcd,weights,initmean,initcov,res,raww,Hopt]=unimcd(y,h)
%UNIMCD computes the MCD estimator of a univariate data set. This
% estimator is given by the subset of h observations with smallest variance.
% The MCD location estimate is then the mean of those h points,
% and the MCD scale estimate is their standard deviation. Based on these raw estimates,
% a reweigthing step is applied as in the FASTMCD algorithm (see mcdcov.m).
% We recommend to use the function mcdcov.m, which calls unimcd.m.
%
% The MCD method was introduced in:
%
% Rousseeuw, P.J. (1984), "Least Median of Squares Regression,"
% Journal of the American Statistical Association, Vol. 79, pp. 871-881.
%
% The algorithm to compute the univariate MCD is described in
%
% Rousseeuw, P.J., Leroy, A., (1988), "Robust Regression and Outlier Detection,"
% John Wiley, New York.
%
% This function is part of LIBRA: the Matlab Library for Robust Analysis,
% available at:
% http://wis.kuleuven.be/stat/robust.html
%
% Written by: Katrien Van Driessen
%Revision by S. Verboven
%Last update
ncas=length(y);
len=ncas-h+1;
if len==1
tmcd=mean(y);
smcd=sqrt(var(y));
weights=ones(length(y),1);
else
[y,I]=sort(y);
ay(1)=sum(y(1:h));
for samp=2:len
ay(samp)=ay(samp-1)-y(samp-1)+y(samp+h-1);
end
ay2=ay.^2/h;
sq(1)=sum(y(1:h).^2)-ay2(1);
for samp=2:len
sq(samp)=sq(samp-1)-y(samp-1)^2+y(samp+h-1)^2-ay2(samp)+ay2(samp-1);
end
sqmin=min(sq);
ii=find(sq==sqmin);
Hopt = I(ii:ii+h-1);
ndup=length(ii);
slutn(1:ndup)=ay(ii);
initmean=slutn(floor((ndup+1)/2))/h; %initial mean
initcov=sqmin/(h-1); %initial variance
% calculating consistency factor
res=(y-initmean).^2/initcov;
sortres=sort(res);
factor=sortres(h)/chi2inv(h/ncas,1);
initcov=factor*initcov;
res=(y-initmean).^2/initcov; %raw_robdist^2
quantile=chi2inv(0.975,1);
weights=res<quantile;
raww=weights; %raw_weights
%reweighting procedure
if size(weights,1)~=size(y,1)
weights=weights';
end
tmcd=sum(y.*weights)/sum(weights);
smcd=sqrt(sum((y-tmcd).^2.*weights)/(sum(weights)-1));
[I2,J] = sort(I);
weights=weights(J); %rew_weights
end