function result = cdq(x,y,c,varargin)
%CDQ computes Censored Depth Quantiles for regression, as described in
%
% Debruyne, M., Hubert, M., Portnoy, S., Vanden Branden, K. (2008),
% "Censored depth quantiles",
% Computational Statistics and Data Analysis, 52, 1604-1614.
%
% Required input arguments:
% x : Data matrix of explanatory variables (also called 'regressors').
% Rows of x represent observations, and columns represent variables.
% If the regression should contain an intercept, the x-matrix should
% contain a column of ones.
% y : A vector with n elements that contains the response variables.
% c : A vector that contains the indices of the censored observations.
%
% Optional input arguments:
% nrCan : Number of candidate hyperplanes to compute objective
% function at first grid point. The default is 500.
% nrCom : The objective function at each grid point after the first one is optimized over
% a set of nrCom hyperplanes having p-1 points in common with
% the previously fitted hyperplane. The default is 100.
% If nrCom = 0, the objective function at each grid point after the first one is optimized over
% a fixed set of 'nrCan' hyperplanes.
% nrLam : Number of hyperplanes used in the computation of the regression depth of each hyperplane.
% The default is 500.
% grid : Vector containing the gridpoints. The default is 0.05:0.05:0.95
% maxit : Maximum number of iterations at each grid point. The default
% is 4.
%
% I/O: result=cdq(x,y,c,'nrCan',nrCan,'nrCom',nrCom,'nrLam',nrLam,'grid',grid,'maxit',maxit);
% 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 of CDQ is a (g x p) matrix containing the regression quantiles,
% with g the number of grid points and p the number of regression parameters.
%
% This function is part of LIBRA: the Matlab Library for Robust Analysis,
% available at:
% http://wis.kuleuven.be/stat/robust.html
%
% Written by Michiel Debruyne
% Last Update: April 27, 2007.
%Handle defaults and optional user inputs.
nrCan=500;
nrCom=100;
nrLam=500;
grid=0.05:0.05:0.95;
maxit=4;
default=struct('nrCan',nrCan,'nrCom',nrCom,'nrLam',nrLam,'grid',grid,'maxit',maxit);
list=fieldnames(default);
options=default;
IN=length(list);
i=1;
%
if nargin==1
error('At least 3 inputs are required, see help cdq.')
elseif nargin==2
error('At least 3 inputs are required, see help cdq.')
end
if nargin>2
nrCan=round(options.nrCan);
nrLam=round(options.nrLam);
nrCom=floor(options.nrCom);
grid=sort(options.grid);
maxit=round(options.maxit);
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('nrCan',varargin{j})
nrCan=varargin{j+1};
if nrCan<1
error('nrCan should be positive, see help cdq.')
end
end
end
if rem(j,2)~=0 % fieldnames are placed on odd index
if strcmp('nrCom',varargin{j})
nrCom=varargin{j+1};
if nrCom<1
error('nrCom should be positive, see help cdq.')
end
end
end
if rem(j,2)~=0 % fieldnames are placed on odd index
if strcmp('nrLam',varargin{j})
nrLam=varargin{j+1};
if nrLam<1
error('nrLam should be positive, see help cdq.')
end
end
end
if rem(j,2)~=0 % fieldnames are placed on odd index
if strcmp('grid',varargin{j})
grid=varargin{j+1};
if size(grid,1)>1 & size(grid,2)>1
error('grid should be a vector containing the grid points, not a matrix. See help cdq.')
elseif size(grid,1)==1 & size(grid,2)==1
error('grid should be a vector containing the grid points, not a number. See help cdq.')
elseif sum(grid>=1)+sum(grid<=0)
error('All entries of grid should be between 0 and 1. See help cdq.')
end
end
end
if rem(j,2)~=0 % fieldnames are placed on odd index
if strcmp('maxit',varargin{j})
maxit=varargin{j+1};
if maxit<2
error('maxit should at least be 2, see help cdq.')
end
end
end
end
end
%Start main program
n=length(y);
p=size(x,2);
L=length(grid);
%Construct sets of hyperplanes.
M=nrCan+nrLam;
for i=1:M
ok=0;
while (ok==0)
a=transpose(willcomb(n,p));
ppu=x(a,:);
if (rank(ppu)==p)
ok=1;
end
end
beti=ppu\y(a);
betakan(:,i)=beti;
end
setCan=betakan(:,1:nrCan);
setLam=betakan(:,nrCan+1 : nrCan+nrLam);
%tauh contains the tau-hats, defined for every crossed censored observation.
tauh=[];
crossedobs=[];
tau=grid(1);
beta(1,1)=tau;
%Determine the first estimate
betaPrev=maxObj(x,y,tau,c,crossedobs,tauh,setCan,setLam);
tauPrev=tau;
beta(2:(p+1),1)=betaPrev;
l=2;
iterations=0;
while l<=L
tau=grid(l);
if nrCom~=0 & iterations==0
setCom=bepKanU(x,y,betaPrev,nrCom);
[betaNext,mm]=maxObj(x,y,tau,c,crossedobs,tauh,setCom,setLam);
elseif nrCom==0
betaNext=maxObj(x,y,tau,c,crossedobs,tauh,setCan,setLam);
else
[betaNext,mm]=maxObj(x,y,tau,c,crossedobs,tauh,setCom,setLam);
end
residuals=y-x*betaNext;
%Good observations were crossed and still are.
resn=residuals(crossedobs)<=0;
goodObs=crossedobs(resn);
%Bad observations were crossed but are not anymore.
badObs=crossedobs(residuals(crossedobs)>0);
%New ones are the ones that are crossed for the first time.
notCross=comp(crossedobs,c);
newObs=notCross(residuals(notCross)<=0);
if iterations<maxit & ((length(badObs)>0)|(length(newObs)>0))
%Retain the weights of the good ones.
tauh=tauh(resn);
crossedobs=goodObs;
%Add the new ones.
nrNew=length(newObs);
crossedobs=[crossedobs newObs];
%Their weight equals the previous grid point.
tauh=[tauh repmat(tauPrev,[1,nrNew])];
iterations=iterations+1;
elseif iterations==maxit
%No solution is found (zero in output).
beta(:,end+1)=[tau; zeros(p,1)];
iterations=0;
residuals=y-x*betaNext;
if ~isempty(tauh)
tauh=tauh(tauh~=tauPrev);
crossedobs=crossedobs(1:length(tauh));
end
l=l+1;
%BetaNext is a good solution.
else
beta(:,end+1)=[tau;betaNext];
iterations=0;
betaPrev=betaNext;
tauPrev=tau;
l=l+1;
end
end
result=beta(2:end,:)';
%function computing the maximum of the objective function
function[beta,maxobjf]=maxObj(x,y,tau,c,crossedobs,tauh,cand,setLam)
n=length(y);
nrCan=size(cand,2);
nrLam=size(setLam,2);
nc=comp(c,1:n);
i=1;
maxobjf=0;
while i<=nrCan
betai=cand(:,i);
ri=y-x*betai;
minobjBetai=10^20;
j=1;
while j<=nrLam & minobjBetai>maxobjf
%Take next direction.
setLamj=betai-setLam(:,j);
if ~isempty(find(abs(setLamj)>0.0000001))
%Project.
proj=x*setLamj;
notCross=comp(crossedobs,c);
Kcomp=[notCross nc];
%Part objective function observations in Kcomp.
riKcompos=ri(Kcomp)>=-0.0000001; projKcompos=proj(Kcomp)>0.0000001;
riKcompne=ri(Kcomp)<=0.0000001; projKcompne=proj(Kcomp)<-0.0000001;
obj1V=tau*sum(riKcompos.*projKcompne)+(1-tau)*sum(riKcompne.*projKcompos);
obj1G=tau*sum(riKcompos.*projKcompos)+(1-tau)*sum(riKcompne.*projKcompne);
if ~isempty(crossedobs)
%Part objective function for pseudo-observations in crossed censured observations.
xcl=x(crossedobs,:)*setLamj; ric=ri(crossedobs); xclp=xcl>0.0000001;
xcln=xcl<-0.0000001; ricn=ric<=0.0000001;ricp=ric>=-0.0000001;
proposResnegV=tauh(xclp & ricn);
pronegResposV=tauh(xcln & ricp);
obj2V=tau*sum((tau-pronegResposV)./(1-pronegResposV))+(1-tau)*sum((tau-proposResnegV)./(1-proposResnegV));
proposResnegG=tauh(xcln & ricn);
pronegResposG=tauh(xclp & ricp);
obj2G=tau*sum((tau-pronegResposG)./(1-pronegResposG))+(1-tau)*sum((tau-proposResnegG)./(1-proposResnegG));
%Part objective function for pseudo-observations at infinity.
pronegV=tauh(xcl<-0.0000001);
obj3V=tau*sum(1-(tau-pronegV)./(1-pronegV));
pronegG=tauh(xcl>0.0000001);
obj3G=tau*sum(1-(tau-pronegG)./(1-pronegG));
else
obj2V=0;
obj2G=0;
obj3V=0;
obj3G=0;
end
%Paste 3 parts.
objG=obj1G+obj2G+obj3G;
objV=obj1V+obj2V+obj3V;
minobjBetai=min([objG,objV,minobjBetai]);
end
j=j+1;
end
if minobjBetai>maxobjf
beta=betai;
maxobjf=minobjBetai;
end
i=i+1;
end
% ------------------------------------------------------------------------
function[thetann] = exchangeob(theta,Minv,zu,yu,zi,yi)
zit=transpose(zi);
u=Minv*zi;
w=-1/(1+zit*u);
thetan=theta-(yi-zit*theta)*w*u;
Mninv=Minv+w*u*transpose(u);
zut=transpose(zu);
un=Mninv*zu;
wn=-1/(1-zut*un);
thetann=thetan+(yu-zut*thetan)*wn*un;
% ------------------------------------------------------------------------
function [nr] = willcomb(n,p)
he=randperm(n);
nr=he(1:p);
% ------------------------------------------------------------------------
function[zComp]=comp(z,c)
n=length(z);
i=1;
while i<=n
c=c(z(i)~=c);
i=i+1;
end
zComp=c;
% ------------------------------------------------------------------------
function[setHyp]=bepKanU(x,y,betaNext,M)
betaNext;
n=size(x,1);
p=size(x,2);
residuals=y-x*betaNext;
a=abs(residuals)<=0.00001;
blq=cumsum(a);
if blq(end)>p
a(blq>p)=0;
end
xa=x(a,:);
ya=y(a,1);
xb=x(a==0,:);
yb=y(a==0,1);
aa=size(xa);
setHyp=zeros(p,M);
setHyp(:,1)=betaNext;
minv=inv(transpose(xa)*xa);
for i=2:M
j=ceil(rand(1,1)*p);
k=ceil(rand(1,1)*(n-p));
bb=size(xb);
aa=size(xa);
if rank([xa(1:(j-1),:);xa((j+1):end,:);xb(k,:)])==p
gre=exchangeob(betaNext,minv,transpose(xa(j,:)),ya(j),transpose(xb(k,:)),yb(k,1));
setHyp(:,i)=gre;
else
setHyp(:,i)=betaNext;
end
end