ltsregres.m

function [rew,raw] = ltsregres(x,y,varargin)

%LTSREGRES carries out least trimmed squares (LTS) regression, introduced in
%   
%   Rousseeuw, P.J. (1984), "Least Median of Squares Regression,"
%   Journal of the American Statistical Association, Vol. 79, pp. 871-881.
% 
% The LTS regression method minimizes the sum of the h smallest squared 
% residuals, where h must be at least half the number of observations. The 
% default value of h is roughly 0.75n (n is the total number of observations), 
% but the user may choose any value between n/2 and n.
%
% To compute the LTS estimator, the FAST-LTS algorithm is used.
% Reference: 
%  Rousseeuw, P.J. and Van Driessen, K. (2006), 
%  "Computing LTS Regression for Large Data Sets", Data Mining and
%  Knowledge Discovery, 12, 29-45.
%
%  Also available at http://www.agoras.ua.ac.be/
%
% The LTS regression method is intended for continuous variables, and assumes 
% that the number of observations n is at least 2 times the number of 
% regression coefficients p. If p is too large with respect to n, it is  
% better to first reduce p by variable selection or principal components 
% (see rpcr.m, rsimpls.m). The response variable should be univariate, otherwise 
% robust multivariate regression should be performed (see mcdregres.m).
%
% The LTS is a robust method in the sense that the estimated regression
% fit is not unduly influenced by outliers in the data, even if there are
% several outliers. Due to this robustness, we can detect outliers by their
% large LTS residuals.
%
% Required input arguments: 
%    x : Data matrix of explanatory variables (also called 'regressors'). 
%        Rows of x represent observations, and columns represent variables.  
%        Missing values (NaN's) and infinite values (Inf's) are allowed, since observations (rows) 
%        with missing or infinite values will automatically be excluded from the computations.
%    y:  A vector with n elements that contains the response variables.
%        Missing values (NaN's) and infinite values (Inf's) are allowed, since observations (rows) 
%        with missing or infinite values will automatically be excluded from the computations.
%
% Optional input arguments:           
%   intercept : If 1, a model with constant term will be fitted (default),
%               if 0, no constant term will be included. 
%   intadjust : If 1, the intercept adjustment will be applied in each step
%               of the algorithm. These calculations need substantially more 
%               computation time than intadjust=0, which is the default value.
%           h : The number of observations that have determined the least 
%               trimmed squares estimator. Any value between n/2 and n may be specified.
%       alpha : (1-alpha) measures the fraction of outliers the algorithm should 
%               resist. Any value between 0.5 and 1 may be specified. (default = 0.75)
%      ntrial : Number of initial subsets drawn. Its default value is 500.
%       plots : If equal to one, a menu is shown which allows to draw several plots,
%               such as residual plots and a regression outlier map. (default)
%               If the input argument 'classic' is equal to one, the classical
%               plots are drawn as well.
%               If 'plots' is equal to zero, all plots are suppressed.
%               See also makeplot.m
%     classic : If equal to one, classical least squares regression will be performed,
%               see ols.m (default = 0).
%
% Input arguments for advanced users:
%     Hsets : Instead of random trial h-subsets (default, Hsets = []), Hsets makes it possible to give certain
%             h-subsets as input. Hsets is a matrix that contains the indices of the observations of one
%             h-subset as a row.
%
% I/O:  result=ltsregres(x,y,'plots',0,'intercept',0)
%       [rew,raw] = ltsregres(x,y)
%  The user should only give the input arguments that have to change their default value.
%  The name of the input arguments needs to be followed by their value.
%  The order of the input arguments is of no importance.
%
% The output consists of two structures 'rew' and 'raw' containing the 
% following fields:
%
%    raw.coefficients : Vector of raw LTS coefficient estimates (including the 
%                       intercept, when options.intercept=1).
%          raw.fitted : Vector like y containing the raw fitted values of the response.
%             raw.res : Vector like y containing the raw residuals from the regression.
%           raw.scale : Scale estimate of the raw residuals.
%       raw.objective : Objective function of the LTS regression method, i.e. the sum 
%                       of the h smallest squared raw residuals.
%              raw.wt : Vector like y containing weights that can be used in a weighted
%                       least squares. These weights are 1 for points with reasonably
%                       small raw residuals, and 0 for points with large raw residuals.
%
% 	        rew.slope : Vector of the slope coefficients obtained after reweighting.
%             rew.int : The intercept.
%          rew.fitted : Vector like y containing the fitted values of the response
%                       after reweighting.
%             rew.res : Vector like y containing the residuals from the weighted
%                       least squares regression.
%           rew.scale : Scale estimate of the reweighted residuals.
%        rew.rsquared : Robust version of R squared. This is 1 minus the fraction:
%                       (sum of the quan smallest squared residuals) over (sum of 
%                       the quan smallest (y-loc)^2), where the denominator
%                       is minimized over loc. Note that loc is not subtracted from
%                       y if intercept = 0 in the call to ltsregres.	
%               rew.h : The number of observations that have determined the LTS estimator, 
%                       i.e. the value of h. 
%       rew. Hsubsets : A structure that contains Hopt and Hfreq:
%                        Hopt  : The subset of h points whose covariance matrix has minimal determinant, 
%                                ordered following increasing robust distances.
%                        Hfreq : The subset of h points which are the most frequently selected during the whole
%                                algorithm.
%           rew.alpha : (1-alpha) measures the fraction of outliers the algorithm should 
%                       resist. 
%              rew.rd : The robust distances for the observations of the design matrix, 
%                       based on the MCD estimator (mcdcov.m)  
%            rew.resd : Vector like y containing the standardized residuals
%                       from the weighted least squares regression.
%         rew.weights : Vector like y containing weights that have been used in a weighted
%                       least squares. These weights are 1 for points with reasonably
%                       small raw residuals, and 0 for points with large raw residuals.
%          rew.cutoff : Structure which contains cutoff values for the robust distances computed by mcdcov.m,
%                       and for the standardized residuals. 
%            rew.flag : Vector like y containing flags based on the reweighted regression.
%                       These flags determine which observations can be considered as 
%                       outliers.
%          rew.method : Character string naming the method (Least Trimmed Squares).	
%           rew.class : 'LTS'
%         rew.classic : If the input argument 'classic' equals 1, this structure contains the 
%                       results of a classical least squares regression
%               rew.X : If x is univariate, same as the input x in the call to ltsregres,
%                       without rows containing missing or infinite values. 
%               rew.y : If x is univariate, same as the input y in the call to ltsregres,
%                       without rows containing missing or infinite values.

