cvMcd.m

function result = cvMcd(data,kmax,resMCD,h)

%CVMCD calculates the robust cross-validated PRESS (predicted residual error sum of squares)
% curve for the MCD method in a fast way. 
%
% Input arguments: 
%   data   : the full data set
%   kmax   : the maximal number of components to be considered (mostly kmax = p).
%   resMCD : the result of mcdcov(data,'plots',0,'factor',1)
%   h      : the quantile used in MCD.
%
% output:
%   result.press     : vector of length kmax with the press values 
%   result.weights   : the weights for all observations
%
% This function is part of LIBRA: the Matlab Library for Robust Analysis,
% available at: 
%              http://wis.kuleuven.be/stat/robust.html
%
% Written by Sanne Engelen 
% Last Update: 01/07/2004

% Some initialisations:
n = size(data,1);
p = size(data,2);
r = rank(data);
Pk = [];
Lk = [];
teller_if_lus = 0;

if nargin < 4
    alfa = 0.75;
    h=floor(2*floor((n+p+1)/2)-n+2*(n-floor((n+p+1)/2))*alfa);
end

outWeights = weightscvMcd(data,r,kmax,resMCD,h);

w_min = outWeights.w_min;

Hopt = resMCD.Hsubsets.Hopt;
inputH0.H0 = Hopt;
Tfull = mean(data(Hopt,:));
Sfull = cov(data(Hopt,:));

for i = 1:n 
    % deciding which index should be removed from H0.
    inputH0.same = 0;
    if isempty(find(inputH0.H0 == i))
        inputH0.j = h;
        if teller_if_lus >= 1
            inputH0.same = 1;
        end
        teller_if_lus = teller_if_lus + 1;
    else
        inputH0.j = find(inputH0.H0 == i);
    end
      
    % assigning the input variables:
    inputFull.T = Tfull;
    inputFull.S = Sfull;
    
    if ~inputH0.same
        res = removeObsMcd(data,i,inputH0,inputFull);
    end
    if (isempty(find(inputH0.H0 == i))) & (teller_if_lus == 1)
        resfixed = res;
    end
    if isempty(find(inputH0.H0 == i)) & (teller_if_lus ~= 1)
        res = resfixed;
    end

    P_min_i = res.P_min_i;
    L_min_i = res.L_min_i;
    mu_min_i = res.mu_min_i;
    
     for k = 1:kmax
        clear Pk Lk;        
        Pk = P_min_i(:,1:k);
        Lk = L_min_i(1:k,1:k);
        Xhoedk_min_i(i,(k-1)*p + 1:k*p) = (data(i,:) - mu_min_i)*Pk*Pk' + mu_min_i; 
        
        if k~=r
            odk(i,k) = norm(data(i,:) - Xhoedk_min_i(i,(k-1)*p + 1:k*p));
        else
            odk(i,k) = 0;
        end
    end
end

for k = 1:kmax
    press_min(k) = 1/sum(w_min)*w_min*odk(:,k).^2;
end

result.press = press_min;
result.weights = outWeights;

%----------------------------------------------------------------------------------
function out = weightscvMcd(data,r,kmax,resMCD,h)

% computes the weights used to calculate the robust PRESS values.
% 
% input: 
%   data   : the whole data
%       r  : the rank of the data
%   kmax   : the maximal number of components to be considered
%   resMCD : the result of mcdcov(data,'plots',0,'factor',1)
%   h      : the number of observations on which the computations are based.
%
% output: 
%   out.w_min    : the weights computed by taken the minimum over all k

% Some initialisations:
n = size(data,1);
p = size(data,2);
Pk = [];
Lk = [];
Tik = [];

[P,L] = eig(resMCD.cov);
[L,I] = greatsort(diag(L));
P = P(:,I);

for k = 1:kmax
    Pk = P(:,1:k);
    if h==n
        Lk=L(1:k);
    else
    Lk = chi2inv(h/n,k)/chi2inv(h/n,kmax/2)*L(1:k);% with correction for the factor
    end
    muk = resMCD.center;
    Xhoedk(:,(k-1)*p + 1:k*p) = (data - repmat(muk,n,1))*Pk*Pk' + repmat(muk,n,1); 
    Tk = (data - repmat(muk,n,1))*Pk;
    
    for i =1:n
        % defining the sd for the observation that is left:
        sdk(i,k) = sqrt(mahalanobis(Tk(i,:),zeros(1,k),'cov',diag(Lk)));
    
        % defining the od for the observation that is left:
        if k~=r
            odk(i,k) = norm(data(i,:) - Xhoedk(i,(k-1)*p + 1:k*p));
        else
            odk(i,k) = 1;
        end
    end
    
    % defining weights for odk and sdk:
    if k~=r
        [m,s]=unimcd(odk(:,k).^(2/3),h);
        cutoff(k)=sqrt(norminv(0.975,m,s).^3); 
        wod(:,k) = (odk(:,k) <= cutoff(k));
    else
        cutoff(k)= 0; 
        wod(:,k) = 1;
    end
    wsd(:,k) = (sdk(:,k) <= sqrt(chi2inv(0.975,k)));
end

% determine the weights for every observation:
wk = wsd & wod;

if size(wk,1) == 1 | size(wk,2) == 1
    w_min = wk';
else
    w_min = min(wk');
end

out.w_min = w_min;