function result = rsimca(x,group,varargin)

%RSIMCA performs a robust version of the SIMCA method. This is a classification
% method on a data matrix x with a known group structure. On each group a 
% robust PCA analysis (ROBPCA) is performed. Afterwards a classification
% rule is developped to determine the assignment of new observations. Since RSIMCA
% depends on the ROBPCA method (see robpca.m) it is able to deal with high-dimensional data.
%
% The RSIMCA method is described in:
%    K. Vanden Branden and M. Hubert (2005), 
%    Robust classification in high dimensions based on the SIMCA method,
%    Chemometrics and Intelligent Laboratory Systems, 79, 10-21.
%
% Required input arguments:
%          x : training data set (matrix of size n by p).
%      group : column vector containing the group numbers of the training
%              set x. For the group numbers, any strict positive integer is
%              allowed.
%
% Optional input arguments:
%          alpha : (1-alpha) measures the fraction of outliers (in each group) the algorithm should 
%                  be able to resist. Any value between 0.5 and 1 may be specified. (default = 0.75) 
%                  If k groups are present in the data, per group a different value for alpha may be specified 
%                  (default = alpha = [0.75, ... , 0.75]).
%              k : Is a vector with size equal to the number of groups, or otherwise 0. It represents the number
%                  of components to be retained in each group. (default = 0).
%           kmax : Maximal number of principal components to compute (default = 10).
%                  If k is provided, kmax does not need to be specified, unless k is larger
%                  than 10.    
%          scree : If equal to one, a scree plot is drawn for each group. (default = 0).
%          press : If equal to one, a plot of robust press-values is drawn for each group.
%                  If k is given as input, the default value is 0, else the default value is one.
%         method : Indicates which classification rule is wanted. `1' results in rule (R1)
%                  based on the scaled orthogonal and score distances. `2' corresponds with
%                  (R2) based on the squared scaled orthogonal and score distances. Default is 2. 
%          gamma : Represents the value(s) used in the classification rule: weight gamma is given to the od's,
%                  weight (1-gamma) to the sd's. (default = 0.5).
%     misclassif : String which indicates how to estimate the probability of
%                  misclassification. It can be based on the training data ('training'), 
%                  a validation set ('valid'), or cross-validation ('cv'). Default is 'training'.
% membershipprob : Vector which contains the membership probability of each
%                  group (sorted by increasing group number). These values are used to
%                  obtain the total misclassification percentage. If no priors are given, they are 
%                  estimated as the proportions of regular observations in the training set. 
%          valid : If misclassif was set to 'valid', this field should contain 
%                  the validation set (a matrix of size m by p).
%     groupvalid : If misclassif was set to 'valid', this field should contain the group numbers
%                  of the validation set (a column vector).
%     predictset : Contains a new data set (a matrix of size mp by p) from which the 
%                  class memberships are unknown and should be predicted.  
%          plots : If equal to 1, one figure is created with the training data and the
%                  boundaries for each group. This plot is
%                  only available for trivariate (or smaller) data sets. For technical reasons, a maximum 
%                  of 6 groups is allowed. Default is one.
%    plotsrobpca : If equal to one, a robust score diagnostic plot is
%                  drawn (default). If the input argument 'classic' is equal to one, 
%                  the classical plots are drawn as well.
%                  If 'plots' is equal to zero, this plot is suppressed.
%                  See also makeplot.m
%          labsd : The 'labsd' observations with largest score distance are
%                  labeled on the diagnostic plot. (default = 3)
%          labod : The 'labod' observations with largest orthogonal distance are
%                  labeled on the diagnostic plot. default = 3) 
%            mcd : If equal to one: when the number of variables is sufficiently small,
%                  the loadings are computed as the eigenvectors of the MCD covariance matrix, 
%                  hence the function 'mcdcov.m' is automatically called. The number of 
%                  principal components is then taken as k = rank(x). (default)
%                  If equal to zero, the robpca algorithm is always applied.
%        classic : If equal to one, the classical SIMCA analysis will be performed
%                  (see also csimca.m). (default = 0)
%        compare : If equal to one, the classical SIMCA analysis will be performed
%                  with the same weights and the same priors as the robust analysis
%                  has been performed. This is especially useful to compare the robust
%                  and classical result on the same data with the same priors. (default = 0)
%         
%
% I/O: result=rsimca(x,group,'alpha',0.5,'method',1,'misclassif','training',...
%                  'membershipprob',proportions,'valid',y,'groupvalid',groupy,'plots',0,'classic',0);
%
%  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.
%
% Examples: out=rsimca(x,group,'method','1')
%           out=rsimca(x,group,'plots',0)
%           out=rsimca(x,group,'valid',y,'groupvalid',groupy)
%
% The output is a structure containing the following fields:
%        result.assignedgroup : If there is a validation set, this vector contains the assigned group numbers
%                               for the observations of the validation set. Otherwise it contains the
%                               assigned group numbers of the original observations based on the discriminant rules.
%                  result.pca : A cell containing the results of the different ROBPCA analysis on the training sets.
%               result.method : String containing the method used to obtain
%                               the discriminant rules (either 1 for 'R1' or 2 for 'R2'). This
%                               corresponds to the input argument method. 
%            result.flagtrain : Observations from the training set whose score distance and/or orthogonal distance
%                               exceeds a certain cut-off value can be considered as outliers and receive a flag equal 
%                               to zero. The regular observations receive a flag 1. (See also robpca.m)
%            result.flagvalid : Observations from the validation set whose score distance and/or orthogonal distance 
%                               exceeds a certain cut-off value can be considered as outliers and receive a
%                               flag equal to zero. The regular observations receive a flag 1. 
%                               If there is no validation set, this field is equal to zero.
%         result.grouppredict : If there is a prediction set, this vector contains the assigned group numbers
%                               for the observations of the prediction set. 
%          result.flagpredict : Observations from the new data set (predict) whose robust distance (to the center of their group)
%                               exceeds a certain cut-off value can be considered as overall outliers and receive a
%                               flag equal to zero. The regular observations receive a flag 1. 
%                               If there is no prediction set, this field is
%                               equal to zero.
%       result.membershipprob : A vector with the membership probabilities. If no priors are given, they are estimated 
%                               as the proportions of regular observations in the training set.
%           result.misclassif : String containing the method used to estimate the misclassification probabilities
%                               (same as the input argument misclassif).
%     result.groupmisclasprob : A vector containing the misclassification probabilities for each group.
%       result.avemisclasprob : Overall probability of misclassification (weighted average of the misclassification
%                               probabilities over all groups).
%                result.class : 'RSIMCA'
%              result.classic : If the input argument 'classic' is equal to one, this structure
%                               contains results of the classical SIMCA analysis. 
%              result.compare : If the input argument 'compare' is equal to one, this strucuture 
%                               contains results for the classical SIMCA analysis with the same weights
%                               and priors as in the robust analysis.
%                    result.x : The training data set (same as the input argument x) (only in output when p<=3).
%                result.group : The group numbers of the training set (same as the input argument group) (only in output when p<=3).
%
% This function is part of LIBRA: the Matlab Library for Robust Analysis,
% available at: 
%              http://wis.kuleuven.be/stat/robust.html
%
% Written by Karlien Vanden Branden  
% Last Update: 05/07/2005
%