%
% This function is part of LIBRA: the Matlab Library for Robust Analysis,
% available at: 
%              http://wis.kuleuven.be/stat/robust.html
%
% Version 22/12/2000, 
% Written by Katrien Van Driessen and Randy Brenkers
% Revisions by Sabine Verboven, Sanne Engelen, Nele Smets
% Last update: 06/07/2004,  03/07/2006 (profiler tips --> & -> &&, | --> ||)
% Last revision:  03/07/2006

%MATLAB: to avoid singularMatrix in line 538 
warning off 
%The maximum value for p = number of variables
pmax=20;
%The maximum value for n = number of observations
nmax=50000;
%To change the number of subdatasets and their size, the values of maxgroup
%and nmini can be changed
maxgroup=5;
nmini=300;
%The number of iteration steps in stages 1,2 and 3 can be changed
% by adapting the parameters csteps1, csteps2, and csteps3.
csteps1=2;
csteps2=2;
csteps3=100;

pmax1=pmax+1;
pmax2=pmax*pmax;
nvm11=pmax*pmax1;
nvm12=pmax1*pmax1;
km10=10*maxgroup;
nmaxi=nmini*maxgroup;
maxmini=fix(((3*nmini-1)/2)+1);
% dtrial : number of subsamples if not all (p+1)-subsets will be considered. 
dtrial=500;

[n,p]=size(x);
[m,q]=size(y);
if q~=1
    if m==1
        y=y';
    else
        error('y is not one-dimensional.');
    end
end
na.x=~isfinite(x*ones(p,1));
na.y=~isfinite(y);
if size(na.x,1)~=size(na.y,1)
    error('Number of observations in x and y not equal.');
end
% Observations with missing or infinite values are ommitted. 
ok=~(na.x|na.y);
x=x(ok,:);
y=y(ok,:);
dx=size(x);
dy=size(y,1);
n=dx(1);
% Some checks are now performed.
if n == 0
    error('All observations have missing values!');
end
if n > nmax
    error(['The program allows for at most ' int2str(nmax) ' observations.']);
end
%internal variables and default values
seed=0;
alfa=0.75;
hdef=quanf(alfa,n,p+1); %default value of h
hmin=quanf(0.5,n,p+1); %minmal value of h
default=struct('intercept',1,'intadjust',0,'alpha',alfa,'h',hdef,'plots',1,'ntrial',dtrial,'classic',0,'Hsets',[]);
list=fieldnames(default);
options=default;
IN=length(list);
i=1;
counter=1;
%Reading optional inputarguments
if nargin > 3
    %
    % placing inputfields in array of strings
    %
    for j=1:nargin-2
        if rem(j,2)~=0
            chklist{i}=varargin{j};
            i=i+1;
        end
    end 
    dummy=sum(strcmp(chklist,'h')+2*strcmp(chklist,'alpha'));
    switch dummy
        case 0 % defaultvalues should be taken
            alfa=options.alpha;
            h=options.h;
        case 3
            error('Both inputarguments alpha and h are provided. Only one is required.')
    end
    %
    % Checking which default parameters have to be changed
    % and keep them in the structure 'options'.
    %
    while counter<=IN 
        index=strmatch(list(counter,:),chklist,'exact');
        if ~isempty(index) % in case of similarity
            for j=1:nargin-3 % searching the index of the accompanying field
                if rem(j,2)~=0 % fieldnames are placed on odd index
                    if strcmp(chklist{index},varargin{j})
                        I=j;
                    end
                end
            end
            options=setfield(options,chklist{index},varargin{I+1});
            index=[];
        end
        counter=counter+1;
    end
    if dummy==1% checking inputvariable h
        if options.h < hmin                                                      
            error(['The LTS must cover at least ' int2str(hmin) ' observations.'])
        elseif options.h > n
            error('h is greater than the number of non-missings and non-infinites.')
        elseif options.h < p
            error(['h should be larger than the dimension ' int2str(p) '.'])
        elseif options.h==0
            options.h=h;
        end
        options.alpha=options.h/n;
    elseif dummy==2
        if options.alpha < 0.5
            options.alpha=0.5;
            mess=sprintf(['Attention : Alpha should be larger than 0.5. \n',...
                    'It is set to 0.5.']);
            disp(mess)
        end
        if options.alpha > 1
            options.alpha=0.75;
            mess=sprintf(['Attention : Alpha should be smaller than 1.\n',...
                    'It is set to 0.75.']);
            disp(mess)
        end
        if options.alpha == 0
            options.alpha=0.75;
        end
        options.h=quanf(options.alpha,n,p);
    end
