function [Bk,intk,sigmayykmaxrew_k,sigmattkmaxrew_k,Hopt] = extractmcdregres(resMcd,T,y,kmax,n,q,k,h,cutoffWeights)
%EXTRACTMCDREGRES is an auxiliary function for cross-validation with RPCR and RSIMPLS.
% It extracts mcd-regression coefficients for a certain k based on
% the raw estimate of the slope and intercept of the regression with kmax coefficients.
% Input arguments:
% resMCD : the results of the mcdregression on the data minus observation i for kmax components.
% kmax : the maximal number of components
% T : the k scores from ROBPCA
% k : the number of components.
%
% Output :
% Bk : the slope for k components.
% intk : the intercept for k components.
% Hopt : remark that only for k < kmax a hopt is provided. For kmax Hopt needs to be calculated outside this function.
% sigmayykmaxrew_k and sigmattkmaxrew_k : estimates of covariance matrices that are needed in other functions.
%
% 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: 08/10/2004
geg = [T,y];
j = k;
mutkmaxraw_k = resMcd.Mu(1:j);
muykmaxraw_k = resMcd.Mu((kmax+1):(kmax + q));
sigmattkmaxraw_k = resMcd.Sigma(1:j,1:j);
sigmaytkmaxraw_k = resMcd.Sigma((kmax+1):(kmax + q),1:j);
sigmatykmaxraw_k = resMcd.Sigma(1:j,(kmax + 1):(kmax + q));
sigmayykmaxraw_k = resMcd.Sigma((kmax + 1):(kmax + q),(kmax + 1):(kmax + q));
% perform csteps:
Sigmaraw_k = [sigmattkmaxraw_k, sigmatykmaxraw_k ; sigmaytkmaxraw_k, sigmayykmaxraw_k];
Muraw_k = [mutkmaxraw_k;muykmaxraw_k]';
prevdet = 0;
if j ~= kmax
geg_k = [T(:,1:j),y];
for noCsteps = 1:10
[sortMD,indsm] = sort(libra_mahalanobis(geg_k,Muraw_k,'cov',Sigmaraw_k));
Sigmaraw_k = cov(geg_k(indsm(1:h),:));
Muraw_k = mean(geg_k(indsm(1:h),:));
obj = det(Sigmaraw_k);
if noCsteps >= 20 |abs(obj - prevdet) < 10^(-4)
if noCsteps >= 20
disp('no convergence in csteps')
end
break
end
prevdet = obj;
end
Hopt = indsm(1:h);
end
mukmax_k = Muraw_k;
mutkmax_k = Muraw_k(1:j)';
muykmax_k = Muraw_k((j+1):(j + q))';
sigmakmax_k = Sigmaraw_k;
sigmattkmax_k = Sigmaraw_k(1:j,1:j);
sigmaytkmax_k = Sigmaraw_k((j+1):(j + q),1:j);
sigmatykmax_k = Sigmaraw_k(1:j,(j + 1):(j + q));
sigmayykmax_k = Sigmaraw_k((j + 1):(j + q),(j + 1):(j + q));
sigmattinvkmax_k = inv(sigmattkmax_k);
Braw_k = sigmattinvkmax_k*sigmatykmax_k;
intraw_k = (muykmax_k - (sigmaytkmax_k*sigmattinvkmax_k*mutkmax_k))';
% reweighting:
rewweights=zeros(n,1);
rewfitted=T(:,1:j)*Braw_k + repmat(intraw_k,n,1);
rewresid=y-rewfitted;
rewE=sigmayykmax_k-Braw_k'*sigmattkmax_k*Braw_k;
for noObs=1:n
if (sqrt(rewresid(noObs,1:q)*inv(rewE)*rewresid(noObs,1:q)')) <= cutoffWeights
rewweights(noObs)=1;
end
end
rewclasscov=cov(geg(rewweights==1,:));
rewclasscenter=mean(geg(rewweights==1,:));
mukmaxrew_k = rewclasscenter;
mutkmaxrew_k = rewclasscenter(1:j)';
muykmaxrew_k = rewclasscenter((j+1):(j + q))';
sigmakmaxrew_k = rewclasscov;
sigmattkmaxrew_k = rewclasscov(1:j,1:j);
sigmaytkmaxrew_k = rewclasscov((j+1):(j + q),1:j);
sigmatykmaxrew_k = rewclasscov(1:j,(j + 1):(j + q));
sigmayykmaxrew_k = rewclasscov((j + 1):(j + q),(j + 1):(j + q));
sigmattinvkmaxrew_k = inv(sigmattkmaxrew_k);
Bk = sigmattinvkmaxrew_k*sigmatykmaxrew_k;
intk = (muykmaxrew_k - (sigmaytkmaxrew_k*sigmattinvkmaxrew_k*mutkmaxrew_k))';%
rewclasscov=cov(geg(rewweights==1,:));
rewclasscenter=mean(geg(rewweights==1,:));
mukmaxrew_k = rewclasscenter;
mutkmaxrew_k = rewclasscenter(1:j)';
muykmaxrew_k = rewclasscenter((j+1):(j + q))';
sigmakmaxrew_k = rewclasscov;
sigmattkmaxrew_k = rewclasscov(1:j,1:j);
sigmaytkmaxrew_k = rewclasscov((j+1):(j + q),1:j);
sigmatykmaxrew_k = rewclasscov(1:j,(j + 1):(j + q));
sigmayykmaxrew_k = rewclasscov((j + 1):(j + q),(j + 1):(j + q));
sigmattinvkmaxrew_k = inv(sigmattkmaxrew_k);
Bk = sigmattinvkmaxrew_k*sigmatykmaxrew_k;
intk = (muykmaxrew_k - (sigmaytkmaxrew_k*sigmattinvkmaxrew_k*mutkmaxrew_k))';