if nargin<2
    error('There are too few input arguments.')
end

% assigning default-values
[n,p]=size(x);
if size(group,1)~=1
    group=group';
end
if n ~= length(group)
    error('The number of observations is not the same as the length of the group vector!')
end
g=group;
counts=tabulate(g); %contingency table (outputmatrix with 3 colums): value - number - percentage 
[lev,levi,levj]=unique(g);
if ~all(counts(:,2)) %some groups have zero values, omit those groups
    disp(['Warning: group(s) ', num2str(counts(counts(:,2)==0,1)'), 'are empty']);
    empty=counts(counts(:,2)==0,:);
    counts=counts(counts(:,2)~=0,:);
else
    empty=[];
end
ng=size(counts,1);
proportions = zeros(ng,1);
y=0; %initial values of the validation data set and its groupsvector
groupy=0;
labsd = 3;
labod = 3;
plotsrobpca = 0;
press = 1;
scree = 0;
counter=1;
gamma = 0.5;
k = zeros(ng,1);
for iClass = 1:ng
    r(iClass,1) = rank(x(find(group == iClass),:));
    kmax(iClass,1) = min([10,floor(counts(iClass,2)/2),r(iClass)]);
end
default=struct('alpha',0.75,'k',k,'kmax',kmax,'scree',scree,'press',press,'method',2,...
    'gamma',0.5,'misclassif','training','membershipprob',proportions,'valid',y,...
    'groupvalid',groupy,'plots',1,'plotsrobpca',plotsrobpca,'labsd',labsd,'labod',labod,...
    'mcd',0,'classic',0,'compare',0,'predictset',[]);
list=fieldnames(default);
options=default;
IN=length(list);
i=1;
%reading the user's input 
if nargin>2
    %
    %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 
    %
    %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-2 %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
end

%Checking gamma
gamma = options.gamma;
if any(gamma>1) | any(gamma<0)
    error('An inappropriate number for gamma is given. A correct value lies between 0 and 1.');
end

%Checking prior (>0 )
prior=options.membershipprob;
if size(prior,1)~=1
    prior=prior';
end
epsilon=10^-4;
if sum(prior) ~=0 & (any(prior < 0) | (abs(sum(prior)-1)) > epsilon)
   error('Invalid membership probabilities.')
end
if length(prior)~=ng
    error('The number of membership probabilities is not the same as the number of groups.')
end

%%%%%%%%%%%%%%%%%%MAIN FUNCTION %%%%%%%%%%%%%%%%%%%%%
%Checking if a validation set is given
if strmatch(options.misclassif, 'valid','exact') 
    if options.valid==0
        error(['The misclassification error will be estimated through a validation set',...
                'but no validation set is given!'])
    else
        validx = options.valid;
        validgroup = options.groupvalid;
        if size(validx,1)~=length(validgroup)
            error('The number of observations in the validation set is not the same as the length of its group vector!')
        end
        if size(validgroup,1)~=1
            validgroup = validgroup';
        end
        countsvalid=tabulate(validgroup);
        countsvalid=countsvalid(countsvalid(:,2)~=0,:);
        if size(countsvalid,1)==1
            error('The validation set must contain observations from more than one group!')
        elseif any(ismember(empty,countsvalid(:,1)))
            error(['Group(s) ' ,num2str(empty(ismember(empty,countsvalid(:,1)))), 'was/were empty in the original dataset.'])
        end
    end
elseif options.valid~=0
    validx = options.valid;
    validgroup = options.groupvalid;
    if size(validx,1) ~= length(validgroup)
        error('The number of observations in the validation set is not the same as the length of its group vector!')
    end
    if size(validgroup,1)~=1
        validgroup = validgroup';
    end
    options.misclassif='valid';
    countsvalid=tabulate(validgroup);
    countsvalid=countsvalid(countsvalid(:,2)~=0,:);
    if size(countsvalid,1)==1
        error('The validation set must contain more than one group!')
    elseif any(ismember(empty,countsvalid(:,1)))
        error(['Group(s) ' , num2str(empty(ismember(empty,countsvalid(:,1)))), ' was/were empty in the original dataset.'])
    end
end

if length(options.alpha) == 1
    alfa = ones(ng,1)*options.alpha;
elseif length(options.alpha) ~= ng
    error('The length of alpha does not correspond with the number of groups.');
else
    alfa = options.alpha;
end
model.counts = counts(:,2);
model.x = x;
model.group = group;

%Checking input variable k & kmax:
if sum(options.kmax>0)~=ng
    mess=sprintf(['Attention (rsimca.m): The value for kmax is incorrect, kmax = ',num2str(options.kmax'),...
            '\n is smaller than 0. kmax is set to ',num2str(kmax')]);
    disp(mess)
end
if sum(k<=options.kmax) ~= ng
    error('The value for k is set too large.');
end
if sum(options.k ~= 0) == ng 
    press = 0;
    scree = 0;
else
    press = options.press;
    scree = options.scree;
end

labsd = floor(max(0,min(options.labsd,n)));
labod = floor(max(0,min(options.labod,n)));

%RSIMCA:
%PRINCIPAL COMPONENT ANALYSIS
%Perform ROBPCA on each group separately:
%   a) if k is not given: decide on the optimal number of components using CV.
%   b) if k is given: perform robpca with the optimal number of components.

for iClass = 1:ng
    indexgroup = find(g==iClass);
    groupi = x(indexgroup,:); 
    if press == 1 & scree == 1
        disp(['A press curve and a scree plot are drawn for group ',num2str(iClass),'.'])
    elseif press == 1 & scree == 0
        disp(['A press curve is drawn for group ',num2str(iClass),'.'])
    elseif press == 0 & scree == 1
        disp(['A scree plot is drawn for group ',num2str(iClass),'.'])
    end
    model.result{iClass} = robpca(groupi,'k',options.k(iClass),'kmax',options.kmax(iClass),'plots',options.plotsrobpca,...
    'alpha',alfa(iClass),'scree',scree,'press',press,'labsd',labsd,'labod',labod,'mcd',options.mcd);
 
    model.flag(1,indexgroup) = model.result{iClass}.flag.all;
    model.k(iClass) = model.result{iClass}.k;
end

%CLASSIFICATION
%Discriminant rule based on the training set x
[odsc,sdsc] = testmodel(model,x,g);
finalgrouptrain = assigngroup(odsc,sdsc,options.method,ng,gamma); 
result1.weights = model.flag;
if sum(prior) == 0
    for iClass=1:ng
        result1.ingroup(iClass) = sum((group(result1.weights == 1) == repmat(lev(iClass),1,sum(result1.weights))));
    end
    result1.prior = result1.ingroup./sum(result1.ingroup)';
else
    result1.prior = prior;
end

%Compute scaled orthogonal and scaled score distances for the validation set
if strmatch(options.misclassif,'valid','exact') 
    [odsc,sdsc] = testmodel(model,validx);
    finalgroup = assigngroup(odsc,sdsc,options.method,ng,gamma);
elseif strmatch(options.misclassif,'cv','exact')  %use cv
    [odsc,sdsc] = leave1out(model);
    finalgroup = assigngroup(odsc,sdsc,options.method,ng,gamma);
end

switch options.misclassif
case 'valid'
    [v,vi,vj]=unique(validgroup);
    odscgroup = [];
    sdscgroup = [];
    for iClass = 1:ng
        indexgroup = find(validgroup == iClass);
        odscgroup = [odscgroup;odsc(indexgroup,iClass)];
        sdscgroup = [sdscgroup;sdsc(indexgroup,iClass)];
    end    
    weightsvalid=zeros(1,length(odscgroup)); 
    weightsvalid(((odscgroup <= 1) & (sdscgroup <= 1)))=1;
    for igamma = 1:length(gamma)
        for iClass=1:ng 
            misclas(iClass)=sum((validgroup(weightsvalid==1)==finalgroup(weightsvalid==1,igamma)') & ...
                (validgroup(weightsvalid==1)==repmat(v(iClass),1,sum(weightsvalid))));
            ingroup(iClass) = sum((validgroup(weightsvalid == 1) == repmat(v(iClass),1, sum(weightsvalid))));
        end
        misclas = (1 - (misclas./ingroup));
        misclasprobpergroup(igamma,:)=misclas;
        misclas=misclas.*result1.prior;
        misclasprob(igamma)=sum(misclas);    
    end
case 'training'
    for igamma = 1:length(gamma)
        for iClass = 1:ng
            result1.misclas(iClass) = sum((group(result1.weights==1)==finalgrouptrain(result1.weights==1,igamma)')& ...
                (group(result1.weights==1)==repmat(lev(iClass),1,sum(result1.weights))));
            result1.ingroup(iClass) = sum((group(result1.weights == 1) == repmat(lev(iClass),1,sum(result1.weights))));
        end
        misclas = (1 - (result1.misclas./result1.ingroup));
        misclasprobpergroup(igamma,:) = misclas;
        misclas = misclas.*result1.prior;
        misclasprob(igamma) = sum(misclas);
        weightsvalid=0;%only available with validation set
    end
    finalgroup = finalgrouptrain;
case 'cv' 
    for igamma = 1:length(gamma)
        for iClass=1:ng 
            misclas(iClass)=sum((group(result1.weights==1)==finalgroup(result1.weights==1,igamma)') & ...
                (group(result1.weights==1)==repmat(lev(iClass),1,sum(result1.weights))));
            ingroup(iClass) = sum((group(result1.weights == 1) == repmat(lev(iClass),1,sum(result1.weights))));
        end
        misclas = (1 - (misclas./ingroup));
        misclasprobpergroup(igamma,:)=misclas;
        misclas=misclas.*result1.prior;
        misclasprob(igamma)=sum(misclas);    
    end
    weightsvalid=0; %only available with validation set
end

if ~isempty(options.predictset)
    [odscpredict,sdscpredict] = testmodel(model,options.predictset);
    finalgrouppredict = assigngroup(odscpredict,sdscpredict,options.method,ng,gamma)';
    weightspredict = (max((odscpredict <= 1) & (sdscpredict <= 1),[],2))';
else
    finalgrouppredict = 0;
    weightspredict = 0;
end   
    
if options.classic
    classicout=csimca(x,g,'k',model.k,'scree',scree,'press',press,'method',options.method,'gamma',gamma,...
        'misclassif',options.misclassif,'membershipprob',prior,'valid',options.valid,...
        'groupvalid',options.groupvalid,'plots',0,'plotspca',0,'labsd',labsd,'labod',labod,...
        'predictset',options.predictset);
else
     classicout=0;
end

if options.compare
   compareout=csimca(x,g,'k',model.k,'scree',scree,'press',press,'method',options.method,'gamma',gamma,...
        'misclassif',options.misclassif,'membershipprob',result1.prior,'valid',options.valid,...
        'groupvalid',options.groupvalid,'plots',0,'plotspca',0,'labsd',labsd,'labod',labod,...
        'weightstrain',result1.weights,'weightsvalid',weightsvalid,'predictset',options.predictset);
else
    compareout = 0;
end

%Output structure
result = struct('assignedgroup',{finalgroup'},'pca',{model.result},'method',options.method,...
    'flagtrain',{model.flag},'flagvalid',weightsvalid,'grouppredict',finalgrouppredict,'flagpredict',weightspredict,...
    'membershipprob',{result1.prior},'misclassif',{options.misclassif},'groupmisclasprob',{misclasprobpergroup},...
    'avemisclasprob',{misclasprob},'class',{'RSIMCA'},'classic',{classicout},'compare',{compareout},'x',x,'group',group);

if size(x,2)>3
    result=rmfield(result,{'x','group'});
end

%Plots:
try
    if options.plots
        makeplot(result);
    end
catch %output must be given even if plots are interrupted 
    %> delete(gcf) to get rid of the menu 
end

%---------------------------------------------
%Leave-One-Out procedure 

function [odsc,sdsc] = leave1out(model)

nClass = length(model.result);

for iClass = 1:nClass
    index = 1;
    indexgroup = find(model.group == iClass);
    teller_if_lus = 0;
    groupi = model.x(indexgroup,:);
    for i = 1:model.counts(iClass)
        groupia = removal(groupi,i,0);
        GRes = model.result;
        H0 = GRes{iClass}.Hsubsets.H0;
        H1 = GRes{iClass}.Hsubsets.H1;
        Hfreq = GRes{iClass}.Hsubsets.Hfreq;
        Hsets = [H0;H1;Hfreq];
        Hsetsmini = RemoveObsHsets(Hsets,i);
        same.value = 0;
        if isempty(find(H0 == i))
            if teller_if_lus >= 1
                same.value = 1;
            end
            teller_if_lus = teller_if_lus + 1;
        end
        interimRes = removeObsRobpca(groupi,i,GRes{iClass}.k,...
            Hsetsmini,same);
        if isempty(find(H0 == i))
            same.res = interimRes;
        end
        GRes{iClass}.M = interimRes.muk_min_i;
        GRes{iClass}.P = interimRes.Pk_min_i;
        GRes{iClass}.L = interimRes.Lk_min_i;
        GRes{iClass}.T = (groupia - ones(model.counts(iClass)-1,1)*GRes{iClass}.M)*GRes{iClass}.P;
        GRes{iClass}.h = GRes{iClass}.h - 1;
        outDist = CompDist(groupia,GRes{iClass});
        GRes{iClass}.cutoff.od = outDist.cutoff.od;
        GRes{iClass}.cutoff.sd = outDist.cutoff.sd;
       
        %Calculate for each class the sd and the od for the observation that was left out:
        for jClass = 1:nClass
            dataicentered = model.x(indexgroup(index),:)-GRes{jClass}.M;
            scorei = dataicentered*GRes{jClass}.P;
            dataitilde = scorei*GRes{jClass}.P';
            sd(indexgroup(index),jClass) = sqrt(scorei*(diag(1./GRes{jClass}.L))*scorei');
            od(indexgroup(index),jClass) = norm(dataicentered-dataitilde);
            if GRes{jClass}.cutoff.od ~= 0  
                odsc(indexgroup(index),jClass) = od(indexgroup(index),jClass)/GRes{jClass}.cutoff.od;
            else
                odsc(indexgroup(index),jClass) = 0;
            end
            sdsc(indexgroup(index),jClass) = sd(indexgroup(index),jClass)/GRes{jClass}.cutoff.sd;
        end
        index = index + 1;
    end
end

%--------------------------------
function [odsc,sdsc] = testmodel(model,validx,validgroup)

%Apply the given model on the test data to obtain different 
%orthogonal distances and score distances.

nClass = length(model.result);
n = size(validx,1);

for jClass = 1:nClass
    for index = 1:n
        out{jClass} = model.result{jClass};
        dataicentered = validx(index,:)-out{jClass}.M;
        scorei = dataicentered*out{jClass}.P;
        dataitilde = scorei*out{jClass}.P';
        sd(index,jClass) = sqrt(scorei*(diag(1./out{jClass}.L))*scorei');
        od(index,jClass) = norm(dataicentered-dataitilde);
        if out{jClass}.cutoff.od ~= 0  
            odsc(index,jClass) = od(index,jClass)/out{jClass}.cutoff.od;
        else
            odsc(index,jClass) = 0;
        end
        sdsc(index,jClass) = sd(index,jClass)/out{jClass}.cutoff.sd; 
    end
end


%-----------------------------------------------------------------------------------------------------------
function Hsets_min_i = RemoveObsHsets(Hsets,i)

% removes the right index from the $h$-subsets in Hsets to 
% obtain (h - 1)-subsets.
% every h-set is put as a row in Hsets.
% i is the index of the observation that is removed from the whole data.

for r = 1:size(Hsets,1)
    if ~isempty(find(Hsets(r,:)== i))
        Hsets_min_i(r,:) = removal(Hsets(r,:),0,find(Hsets(r,:) == i));
    else
        Hsets_min_i(r,:) = Hsets(r,1:(end-1));
    end

    for j = 1:length(Hsets_min_i(r,:))
        if Hsets_min_i(r,j) > i
            Hsets_min_i(r,j) = Hsets_min_i(r,j) - 1;
        end
    end
end

%-------------------------
function result = assigngroup(odsc,sdsc,method,nClass,gamma);

%Obtain the group assignments for given od's and sd's.

if method == 1 
    sd = sdsc;
    od = odsc;
elseif method == 2 
    sd = sdsc.^2;
    od = odsc.^2;
end

for igamma = 1:length(gamma)
    tdist = gamma(igamma).*od + (1-gamma(igamma)).*sd;
    for i = 1:size(od,1)
        result(i,igamma) = find(tdist(i,:) == min(tdist(i,:)));
    end
end

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

function outDist = CompDist(data,out)

% Calculates the score and orthogonal distances 
% input: data : the original data
%        out is a structure that contains the results of the PCA.


[n,p] = size(data);
r = rank(data);
k = out.k;

% Computing distances 
% Robust score distances 
out.sd=sqrt(libra_mahalanobis(out.T,zeros(size(out.T,2),1),'cov',out.L))';
out.cutoff.sd=sqrt(chi2inv(0.975,out.k));
% Orthogonal distances 
XRc=data-repmat(out.M,n,1);
Xtilde=out.T*out.P';
Rdiff=XRc-Xtilde;
out.od = [];
for i=1:n
    out.od(i,1)=norm(Rdiff(i,:));
end
% Robust cutoff-value for the orthogonal distance
if k~=r
    [m,s]=unimcd(out.od.^(2/3),out.h);
    out.cutoff.od = sqrt(norminv(0.975,m,s).^3); 
else
    out.cutoff.od=0;
end

outDist = out;