end
intercept=options.intercept; %if 1 intercept in the model
h=options.h; %number of regular data points on which estimates are based (=[alpha * n])
plots=options.plots; %relevant plots if equal to 1
alfa=options.alpha; %proportion of regular data points
ntrial=options.ntrial; %number of subsets to be taken in the first step
classic=options.classic; %classic least squares regression?
bestobj=inf; %best objective value until 'now' -> goal is to be as small as possible
intadj=options.intadjust; %intercept adjustment needed if set to 1
Hsets = options.Hsets;
if ~isempty(Hsets)
    Hsets_ind = 1;
else
    Hsets_ind = 0;
end

if classic
    res.classic=ols(x,y,'plots',0,'intercept',intercept);
else 
    res.classic=0;
end

if intercept == 1
    dx=dx+[0 1];
    x=cat(2,x,ones(n,1));
end
p=dx(2);
if n < p
   error('Need more observations than variables.');
end

if p > pmax
   error(['The program allows for at most ' int2str(pmax) ' variables.'])
end

rk=rank(x);
if rk < p
   error('x is singular');
end

if h == n 
    res.method='Least Squares Regression.';
    [Q,R]=qr(x,0);
    z=R\(Q'*y);
    raw.coefficients=z;
    residuals=y-x*z;
    raw.res=residuals;
    fitted=x*raw.coefficients;
    raw.fitted=fitted;
    s0=sqrt(sum(residuals.^2)/(n-p));
    if abs(s0) < 1e-7
        weights=abs(residuals)<=1e-7;
        raw.wt=weights;
        raw.scale=0;
        res.scale=0;
        res.coefficients=raw.coefficients;
        raw.objective=0;
    else
        sor=sort(residuals.^2);
        raw.objective=sum(sor(1:h));
        raw.scale=s0;
        weights=abs(residuals/s0)<=norminv(0.9875);
        raw.wt=weights;
        [Q,R]=qr(x(weights==1,:),0);
        z=R\(Q'*y(weights==1));
        res.coefficients=z;
        fitted=x*res.coefficients;
        residuals=y-x*z;
        res.scale=sqrt(sum(weights.*(residuals.^2))/(sum(weights)-1));
        weights=abs(residuals/res.scale)<=norminv(0.9875);
    end
    if intercept
        s1=sum(residuals.^2);
        center=mean(y);
        sh=sum((y-center).^2);
        res.rsquared=1-s1/sh;
    else
        s1=sum(residuals.^2);
        sh=sum(y.^2);
        res.rsquared=1-s1/sh;
    end
    if res.rsquared > 1
        res.rsquared=1;
    elseif res.rsquared < 0
        res.rsquared=0;
    end
    res.Hsubsets.Hopt=1:n;
    res.Hsubsets.Hfreq=1:n;
    if abs(s0) < 0
        res.method=strvcat(res.method,'An exact fit was found!');
    end
    stdres=residuals/res.scale;
    cutoff.resd=sqrt(chi2inv(0.975,1));
    raw=struct('coefficients',{raw.coefficients},'fitted',{raw.fitted},'res',{raw.res},'scale',{raw.scale},...
        'objective',{raw.objective},'wt',{raw.wt});
    rew=struct('slope',{res.coefficients(1:p)},'int',{0},'fitted',{fitted},'res',{residuals},...
        'scale',{res.scale},'rsquared',{res.rsquared},'h',{h},'alpha',{alfa},'resd', {stdres},...
        'rd',{NaN},'cutoff',{cutoff},'flag',{NaN},'weights',{raw.wt},'Hsubsets',{res.Hsubsets},...
        'method',{res.method},'class',{'LTS'},'classic',{res.classic},'X',{x},'y',{y});
    if intercept
        rew=setfield(rew,'int',res.coefficients(p));
        rew=setfield(rew,'slope',res.coefficients(1:p-1));
    end
    if plots && classic 
        mcdres=mcdcov(x,'h',h,'plots',0);
        if -log(abs(det(mcdres.cov)))/size(data,2)> 50
            res.rd=NaN;
        else 
            res.rd=mcdres.rd;
        end
        cutoff.rd=mcdres.cutoff.rd;  
        cutoff.md=mcdres.cutoff.md;
        flags=abs(stdres)<=cutoff.resd;
        rew=setfield(rew,'rd', res.rd);
        rew=setfield(rew,'flag',flags);
        rew=setfield(rew,'cutoff',cutoff);
        try
            makeplot(rew,'classic',1)
        catch %output must be given even if plots are interrupted 
            %> delete(gcf) to get rid of the menu 
        end
    elseif plots 
        mcdres=mcdcov(x,'h',h,'plots',0);
        if -log(abs(det(mcdres.cov)))/size(x,2)> 50
            res.rd=NaN;
        else 
            res.rd=mcdres.rd;
        end
        cutoff.rd=mcdres.cutoff.rd;  
        cutoff.md=mcdres.cutoff.md;
        flags=abs(stdres)<=cutoff.resd;
        rew=setfield(rew,'rd', res.rd);
        rew=setfield(rew,'cutoff',cutoff);
        rew=setfield(rew,'flag',flags);
        try
            makeplot(rew)
        catch %output must be given even if plots are interrupted 
            %> delete(gcf) to get rid of the menu 
        end
    end
    return
end

if p < 5
    eps=1e-12;
elseif p <= 8
    eps=1e-14;
else
    eps=1e-16;
end

% standardization of the data
xorig=x;
yorig=y;
data=[x y];
if ~intercept
    datamed=repmat(0,1,p+1);
    datamad=median(abs(data)).*1.4826;
    for i=1:p+1
        if abs(datamad(i)) <= eps
            datamad(i)=sum(abs(data(:,i)));
            datamad(i)=(datamad(i)/n)*1.2533;
            if abs(datamad(i)) <= eps;
                error('The MAD of some variable is zero');
            end
        end   
    end
    x=x./repmat(datamad(1:p),n,1);
    y(:,1)=y(:,1)./datamad(p+1);
else
    datamed=median(data);
    datamed(p)=0;
    datamad(p)=1;
    for i=1:p+1
        if i ~= p
            datamad(i)=median(abs(data(:,i)-datamed(i)))*1.4826;
            if abs(datamad(i)) <= eps
                datamad(i)=sum(abs(data(:,i)-datamed(i)));
                datamad(i)=(datamad(i)/n)*1.2533;
                if abs(datamad(i)) <= eps
                    error('The MAD of some variable is zero');
                end
            end
        end
    end
    x=(x-repmat(datamed(1:p),n,1))./repmat(datamad(1:p),n,1);
    y(:,1)=(y(:,1)-datamed(p+1))./datamad(p+1);
end

res.method='Least Trimmed Squares Regression.';
al=0;
teller = zeros(1,n+1);
if Hsets_ind
    csteps = csteps1;
    inplane = NaN;
    fine = 0;
    part =  0;
    final = 1;
    tottimes = 0;
    nsamp = size(Hsets,1);
    obsingroup = n;
else
    if n >= 2*nmini  
        maxobs=maxgroup*nmini;
        if n >= maxobs
            ngroup=maxgroup;
            group(1:maxgroup)=nmini;
        else
            ngroup=floor(n/nmini);
            minquan=floor(n/ngroup);
            group(1)=minquan;
            for s=2:ngroup
                group(s)=minquan+double(rem(n,ngroup)>=s-1);
            end
        end
        part=1;
        adjh=floor(group(1)*alfa);
        nsamp=floor(ntrial/ngroup);
        minigr=sum(group);
        obsingroup=fillgroup(n,group,ngroup,seed);
        totgroup=ngroup;
    else
        [part,group,ngroup,adjh,minigr,obsingroup]=deal(0,n,1,h,n,n);
        replow=[50,22,17,15,14,zeros(1,45)];
        if n < replow(p)
            al=1;
            perm=[1:p-1,p-1];                 
            nsamp=nchoosek(n,p);
        else
            al=0;
            nsamp=ntrial;
        end
    end
    
    csteps=csteps1;
    [tottimes,fine,final]=deal(0);
    if part
        bobj1=repmat(inf,ngroup,10); 
        bcoeff1=cell(ngroup,10);
        [bcoeff1{:}]=deal(NaN);
    end
end
bcoeff=cell(1,10);
bobj=repmat(inf,1,10);
[bcoeff{:}]=deal(NaN);
seed=0;
coeffs=repmat(NaN,p,1);
while final ~= 2
    if fine || (~part && final)
        if ~Hsets_ind
            nsamp=10;
        end
        if final
            adjh=h;
            ngroup=1;
            if n*p <= 1e+5
                csteps=csteps3;
            elseif n*p <= 1e+6
                csteps=10-(ceil(n*p/1e+5)-2);
            else
                csteps=1;
            end
            if n > 5000
                nsamp=1;
            end
        else
            adjh=floor(minigr*alfa);
            csteps=csteps2;
        end
    end
    for k=1:ngroup
        for i=1:nsamp
            tottimes=tottimes+1;
            prevobj=0;
            if ~Hsets_ind
                if final
                    if ~isinf(bobj(i))
                        z=bcoeff{i};  	
                    else
                        break
                    end
                elseif fine
                    if ~isinf(bobj1(k,i))
                        z=bcoeff1{k,i};
                    else
                        break
                    end
                else
                    z(1,1)=Inf;
                    while abs(z(1,1)) == Inf
                        if ~part
                            if al
                                k=p;
                                perm(k)=perm(k)+1; 
                                while ~(k==1 || perm(k) <= (n-(p-k)))
                                    k=k-1;
                                    perm(k)=perm(k)+1;
                                    for j=(k+1):p
                                        perm(j)=perm(j-1)+1;
                                    end	
                                end
                                index=perm;
                                if ~isempty(find(perm>n)) %to avoid indexproblems
                                  break
                                end
                            else
                                [index,seed]=randomset(n,p,seed);
                            end
                        else
                            [index,seed]=randomset(group(k),p,seed);
                            index=obsingroup{k}(index);
                        end
                        if p > 1
                            z=x(index,:)\y(index,1);
                            %problems arise whenever the subsample contains
                            %equal x-values. To avoid warnings the tests in line
                            %551 and line 591 are adapted to abs(z(1,1)).
                            %However, in the first run the matrix will still
                            %be singular having coefficients [-inf a b ...] or
                            %[inf inf ...] producing the warning. To avoid
                            %this we turned off the warnings in this
                            %function (line 140).
                        elseif x(index,1) ~= 0
                            z(1,1)=y(index,1)/x(index,1);
                        else
                            z(1,1)=x(index,1);
                        end
                    end
                end
                if abs(z(1,1)) ~= Inf 
                    if ~part || final
                        residu=y-x*z;
                    elseif ~fine
                        residu=y(obsingroup{k},1)-x(obsingroup{k},:)*z;
                    else
                        residu=y(obsingroup{totgroup+1},1)-x(obsingroup{totgroup+1},:)*z;
                    end        
                    more1=0;
                    more2=0;
                    nmore=200;
                    nmore2=nmore/2;
                    if intadj %intercept adjustment
                        [sortres,sortind]=sort(residu);
                        if ~part %n<600 
                            [center,cover,loc]=mcduni(sortres,obsingroup,adjh,obsingroup-adjh+1,alfa);
                            z(p)=z(p)+center;
                            residu=residu-center;
                        elseif ~fine %n>600, first step
                            [center,cover,loc]=mcduni(sortres,size(obsingroup{k},2),adjh,size(obsingroup{k},2)-adjh+1,alfa);   
                            z(p)=z(p)+center;
                            residu=residu-center;
                        elseif ~final && size(obsingroup{totgroup+1},2)-adjh <= nmore %fine =  merged set
                            [center,cover,loc]=mcduni(sortres,size(obsingroup{totgroup+1},2),adjh,size(obsingroup{totgroup+1},2)-adjh+1,alfa);
                            z(p)=z(p)+center;
                            residu=residu-center;
                        elseif final & n-adjh <= nmore %final = complete data set
                            [center,cover,loc]=mcduni(sortres,n,adjh,n-adjh+1,alfa);
                            z(p)=z(p)+center;
                            residu=residu-center;
                        else   
                            [sortres1,sortind1]=sort(abs(sortres));
                            [sortres2,sortind2]=sort(sortres(sortind1(1:adjh)));
                            further = 1;
                            if final && (sortind1(sortind2(1))+nmore-nmore2+adjh-1 > n  || sortind1(sortind2(1))-nmore2< 1)
                                [center,cover,loc]=mcduni(sortres,n,adjh,n-adjh+1,alfa);
                                z(p)=z(p)+center;
                                residu=residu-center;
                            elseif ~final && fine && (sortind1(sortind2(1))+nmore-nmore2+adjh-1 > size(obsingroup{totgroup+1},2) || sortind1(sortind2(1))-nmore2 < 1)
                                [center,cover,loc]=mcduni(sortres,size(obsingroup{totgroup+1},2),adjh,size(obsingroup{totgroup+1},2)-adjh+1,alfa);
                                z(p)=z(p)+center;
                                residu=residu-center;
                            else   
                                while further
                                    sortres2(1:adjh+nmore)=sortres(sortind1(sortind2(1))-nmore2:sortind1(sortind2(1))+adjh-1+nmore-nmore2);
                                    [center,cover,loc]=mcduni(sortres2,adjh+nmore,adjh,nmore+1,alfa);
                                    if loc == 1 && ~more1
                                        if ~more2
                                            nmore=nmore2;
                                            nmore2=nmore2+nmore2;
                                            more1=1;
                                            if sortind1(sortind2(1))-nmore2 < 1
                                                further=0;
                                            end
                                        else 
                                            further=0;
                                        end
                                    else
                                        if loc == nmore+1 && ~more2
                                            if ~more1
                                                nmore=nmore2;
                                                nmore2=-nmore2;
                                                more2=1;
                                                if final && sortind1(sortind2(1))+nmore-nmore2+adjh-1 > n  
                                                    further=0;
                                                elseif fine && (sortind1(sortind2(1))+nmore-nmore2+adjh-1 > size(obsingroup{totgroup+1},2) || sortind1(sortind2(1))-nmore2<1)
                                                    further=0;
                                                end
                                            else 
                                                further = 0;
                                                
                                            end
                                        else
                                            if loc == 1 && more1
                                                if ~more2
                                                    nmore2=nmore2+100;
                                                    if sortind1(sortind2(1))-nmore2 < 1
                                                        further=0;
                                                    end
                                                else 
                                                    further = 0;
                                                end
                                            else
                                                if loc == nmore+1 && more2
                                                    if ~more1
                                                        nmore2=nmore2+100;
                                                        if final && sortind1(sortind2(1))+nmore-nmore2+adjh-1 > n
                                                            further=0;
                                                        elseif fine && (sortind1(sortind2(1))+nmore-nmore2+adjh-1 > size(obsingroup{totgroup+1},2) || sortind1(sortind2(1))-nmore2<1)
                                                            further=0;
                                                        end
                                                    else 
                                                        further=0;
                                                    end
                                                else
                                                    further=0;
                                                end
                                            end
                                        end
                                    end
                                end
                                z(p)=z(p)+center;
                                residu=residu-center;
                            end
                        end
                    end
                end
            end
            for j=1:csteps %csteps on the subsets
                tottimes=tottimes+1;
                if ~Hsets_ind
                    if z(1,1)~=inf
                        [sortres,sortind]=sort(abs(residu));
                        if fine && ~final
                            sortind=obsingroup{totgroup+1}(sortind);                    
                        elseif part && ~final
                            sortind=obsingroup{k}(sortind);
                        end
                        obs_in_set=sort(sortind(1:adjh));
                        teller(obs_in_set) = teller(obs_in_set) + 1;
                        teller(end) = teller(end) + 1;
                    end
                else
                    obs_in_set = Hsets(i,:);
                end
                
                if Hsets_ind ||(~Hsets_ind && z(1,1)~=inf)
                    [Q,R]=qr(x(obs_in_set,:),0);
                    z=R\(Q'*y(obs_in_set,1));
                    if ~part || final
                        residu=y-x*z;   
                    elseif ~fine
                        residu=y(obsingroup{k},1)-x(obsingroup{k},:)*z;
                    else
                        residu=y(obsingroup{totgroup+1},1)-x(obsingroup{totgroup+1},:)*z;
                    end
                    more1=0;
                    more2=0;
                    nmore=200;
                    nmore2=nmore/2;
                    if intadj %intercept adjustment
                        [sortres,sortind]=sort(residu);
                        if ~part
                            [center,cover,loc]=mcduni(sortres,obsingroup,adjh,obsingroup-adjh+1,alfa);
                            z(p)=z(p)+center;
                            residu=residu-center;
                        elseif ~fine
                            [center,cover,loc]=mcduni(sortres,size(obsingroup{k},2),adjh,size(obsingroup{k},2)-adjh+1,alfa);   
                            z(p)=z(p)+center;
                            residu=residu-center;
                        elseif ~final && size(obsingroup{totgroup+1},2)-adjh <= nmore
                            [center,cover,loc]=mcduni(sortres,size(obsingroup{totgroup+1},2),adjh,size(obsingroup{totgroup+1},2)-adjh+1,alfa);
                            z(p)=z(p)+center;
                            residu=residu-center;
                        elseif final && n-adjh <= nmore
                            [center,cover,loc]=mcduni(sortres,n,adjh,n-adjh+1,alfa);
                            z(p)=z(p)+center;
                            residu=residu-center;
                        else   
                            [sortres1,sortind1]=sort(abs(sortres));
                            [sortres2,sortind2]=sort(sortres(sortind1(1:adjh)));
                            further = 1;
                            if final && (sortind1(sortind2(1))+nmore-nmore2+adjh-1 > n  || sortind1(sortind2(1))-nmore2< 1)
                                [center,cover,loc]=mcduni(sortres,n,adjh,n-adjh+1,alfa);
                                z(p)=z(p)+center;
                                residu=residu-center;
                            elseif ~final && fine && (sortind1(sortind2(1))+nmore-nmore2+adjh-1 > size(obsingroup{totgroup+1},2) || sortind1(sortind2(1))-nmore2 < 1)
                                [center,cover,loc]=mcduni(sortres,size(obsingroup{totgroup+1},2),adjh,size(obsingroup{totgroup+1},2)-adjh+1,alfa);
                                z(p)=z(p)+center;
                                residu=residu-center;
                            else   
                                while further
                                    sortres2(1:adjh+nmore)=sortres(sortind1(sortind2(1))-nmore2:sortind1(sortind2(1))+adjh-1+nmore-nmore2);
                                    [center,cover,loc]=mcduni(sortres2,adjh+nmore,adjh,nmore+1,alfa);
                                    if loc == 1 && ~more1
                                        if ~more2
                                            nmore=nmore2;
                                            nmore2=nmore2+nmore2;
                                            more1=1;
                                            if sortind1(sortind2(1))-nmore2 < 1
                                                further=0;
                                            end
                                        else 
                                            further=0;
                                        end
                                    else
                                        if loc == nmore+1 && ~more2
                                            if ~more1
                                                nmore=nmore2;
                                                nmore2=-nmore2;
                                                more2=1;
                                                if final && sortind1(sortind2(1))+nmore-nmore2+adjh-1 > n  
                                                    further=0;
                                                elseif fine && (sortind1(sortind2(1))+nmore-nmore2+adjh-1 > size(obsingroup{totgroup+1},2) | sortind1(sortind2(1))-nmore2<1)
                                                    further=0;
                                                end
                                            else further=0; 
                                            end
                                        else
                                            if loc == 1 && more1
                                                if ~more2
                                                    nmore2=nmore2+100;
                                                    if sortind1(sortind2(1))-nmore2 < 1
                                                        further=0;
                                                    end
                                                else further = 0; 
                                                end
                                            else
                                                if loc == nmore+1 && more2
                                                    if ~more1
                                                        nmore2=nmore2+100;
                                                        if final && sortind1(sortind2(1))+nmore-nmore2+adjh-1 > n
                                                            further=0;
                                                        elseif fine && (sortind1(sortind2(1))+nmore-nmore2+adjh-1 > size(obsingroup{totgroup+1},2) || sortind1(sortind2(1))-nmore2<1)
                                                            further=0;
                                                        end
                                                    else further=0; 
                                                    end
                                                else
                                                    further=0;
                                                end
                                            end
                                        end
                                    end
                                end
                                z(p)=z(p)+center;
                                residu=residu-center;
                            end
                        end
                    end
                    sor=sort(abs(residu));
                    obj=sum(sor(1:adjh).^2); %objective function value after the iteration
                    if j >= 2 && obj == prevobj
                        break;
                    end
                    prevobj=obj;
                end %end final
            end
            if ~final 
                if fine ||~part %merged set or n<600
                    if obj < max(bobj)		
                        [bcoeff,bobj]=insertion(bcoeff,bobj,z,obj,1,eps);
                    end
                else
                    if obj < max(bobj1(k,:))	
                        [bcoeff1,bobj1]=insertion(bcoeff1,bobj1,z,obj,k,eps);
                    end
                end
            end
            if final && obj < bestobj
                bestset=obs_in_set;
                bestobj=obj;
                coeffs=z;
            end 
        end    %end for nsamp
    end %end for ngroups
    if part && ~fine
        fine = 1;
    elseif (part && fine && ~final) || (~part && ~final)
        final = 1;
    else
        final = 2;
    end
end %end while,  so final = 2

if p <= 1
    coeffs(1)=coeffs(1)*datamad(p+1)/datamad(1);
else
    coeffs(1:p-1)=(coeffs(1:p-1)*datamad(p+1))'./datamad(1:p-1);
    if ~intercept
        coeffs(p)=coeffs(p)*datamad(p+1)/datamad(p);
    else
        coeffs(p)=coeffs(p)*datamad(p+1);
        coeffs(p)=coeffs(p)-sum(coeffs(1:p-1)'.*datamed(1:p-1));
        coeffs(p)=coeffs(p)+datamed(p+1);
    end      
end
bestobj=bestobj*(datamad(p+1)^2);
x=xorig;
y=yorig;
raw.coefficients=coeffs;
raw.objective=bestobj;
fitted=x*coeffs;
raw.fitted=fitted;
residuals=y-fitted;
raw.residuals=residuals;
sor=sort(residuals.^2);
factor=rawconsfactorlts(h,n);
sh0=sqrt((1/h)*sum(sor(1:h)));
s0=sh0*factor;
cutoff.resd=sqrt(chi2inv(0.975,1));
if abs(s0) < 1e-7
   weights=abs(residuals)<=1e-7;
   raw.wt=weights;
   raw.scale=0;
   res.scale=0;
   res.coefficients=raw.coefficients;
   raw.objective=0;
else
    raw.scale=s0;
    m=2*n/asvarscalekwad(h,n);
    quantile=tinv(0.9875,m);
    weights=abs(residuals/s0)<=quantile;
    raw.wt=weights;
    [Q,R]=qr(x(weights==1,:),0);
    z=R\(Q'*y(weights==1));
    res.coefficients=z;
    fitted=x*res.coefficients;
    residuals=y-fitted;
    res.scale=sqrt(sum(weights.*residuals.^2)/(sum(weights)-1));
    s0=res.scale;
    weights=abs(residuals/res.scale)<=cutoff.resd;
end
res.flag=weights;

res.Hsubsets.Hopt = bestset;
[telobs,indobs] = greatsort(teller(1:(end - 1)));
res.Hsubsets.Hfreq = indobs(1:(h));
if size(res.Hsubsets.Hfreq,2) == 1
    res.Hsubsets.Hfreq = res.Hsubsets.Hfreq';
end

if intercept
    yw=y(raw.wt==1);
    cyw=mcenter(yw);
    sres=sum(residuals(raw.wt==1).^2);
    cwy2=sum(cyw.^2);
    res.rsquared=1-sres/cwy2;
else
    sor=sort(residuals.^2);
    s1=sum(sor(1:h));
    sor=sort(y.^2);
    sh=sum(sor(1:h));
    res.rsquared=1-(s1/sh);
end
if res.rsquared > 1
   res.rsquared=1;
elseif res.rsquared < 0
    res.rsquared=0;
end

if abs(s0) < 1e-7
    res.method=strvcat(res.method,'An exact fit was found!');
    res.Hsubsets.Hopt=1:n;
    res.Hsubsets.Hfreq=1:n;
    disp('Exact fit was encountered')
end   

if ~intercept
    data=x;
else 
    data=x(:,1:p-1);
end
%calculating residual distances : in case of a univariate analysis they are
%equal to the standardized residuals.
stdres=residuals/res.scale;
cutoff.resd=sqrt(chi2inv(0.975,1));

%assigning ouput
raw=struct('coefficients',{raw.coefficients},'fitted',{raw.fitted},'res',{raw.residuals},'scale',{raw.scale},...
    'objective',{raw.objective},'wt',{raw.wt});
rew=struct('slope',{res.coefficients(1:p)},'int',{0},'fitted',{fitted},'res',{residuals},...
    'scale',{res.scale},'rsquared',{res.rsquared},'h',{h},'Hsubsets',{res.Hsubsets},'alpha',{alfa},'rd',{0},'resd', {stdres},...
    'cutoff',{cutoff},'flag',{res.flag},...
    'method',{res.method},'class',{'LTS'},'classic',{res.classic},'X',{data},'y',{y});

if intercept
    rew=setfield(rew,'int',res.coefficients(p));
    rew=setfield(rew,'slope',res.coefficients(1:p-1));
end

if isfield(rew,'X') && ((size(x,2)-intercept)~=1 || size(y,2)~=1)
    rew=rmfield(rew,{'X','y'});
end
warning on
if plots && classic
    mcdres=mcdcov(data,'h',h,'plots',0);
    if -log(abs(det(mcdres.cov)))/size(data,2)> 50
        res.rd=NaN;
    else
        res.rd=mcdres.rd;
    end
    cutoff.rd=mcdres.cutoff.rd; 
    cutoff.md=mcdres.cutoff.md;
    rew=setfield(rew,'cutoff',cutoff);
    rew=setfield(rew,'rd', res.rd);
    try
        makeplot(rew,'classic',1)
    catch %output must be given even if plots are interrupted 
        %> delete(gcf) to get rid of the menu 
    end
elseif plots
    mcdres=mcdcov(data,'h',h,'plots',0);
    if -log(abs(det(mcdres.cov)))/size(data,2)> 50
        res.rd=NaN;
    else
        res.rd=mcdres.rd;
    end
    cutoff.rd=mcdres.cutoff.rd; 
    cutoff.md=mcdres.cutoff.md;
    rew=setfield(rew,'cutoff',cutoff);
    rew=setfield(rew,'rd', res.rd);
    try
        makeplot(rew)
    catch %output must be given even if plots are interrupted 
        %> delete(gcf) to get rid of the menu 
    end
end
if rew.rd==0
    rew=rmfield(rew,'rd');
end
% --------------------------------------------------------------------

function obsingroup = fillgroup(n,group,ngroup,seed)

% Creates the subdatasets.

obsingroup=cell(1,ngroup+1);
jndex=0;
for k = 1:ngroup
   for m = 1:group(k)
      [random,seed]=uniran(seed);
      ran=floor(random*(n-jndex)+1);
      jndex=jndex+1;
      if jndex == 1
         index(1,jndex)=ran;
         index(2,jndex)=k;
      else
         index(1,jndex)=ran+jndex-1;
         index(2,jndex)=k;
         ii=min(find(index(1,1:jndex-1) > ran-1+[1:jndex-1]));
         if length(ii)
            index(1,jndex:-1:ii+1)=index(1,jndex-1:-1:ii);
            index(2,jndex:-1:ii+1)=index(2,jndex-1:-1:ii);
            index(1,ii)=ran+ii-1;
            index(2,ii)=k;
         end
      end
   end
   obsingroup{k}=index(1,index(2,:)==k);
   obsingroup{ngroup+1}=[obsingroup{ngroup+1},obsingroup{k}];
end

% --------------------------------------------------------------------

function [ranset,seed] = randomset(tot,nel,seed)

for j = 1:nel
   [random,seed]=uniran(seed);
   num=floor(random*tot)+1;
   if j > 1
      while any(ranset==num)
         [random,seed]=uniran(seed);
         num=floor(random*tot)+1;
      end
   end
   ranset(j)=num;   
end


% --------------------------------------------------------------------

function [output] = replow(k,pmax)

replow=[500,50,22,17,15,14];
help=zeros(1,pmax-5);
replow=[replow help];
output=replow(k);

% --------------------------------------------------------------------

function [initmean,initcov,iloc]=mcduni(y,ncas,h,len,alfa)

% The exact MCD algorithm for the univariate case. 

y=sort(y);
ay(1)=sum(y(1:h));
factor=rawconsfactormcd(h,ncas);

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);
ndup=length(ii);
slutn(1:ndup)=ay(ii);
initmean=slutn(floor((ndup+1)/2))/h;
initcov=factor^2*sqmin/h;
iloc=ii(1);

% -----------------------------------------------------------------------------

function [bestmean,bobj]=insertion(bestmean,bobj,z,obj,row,eps)

insert=1;
equ=find(obj==bobj(row,:));
for j=equ
   if (z==bestmean{row,j}) 
      insert=0;
   end
end
if insert 
   ins=min(find(obj < bobj(row,:))); 
   if ins==10
      bestmean{row,ins}=z;
      bobj(row,ins)=obj;
   else
      [bestmean{row,ins+1:10}]=deal(bestmean{row,ins:9});
      bestmean{row,ins}=z;
      bobj(row,ins+1:10)=bobj(row,ins:9);
      bobj(row,ins)=obj;
   end
end

% --------------------------------------------------------------------------------

function quan=quanf(alfa,n,rk)

quan=floor(2*floor((n+rk+1)/2)-n+2*(n-floor((n+rk+1)/2))*alfa);


%---------------------------------------------------------------

function rawconsfacmcd=rawconsfactormcd(quan,n)

qalpha=chi2inv(quan/n,1);
calphainvers=gamcdf(qalpha/2,1/2+1)/(quan/n);
calpha=1/calphainvers;
rawconsfacmcd=calpha;


%-------------------------------------------------------------

function rawconsfaclts=rawconsfactorlts(quan,n)

rawconsfaclts=(1/sqrt(1-((2*n)/(quan*(1/norminv((quan+n)/(2*n)))))*...
		normpdf(1/(1/(norminv((quan+n)/(2*n)))))));

%--------------------------------------------------------------

function asvar=asvarscalekwad(quan,n)

alfa=quan/n;
alfa=1-alfa;
qalfa=chi2inv(1-alfa,1);
c2=gamcdf(qalfa/2,1/2+1);
c1=1/c2;
c3=3*gamcdf(qalfa/2,1/2+2);
asvar=qalfa*(1-alfa)-c2;
asvar=asvar^2;
asvar=(c3-2*qalfa*c2+(1-alfa)*(qalfa^2))-asvar;
asvar=c1^2*asvar;
%--------------------------------------------------------------