function result=bagplot(x,varargin)
%BAGPLOT draws a bagplot, which is a generalisation of the univariate boxplot
% to bivariate data. The original bagplot is described in
%
% Rousseeuw, P.J., Ruts, I. and Tukey, J.W. (1999),
% "The bagplot: a bivariate boxplot", The American Statistician, 53, 382-387.
%
% The construction of this bagplot is based on the halfspacedepth (see also halfspacedepth.m).
% As the computation of the halfspacedepth is rather time-consuming, it is recommended at large datasets
% to perform the computations on a random subset of the data. The default size of the subset is 200, but this can be
% modified by the user.
%
% The bagplot can also be computed based on the adjusted outlyingness (see also adjustedoutlyingness.m).
% This method is described in:
%
% Hubert, M., and Van der Veeken, S. (2008),
% "Outlier detection for skewed data", Journal of Chemometrics, to appear.
%
% For the up-to-date references, please consult the website:
% http://wis.kuleuven.be/stat/robust.html
%
% The bagplot based on the adjusted outlyingness can be obtained by setting the optional input argument 'type' equal to 'ao'.
%
% Required input arguments:
% x : bivariate data matrix (observations in the rows, variables in the
% columns)
%
% Optional input arguments:
% sizesubset : When drawing the bagplot based on the halfspacedepth, the size of the subset used to perform
% the main computations.
% type : To draw the bagplot based on the halfspacedepth, this parameter should be equal to 'hd' (default).
% To draw the bagplot based on the adjusted outlyingness, it should be set to 'ao'.
% plots : If equal to 1, a bagplot is drawn (default). If equal to zero, no plot is made.
% colorbag : The color of the bag.
% colorfence : The color of the fence.
% databag : If this parameter is 1, the data within the bag are
% plotted. If this parameter is set equal to 0, data points
% are not displayed.
% datafence : If this parameter is 1, the data within the fence are
% plotted. If this parameter is set equal to 0, data points
% are not displayed.
%
% I/O: result=bagplot(x,'sizesubset',200,'type','hd','colorbag',[],'colorfence',[],'databag',1,'datafence',1);
% 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:
% result=bagplot(x,'colorfence',[0.2 0.2 0.5],'databag',0)
% result=bagplot(x,'datafence',0,colorbag',[1 0 0])
%
%
% The output of bagplot is a structure containing
%
% result.center : center of the data. When 'type=hd', this corresponds with the Tukey median. When
% 'type=ao', the point with smallest adjusted outlyingness
% result.type : same of the input parameter type: 'hd' or 'ao'.
% result.flag : is 0 for outliers and equals 1 for regular points
% result.datatype : is 2 for observations in the bag, 1 for the observations in the
% fence and zero for outliers.
% result.depth : when 'type=hd' this is the halfspacedepth of each data point,
% when 'type=ao' the adjusted outlyingness.
%
% Written by Fabienne Verwerft on 25/05/2005
% Update and revision by Stephan Van der Veeken 10/12/2007
% Last revision: 11/02/2008, 25/02/2008
if nargin<1
error('Input argument is undefined')
end
if size(x,1)<10
error('At least 10 datapoints are needed for the calculation')
end
if size(x,2)~=2
error('Data should be 2 dimensional')
end
S=x;
counter=1;
default=struct('colorbag',[0.6 0.6 1],'colorfence',[0.8 0.8 1],'sizesubset',200,...
'databag',1,'datafence',1,'plots',1,'type','hd');
list=fieldnames(default);
result=default;
IN=length(list);
i=1;
%reading the user's input
if nargin>1
%
%placing inputfields in array of strings
%
for j=1:nargin-1
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 'result'.
%
while counter<=IN
index=strmatch(list(counter,:),chklist,'exact');
if ~isempty(index) %in case of similarity
for j=1:nargin-1 %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
result=setfield(result,chklist{index},varargin{I+1});
index=[];
end
counter=counter+1;
end
end
colorbag=result.colorbag;
colorfence=result.colorfence;
databag=result.databag;
datafence=result.datafence;
plots=result.plots;
type=result.type;
sizesubset=result.sizesubset;
switch type
case 'hd'
s1=S(:,1);
s2=S(:,2);
Q=bagp(s1,s2,sizesubset);
bag=Q.interpol;
datatyp=Q.datatyp;
tukm=Q.tukmed;
depth=Q.depth;
case 'ao'
s1=S(:,1);
s2=S(:,2);
n=length(s1);
Q = adjustedoutlyingness(S);
D=[S,Q.outl,Q.flag,(1:n)'];
P=sortrows(D,3);
L=[P(:,1),P(:,2)];
n=size(P,1);
f=floor(n/2);
g=sum(Q.flag);
h=n-g;
bag=L((1:f),:);
hulp=[ones(f,1);2*ones(g-f,1);3*ones(h,1)];
d1=[P,hulp];
d2=sortrows(d1,5);
datatyp=[d2(:,1),d2(:,2),d2(:,6)];
tukm=L(1,:);
depth=Q.outl;
end
i=find(datatyp(:,3)==1);
data1=datatyp(i,1:2);
i=find(datatyp(:,3)==2);
data2=datatyp(i,1:2);
data=[data1;data2];
i=find(datatyp(:,3)>=3);
outl=datatyp(i,1:2);
plak=[data;bag];
k=convhull(plak(:,1),plak(:,2));
whisk=plak(k,1:2);
if plots==1
fill(whisk(:,1),whisk(:,2),colorfence,'LineStyle','none')
hold on %This is essential in order to prevent that only the last executed plotting command is executed
q=convhull(bag(:,1),bag(:,2));
bagq=bag(q,1:2);
fill(bagq(:,1),bagq(:,2),colorbag)%This line should be placed after the command plot(data...) to only %plot %the data outside the bag
axis square
if databag==1
plot(data1(:,1),data1(:,2),'o','MarkerFaceColor','k','MarkerEdgeColor','k','Markersize',4)
end
if datafence==1
plot(data2(:,1),data2(:,2),'o','MarkerFaceColor','k','MarkerEdgeColor','k','Markersize',4)
end
plot(outl(:,1),outl(:,2),'hk','MarkerFaceColor','k','Markersize',8)
plot(tukm(1),tukm(2),'o','MarkerFaceColor','w','MarkerEdgeColor','w','MarkerSize',10);
plot(tukm(1),tukm(2),'+k','Markersize',8)
switch type
case 'ao'
title('bagplot based on adjusted outlyingness')
set(gcf,'Name', 'Bagplot based on adjusted outlyingness', 'NumberTitle', 'off');
case 'hd'
title('bagplot based on halfspacedepth')
set(gcf,'Name', 'Bagplot based on halfspacedepth', 'NumberTitle', 'off');
end
xmin=min(datatyp(:,1));
xmax=max(datatyp(:,1));
small=0.05*(xmax-xmin);
xaxis1=xmin-small;
xaxis2=xmax+small;
ymin=min(datatyp(:,2));
ymax=max(datatyp(:,2));
small=0.05*(ymax-ymin);
xaxis3=ymin-small;
xaxis4=ymax+small;
axis([xaxis1 xaxis2 xaxis3 xaxis4]);
box on;
hold off
end
Datatyp=datatyp(:,3);
flag=zeros(size(S,1),1);
for i=1:(size(S,1))
if Datatyp(i)==3
flag(i)=0;
else
flag(i)=1;
end
end
datatype=zeros(size(S,1),1);
for i=1:(size(S,1))
if Datatyp(i)==3
datatype(i)=0;
else if Datatyp(i)==2
datatype(i)=1;
else
datatype(i)=2;
end
end
end
if type==1
depth=zeros(size(S,1),1);
for i=1:(size(S,1))
depth(i,1)=halfspacedepth(s1(i,1),s2(i,1),s1,s2);
end
end
result=struct('center',tukm,'depth',depth,'datatype',datatype,'flag',flag);
%-----------
function[pins]=rdraw(n,ntot)
pin=unidrnd(ntot,1,n);
pins=sort(pin);
uu=find(not(diff(pins)==0));
pins=pins(uu);
no=length(pins);
if no<n
i=1;
while i<=(n-no)
in=unidrnd(ntot,1,1);
if sum(pins==in)==0
pins(no+i)=in;
i=i+1;
end
end
end
pins=sort(pins);
%----
function [kount,ADK,empty]=isodepth(x,y,d,varargin)
% ISODEPTH is an algoritm that computes the depth region of a bivariate dataset
% corresponding to depth d.
% First, we have to check whether the data points are in general position. If not,
% a very small random number is added to each of the data points until the
% dataset comes in general position. All this is done in the m-file dithering.
% Then all special k-dividers must be found. The coordinates of the vertices
% of the depth region we are looking for are intersection points of these
% special k-dividers. So, consequently, every intersection point in turn has
% to be tested, for example by computing its depth (see halfspacedepth.m),
% to check whether it is a vertex of the depth region.
%
% The ISODEPTH algoritm is described in:
% Ruts, I., Rousseeuw, P.J. (1996),
% "Computing depth contours of bivariate point clouds",
% Computational Statistics and Data Analysis, 23, 153-168.
%
% Required input arguments:
% x : vector containing the first coordinates of all the data
% points
% y : vector containing the second coordinates of all the data
% points
% d : the depth of which the depth region has to be constructed
%
%
% I/O: [kount, ADK, empty]= isodepth(x,y,d);
%
% The output of isodepth is given by
%
% kount : the total number of vertices of the depth region
% ADK : the coordinates of the vertices of the depth region
% empty : logical value (1 if the depth region is empty, 0 if not)
%
% This function is part of the Matlab Library for Robust Analysis,
% available at:
% http://wis.kuleuven.be/stat/robust.html
%
% Last Update: 29/04/2005
n=length(x);
eps=0.0000001;
%
% Checking input
%
if length(x)==1
error('x is not a vector')
elseif not(length(x)==length(y))
error('The vectors x and y must have the same length.')
end
if sum(isnan(x))>=1 || sum(isnan(y))>=1
error('Missing values are not allowed')
end
if sum(x==x(1))==n
error('All data points ly on a vertical line.')
elseif sum(y==y(1))==n
error('All data points ly on a horizontal line.')
else
R=corrcoef(x,y);
if abs(R(1,2))==1
error('All data points are collineair.')
end
end
%
% Check whether the data is in general position. If not, add a very small random
% number to each of the data points.
%
[x,y, Index, angl, ind1,ind2]=dithering(x,y);
%
% Main part
%
if (n==1)&&(d==1)
kount=n;
ADK=[x,y];
empty=0;
return
end
%
if (d>floor(n/2))
kount=0;
ADK=0;
empty=1;
return
end
%
if n<=3
kount=n;
ADK=[x,y];
empty=1;
return
end
%
nrank(Index)=(1:n);
%
% Let the line rotate from zero to angle(1)
%
ncirq=Index;
kount=1;
halt=0;
M=length(angl);
if angl(1)>(pi/2)
L=1;
D1=ind1(L);
IV1=nrank(D1);
D2=ind2(L);
IV2=nrank(D2);
IV=ncirq(IV1);
ncirq(IV1)=ncirq(IV2);
ncirq(IV2)=IV;
IV=IV1;
nrank(D1)=IV2;
nrank(D2)=IV;
%
if ((IV1==d) && (IV2==(d+1)))||((IV2==d) && (IV1==(d+1)))||((IV1==(n-d)) && (IV2==(n-d+1)))||((IV2==(n-d)) && (IV1==(n-d+1)))
if angl(L)<(pi/2)
dum=angl(L)+(pi/2);
else
dum=angl(L)-(pi/2);
end
if (IV1==d && IV2==(d+1))||(IV2==d && IV1==(d+1))
if dum<=(pi/2)
alfa(kount)=angl(L)+pi;
else
alfa(kount)=angl(L);
end
end
if or((IV1==(n-d) && IV2==(n-d+1)),(IV2==(n-d) && IV1==(n-d+1)))
if dum<=(pi/2)
alfa(kount)=angl(L);
else
alfa(kount)=angl(L)+pi;
end
end
kand1(kount)=ind1(L);
kand2(kount)=ind2(L);
D(kount)=sin(alfa(kount))*x(ind1(L))-cos(alfa(kount))*y(ind1(L));
kount=kount+1;
end
halt=1;
end
%
L=2;
stay=1;
% jflag keeps track of which angle we have to test next
while stay==1
stay=0;
kontrol=0;
if (pi<=(angl(L)+(pi/2))) && ((angl(L)-(pi/2))< angl(1))
D1=ind1(L);
IV1=nrank(D1);
D2=ind2(L);
IV2=nrank(D2);
IV=ncirq(IV1);
ncirq(IV1)=ncirq(IV2);
ncirq(IV2)=IV;
IV=IV1;
nrank(D1)=IV2;
nrank(D2)=IV;
%
if ((IV1==d) && (IV2==(d+1)))||((IV2==d) && (IV1==(d+1)))||((IV1==(n-d)) && (IV2==(n-d+1)))||((IV2==(n-d)) && (IV1==(n-d+1)))
if angl(L)<(pi/2)
dum=angl(L)+(pi/2);
else
dum=angl(L)-(pi/2);
end
if (IV1==d && IV2==(d+1))||(IV2==d && IV1==(d+1))
if dum<=(pi/2)
alfa(kount)=angl(L)+pi;
else
alfa(kount)=angl(L);
end
end
if or((IV1==(n-d) && IV2==(n-d+1)),(IV2==(n-d) && IV1==(n-d+1)))
if dum<=(pi/2)
alfa(kount)=angl(L);
else
alfa(kount)=angl(L)+pi;
end
end
kand1(kount)=ind1(L);
kand2(kount)=ind2(L);
D(kount)=sin(alfa(kount))*x(ind1(L))-cos(alfa(kount))*y(ind1(L));
kount=kount+1;
end
kontrol=1;
end
L=L+1;
if kontrol==1
halt=1;
end
if (L==(M+1)) && (kontrol==1)
jflag=1;
stay=2;
end
if not(stay==2)
if ((halt==1)&&(kontrol==0))||(L==(M+1))
stay=3;
else
stay=1;
end
end
end
if not(stay==2)
if (L>1)
jflag=L-1;
else
jflag=M;
end
end
%
halt2=0;
if not(stay==2)
J=0;
%
% If the first switch didnt occur between 0 and the angle angl(1) look for it
% between the following angles.
%
stay2=1;
if (L==M+1) && (kontrol==0)
halt=0;
halt2=0;
J=J+1;
if J==(M+1)
J=1;
end
L=J+1;
if L==(M+1)
L=1;
end
while stay2==1
stay2=0;
kontrol=0;
if (angl(L)+pi/2)<pi
ang1=angl(L)+pi/2;
else
ang1=angl(L)-pi/2;
end
if J==M
jj=1;
if halt2==0
angl(1)=angl(1)+pi;
end
else
jj=J+1;
end
if (angl(J)<=ang1) && (ang1<angl(jj))
if angl(1)>pi
angl(1)=angl(1)-pi;
end
D1=ind1(L);
IV1=nrank(D1);
D2=ind2(L);
IV2=nrank(D2);
IV=ncirq(IV1);
ncirq(IV1)=ncirq(IV2);
ncirq(IV2)=IV;
IV=IV1;
nrank(D1)=IV2;
nrank(D2)=IV;
%
if ((IV1==d) && (IV2==(d+1)))||((IV2==d) && (IV1==(d+1)))||((IV1==(n-d)) && (IV2==(n-d+1)))||((IV2==(n-d)) && (IV1==(n-d+1)))
if angl(L)<(pi/2)
dum=angl(L)+(pi/2);
else
dum=angl(L)-(pi/2);
end
if (IV1==d && IV2==(d+1))||(IV2==d && IV1==(d+1))
if dum<=(pi/2)
alfa(kount)=angl(L)+pi;
else
alfa(kount)=angl(L);
end
end
if or((IV1==(n-d) && IV2==(n-d+1)),(IV2==(n-d) && IV1==(n-d+1)))
if dum<=(pi/2)
alfa(kount)=angl(L);
else
alfa(kount)=angl(L)+pi;
end
end
kand1(kount)=ind1(L);
kand2(kount)=ind2(L);
D(kount)=sin(alfa(kount))*x(ind1(L))-cos(alfa(kount))*y(ind1(L));
kount=kount+1;
end
kontrol=1;
end
if angl(1)>pi
angl(1)=angl(1)-pi;
end
if L==M
L=1;
else
L=L+1;
end
if kontrol==1
halt=1;
end
if (halt==1)&&(kontrol==0)
if halt2==1
stay2=2;
end
if not(stay2==2)
if L>1
jflag=L-1;
else
jflag=M;
end
stay2=0;
end
else
if L==jj
if jj==1
halt2=1;
end
J=J+1;
if J==(M+1)
J=1;
end
L=J+1;
if L==(M+1)
L=1;
end
stay2=1;
else
stay2=1;
end
end
end
end
end
%
if not(stay2==2)
%
% The first switch has occurred. Now start looking for the next ones,
% between the following angles.
%
for i=(J+1):(M-1)
L=jflag;
stay=1;
while stay==1
stay=0;
kontrol=0;
if ((angl(L)+pi/2)<pi)
ang1=angl(L)+pi/2;
else
ang1=angl(L)-pi/2;
end
if (angl(i)<=ang1)&&(ang1<angl(i+1))
D1=ind1(L);
IV1=nrank(D1);
D2=ind2(L);
IV2=nrank(D2);
IV=ncirq(IV1);
ncirq(IV1)=ncirq(IV2);
ncirq(IV2)=IV;
IV=IV1;
nrank(D1)=IV2;
nrank(D2)=IV;
%
if ((IV1==d) && (IV2==(d+1)))||((IV2==d) && (IV1==(d+1)))||((IV1==(n-d)) && (IV2==(n-d+1)))||((IV2==(n-d)) && (IV1==(n-d+1)))
if angl(L)<(pi/2)
dum=angl(L)+(pi/2);
else
dum=angl(L)-(pi/2);
end
if (IV1==d && IV2==(d+1))||(IV2==d && IV1==(d+1))
if dum<=(pi/2)
alfa(kount)=angl(L)+pi;
else
alfa(kount)=angl(L);
end
end
if or((IV1==(n-d) && IV2==(n-d+1)),(IV2==(n-d) && IV1==(n-d+1)))
if dum<=(pi/2)
alfa(kount)=angl(L);
else
alfa(kount)=angl(L)+pi;
end
end
kand1(kount)=ind1(L);
kand2(kount)=ind2(L);
D(kount)=sin(alfa(kount))*x(ind1(L))-cos(alfa(kount))*y(ind1(L));
kount=kount+1;
end
kontrol=1;
end
if kontrol==0
jflag=L;
else
if not(L==M)
L=L+1;
else
L=1;
end
stay=1;
end
end
end
L=jflag;
%
% Finally, look for necessary switches between the last angle and zero.
%
stay=1;
while stay==1
kontrol=0;
stay=0;
if (angl(L)+pi/2)<pi
ang1=angl(L)+pi/2;
else
ang1=angl(L)-pi/2;
end
if (angl(M)<=ang1)&&(ang1<pi)
D1=ind1(L);
IV1=nrank(D1);
D2=ind2(L);
IV2=nrank(D2);
IV=ncirq(IV1);
ncirq(IV1)=ncirq(IV2);
ncirq(IV2)=IV;
IV=IV1;
nrank(D1)=IV2;
nrank(D2)=IV;
%
if ((IV1==d) && (IV2==(d+1)))||((IV2==d) && (IV1==(d+1)))||((IV1==(n-d)) && (IV2==(n-d+1)))||((IV2==(n-d)) && (IV1==(n-d+1)))
if angl(L)<(pi/2)
dum=angl(L)+(pi/2);
else
dum=angl(L)-(pi/2);
end
if (IV1==d && IV2==(d+1))||(IV2==d && IV1==(d+1))
if dum<=(pi/2)
alfa(kount)=angl(L)+pi;
else
alfa(kount)=angl(L);
end
end
if or((IV1==(n-d) && IV2==(n-d+1)),(IV2==(n-d) && IV1==(n-d+1)))
if dum<=(pi/2)
alfa(kount)=angl(L);
else
alfa(kount)=angl(L)+pi;
end
end
kand1(kount)=ind1(L);
kand2(kount)=ind2(L);
D(kount)=sin(alfa(kount))*x(ind1(L))-cos(alfa(kount))*y(ind1(L));
kount=kount+1;
end
kontrol=1;
end
if kontrol==1
if not(L==M)
L=L+1;
else
L=1;
end
stay=1;
end
end
end
num=kount-1; % num is the total number of special k-dividers
%
% Sort the num special k-dividers. Permute kand1, kand2 and D in the same
% way.
%
[alfa,In]=sort(alfa);
kand1=kand1(In);
kand2=kand2(In);
D=D(In);
%
IW1=1;
IW2=2;
Jfull=0;
NDK=0;
stay2=1;
while stay2==1
stay2=0;
ndata=0;
%
% Compute the intersection point.
%
while abs(-sin(alfa(IW2))*cos(alfa(IW1))+sin(alfa(IW1))*cos(alfa(IW2)))<eps
IW2=IW2+1;
ndata=0;
if IW2==(num+1)
IW2=1;
end
end
%
xcord=(cos(alfa(IW2))*D(IW1)-cos(alfa(IW1))*D(IW2))/(-sin(alfa(IW2))*cos(alfa(IW1))+sin(alfa(IW1))*cos(alfa(IW2)));
ycord=(-sin(alfa(IW2))*D(IW1)+sin(alfa(IW1))*D(IW2))/(-sin(alfa(IW1))*cos(alfa(IW2))+sin(alfa(IW2))*cos(alfa(IW1)));
%
% Test whether the intersection point is a data point. If so,
% adjust IW1 and IW2.
%
if or(kand1(IW1)==kand1(IW2),kand1(IW1)==kand2(IW2))
ndata=kand1(IW1);
end
if or(kand2(IW1)==kand1(IW2),kand2(IW1)==kand2(IW2))
ndata=kand2(IW1);
end
if not(ndata==0)
iv=0;
stay=1;
while stay==1
stay=0;
next=IW2+1;
iv=iv+1;
if next==(num+1)
next=1;
end
if not(next==IW1)
if or(ndata==kand1(next),ndata==kand2(next))
IW2=IW2+1;
if (IW2==(num+1))
IW2=1;
end
stay=1;
end
end
end
if iv==(num-1)
kount=1;
ADK=[x(ndata),y(ndata)];
empty=0;
return
end
end
if IW2==num
kon=1;
else
kon=IW2+1;
end
if kon==IW1
kon=kon+1;
end
if kon==(num+1)
kon=1;
end
%
% Test whether the intersection point lies to the left of the special
% k-divider which corresponds to alfa(kon). If so, compute its depth.
%
stay3=1;
stay4=1;
if (sin(alfa(kon))*xcord-cos(alfa(kon))*ycord-D(kon))<=eps
hdep1=halfspacedepth(xcord,ycord,x,y);
if hdep1==d
NDK=1;
else
hdep2=halfspacedepth(xcord-0.000001,ycord-0.000001,x,y);
hdep3=halfspacedepth(xcord+0.000001,ycord+0.000001,x,y);
hdep4=halfspacedepth(xcord-0.000001,ycord+0.000001,x,y);
hdep5=halfspacedepth(xcord+0.000001,ycord-0.000001,x,y);
hdepvector=[hdep1;hdep2;hdep3;hdep4;hdep5];
if (NDK==0)&&(sum(hdepvector>=d)>=1)
NDK=1;
end
if (hdep1<d)&&(hdep2<d)&&(hdep3<d)&&(hdep4<d)&&(hdep5<d)&&(NDK==1)
%
% The intersection point is not the correct one, try the next
% special k-divider.
%
IW2=IW2+1;
if IW2==(num+1)
IW2=1;
end
stay2=1;
end
end
if not(stay2==1)
%
% Store IW2 and IW1 in kornr. If kornr has already been filled,
% check wether we have encountered this intersection point before.
%
if (IW2>IW1)&&(Jfull==0)
kornr(IW1:(IW2-1),1)=kand1(IW1);
kornr(IW1:(IW2-1),2)=kand2(IW1);
kornr(IW1:(IW2-1),3)=kand1(IW2);
kornr(IW1:(IW2-1),4)=kand2(IW2);
else
if IW2>IW1
i=IW1;
stay3=1;
while stay3==1
if (kornr(i,1)==kand1(IW1))&&(kornr(i,2)==kand2(IW1))&&(kornr(i,3)==kand1(IW2))&&(kornr(i,4)==kand2(IW2))
stay3=0;
else
m1=(y(kornr(i,2))-y(kornr(i,1)))/(x(kornr(i,2))-x(kornr(i,1)));
m2=(y(kornr(i,4))-y(kornr(i,3)))/(x(kornr(i,4))-x(kornr(i,3)));
if not(m1==m2)
xcord1=(m1*x(kornr(i,1))-y(kornr(i,1))-m2*x(kornr(i,3))-y(kornr(i,3)))/(m1-m2);
ycord1=(m2*(m1*x(kornr(i,1))-y(kornr(i,1)))-m1*(m2*x(kornr(i,3))-y(kornr(i,3))))/(m1-m2);
end
if (abs(xcord1-xcord)<eps)&&(abs(ycord1-ycord)<eps)
stay3=0;
end
if stay3==1
kornr(i,1)=kand1(IW1);
kornr(i,2)=kand2(IW1);
kornr(i,3)=kand1(IW2);
kornr(i,4)=kand2(IW2);
end
end
if stay3==1
i=i+1;
if i==IW2
stay3=2;
i=i-1;
end
end
end
else
Jfull=1;
kornr(IW1:num,1)=kand1(IW1);
kornr(IW1:num,2)=kand2(IW1);
kornr(IW1:num,3)=kand1(IW2);
kornr(IW1:num,4)=kand2(IW2);
i=1;
stay4=1;
if IW2==1
stay4=3;
end
while stay4==1
if (kornr(i,1)==kand1(IW1))&&(kornr(i,2)==kand2(IW1))&&(kornr(i,3)==kand1(IW2))&&(kornr(i,4)==kand2(IW2))
stay4=0;
else
m1=(y(kornr(i,2))-y(kornr(i,1)))/(x(kornr(i,2))-x(kornr(i,1)));
warning off MATLAB:divideByZero
m2=(y(kornr(i,4))-y(kornr(i,3)))/(x(kornr(i,4))-x(kornr(i,3)));
warning off MATLAB:divideByZero
if not(m1==m2)
xcord1=(m1*x(kornr(i,1))-y(kornr(i,1))-m2*x(kornr(i,3))-y(kornr(i,3)))/(m1-m2);
ycord1=(m2*(m1*x(kornr(i,1))-y(kornr(i,1)))-m1*(m2*x(kornr(i,3))-y(kornr(i,3))))/(m1-m2);
end
if (abs(xcord1-xcord)<=eps)&&(abs(ycord1-ycord)<=eps)
stay4=0;
end
if stay4==1
kornr(i,1)=kand1(IW1);
kornr(i,2)=kand2(IW1);
kornr(i,3)=kand1(IW2);
kornr(i,4)=kand2(IW2);
end
end
if stay4==1
i=i+1;
if i==IW2
i=i-1;
stay4=2;
end
end
end
end
end
end
elseif (stay3>0)&&(stay4>0)&¬(stay2==1)
%
% The intersection point is not the correct one, try the next
% special k-divider.
%
IW2=IW2+1;
if IW2==(num+1)
IW2=1;
end
stay2=1;
end
%
% Look for the next vertex of the convex figure.
%
if (stay3>0)&&(stay4>0)&¬(stay2==1)
IW1=IW2;
IW2=IW2+1;
if IW2==(num+1)
IW2=1;
end
stay2=1;
end
end
%
% Scan kornr and ascribe the coordinates of the vertices to the variable
% ADK.
%
kount=0;
%
if NDK==0
%
% The requested depth region is empty
%
ADK=0;
empty=1;
return
else
empty=0;
end
%
i=1;
E1=y(kornr(i,2))-y(kornr(i,1));
F1=x(kornr(i,1))-x(kornr(i,2));
G1=x(kornr(i,1))*(y(kornr(i,2))-y(kornr(i,1)))-y(kornr(i,1))*(x(kornr(i,2))-x(kornr(i,1)));
E2=y(kornr(i,4))-y(kornr(i,3));
F2=x(kornr(i,3))-x(kornr(i,4));
G2=x(kornr(i,3))*(y(kornr(i,4))-y(kornr(i,3)))-y(kornr(i,3))*(x(kornr(i,4))-x(kornr(i,3)));
xcord(i)=(-F2*G1+F1*G2)/(E2*F1-E1*F2);
ycord(i)=(-E2*G1+E1*G2)/(E1*F2-E2*F1);
DK(i,:)=[xcord(i),ycord(i)];
Juisteind(i)=i;
xcord1=xcord(i);
ycord1=ycord(i);
xcordp=xcord(i);
ycordp=ycord(i);
kount=kount+1;
i=i+1;
%
while not(i==num+1)
if (kornr(i,1)==kornr(i-1,1))&&(kornr(i,2)==kornr(i-1,2))&&(kornr(i,3)==kornr(i-1,3))&&(kornr(i,4)==kornr(i-1,4))
i=i+1;
else
if (kornr(i,1)==kornr(1,1))&&(kornr(i,2)==kornr(1,2))&&(kornr(i,3)==kornr(1,3))&&(kornr(i,4)==kornr(1,4))
pp=find(not(Juisteind==0));
ADK=DK(pp,:);
empty=0;
return
else
E1=y(kornr(i,2))-y(kornr(i,1));
F1=x(kornr(i,1))-x(kornr(i,2));
G1=x(kornr(i,1))*(y(kornr(i,2))-y(kornr(i,1)))-y(kornr(i,1))*(x(kornr(i,2))-x(kornr(i,1)));
E2=y(kornr(i,4))-y(kornr(i,3));
F2=x(kornr(i,3))-x(kornr(i,4));
G2=x(kornr(i,3))*(y(kornr(i,4))-y(kornr(i,3)))-y(kornr(i,3))*(x(kornr(i,4))-x(kornr(i,3)));
xcord(i)=(-F2*G1+F1*G2)/(E2*F1-E1*F2);
ycord(i)=(-E2*G1+E1*G2)/(E1*F2-E2*F1);
if ((abs(xcord(i)-xcordp)<eps)&&(abs(ycord(i)-ycordp)<eps))||((abs(xcord(i)-xcord1)<eps)&&(abs(ycord(i)-ycord1)<eps))
i=i+1;
else
xcordp=xcord(i);
ycordp=ycord(i);
DK(i,:)=[xcord(i),ycord(i)];
Juisteind(i)=i;
kount=kount+1;
i=i+1;
end
end
end
end
%
% Delete all the empty spaces in the matrix DK. The result is ADK.
%
pp=find(not(Juisteind==0));
ADK=DK(pp,:);
%-------
function [tukmed]=halfmed(x,y,varargin)
% HALFMED is an algoritm that computes the Tukey median of a two-dimensional
% dataset. First, we have to check whether the data points are in general position.
% If not, a random number is added to each of the data points until the dataset
% comes in general position. All this is done in the m-file dithering. Then
% the deepest depth region is constructed. One can accelerate this search
% by giving the optional argument kstar, which is usually the maximal halfspace
% depth of the data points. This advantage is used in the functie bagp.m to
% quicken the computations. The Tukey median is the center of gravity of
% this deepest depth region.
%
% The HALFMED algoritm is described in:
% Rousseeuw, P.J., Ruts, I. (1998),
% "Constructing the bivariate Tukey median",
% Statistica Sinica, 8, 827-839.
%
% Required input arguments:
% x : vector containing the first coordinates of all the data
% points
% y : vector containing the second coordinates of all the data
% points
%
% Optional input argument:
% kstar : an integer between ceil(n/3) and floor(n/2)
% One can also use the maximum halfspace depth of the data
% points. (default = 0)
%
%
% I/O: result=halfmed(x,y,'kstar',0);
% The name of the input arguments needs to be followed by their value.
%
% The output of halfmed is given by
%
% result : the coordinates of the Tukey median
%
% This function is part of the Matlab Library for Robust Analysis,
% available at:
% http://wis.kuleuven.be/stat/robust.html
%
% Last Update: 29/04/2005
xsum=0;
ysum=0;
tukmed=0;
n=length(x);
%
% Checking input
%
if length(x)==1
error('x is not a vector')
elseif not(length(x)==length(y))
error('The vectors x and y must have the same length.')
end
if sum(isnan(x))>=1 || sum(isnan(y))>=1
error('Missing values are not allowed')
end
if sum(x==x(1))==n
error('All data points ly on a vertical line.')
elseif sum(y==y(1))==n
error('All data points ly on a horizontal line.')
else
R=corrcoef(x,y);
if abs(R(1,2))==1
error('All data points are collineair.')
end
end
%
% Looking for optional argument
%
if nargin>2
if strcmp(varargin{1},'kstar')
kstar=varargin{2};
if length(kstar)>1
error('kstar must be a number, not a vector.')
end
else
error('Only kstar can be provided as an optional argument')
end
else
kstar=0;
end
if kstar>=floor(n/2)
error('kstar must be smaller than floor(n/2) because the depth region corresponding with kstar is empty')
end
%
% Check whether the data are in general position. If not, add a very small random
% number to each of the data points.
%
[x,y, Index, angl, ind1,ind2]=dithering(x,y);
%
%Calculation of the Tukey median
%
if n<=3
xsum=sum(x);
ysum=sum(y);
tukmed=[xsum/n,ysum/n];
return
end
%
if kstar==0
ib=ceil(n/3);
else
ib=kstar;
end
ie=floor(n/2);
stay=1;
while stay==1
le=ie-ib;
if le<0
le=0;
end
if le==0
stay=0;
end
if stay==1
[kou,dk,empty]=isodepth(x,y,ib+ceil(le/2));
if empty==1
ie=ib+ceil(le/2);
end
if empty==0
ib=ib+ceil(le/2);
end
if le==1
stay=0;
end
end
end
[kount,DK,empty]=isodepth(x,y,ib);
xsum=sum(DK(:,1));
if not(DK==0)
ysum=sum(DK(:,2));
else
ysum=0;
end
wx=DK(:,1);
if not(DK==0)
wy=DK(:,2);
else
wy=0;
end
%
% The maximal depth is now ib.
%
%
% Calculation of the center of gravity
%
som=0;
tukmed=0;
if kount>1
wx=wx-(xsum/kount);
wy=wy-(ysum/kount);
for i=1:(kount-1)
som=som+abs(wx(i)*wy(i+1)-wx(i+1)*wy(i));
tukmed=tukmed+[(wx(i)+wx(i+1))*abs(wx(i)*wy(i+1)-wx(i+1)*wy(i)),(wy(i)+wy(i+1))*abs(wx(i)*wy(i+1)-wx(i+1)*wy(i))];
end
som=som+abs(wx(kount)*wy(1)-wx(1)*wy(kount));
tukmed=tukmed+[(wx(kount)+wx(1))*abs(wx(kount)*wy(1)-wx(1)*wy(kount)),(wy(kount)+wy(1))*abs(wx(kount)*wy(1)-wx(1)*wy(kount))];
tukmed=(tukmed/(3*som))+[(xsum/kount),(ysum/kount)];
else
tukmed=[xsum,ysum];
end
%--------------
function result=bagp(x,y,sizesubset)
% BAGP is an algoritm that makes the necessary computations to construct
% the bagplot of a bivariate dataset. This function is especially used in
% bagplot.m. First, the data points are standardized. If the dataset is
% too large (n > 200), then a subset is taken. On the other hand, if the
% dataset is too small (n < 10), then no bag can be constructed. The Tukey
% median (see halfmed.m) is also calculated, so that we can center the data.
% Now, the bag lies between two consecutive depth regions. Therefore, one must
% search for the right depth k. Once the k-th and the (k-1)-th depth regions
% are constructed (see isodepth), the bag is the interpolation between the
% vertices of the 2 depth regions. Finally, all the data points are caracterized
% by their position compared to the bag with an integer 0,1,2 or 3.
%
% Required input arguments:
% x : vector containing the first coordinates of all the data
% points
% y : vector containing the second coordinates of all the data
% points
%
% I/O: result = bagp(x,y);
%
% The output of bagp is a structure containing
%
% result.tukmed : Tukey median
% result.datatyp : is 2 for observations in the bag, 1 for
% observations in the fence and 0 for outliers
% result.interpol : The coordinates of the bag.
% result.depth : The halfspacedepth of each observation
%
% This function is part of the Matlab Library for Robust Analysis,
% available at:
% http://wis.kuleuven.be/stat/robust.html
%
% Last Update: 29/04/2005
n=length(x);
eps=0.0000001;
nointer=0;
ntot=n;
%
% Checking input
%
if length(x)==1
error('x is not a vector')
elseif not(length(x)==length(y))
error('The vectors x and y must have the same length.')
end
if sum(isnan(x))>=1 || sum(isnan(y))>=1
error('Missing values are not allowed')
end
if sum(x==x(1))==n
error('All data points ly on a vertical line.')
elseif sum(y==y(1))==n
error('All data points ly on a horizontal line.')
else
R=corrcoef(x,y);
if abs(R(1,2))==1
error('All data points are collineair.')
end
end
%
% Standardization of the data set
%
xmean=mean(x);
ymean=mean(y);
xdev=sqrt(var(x));
ydev=sqrt(var(y));
if xdev>eps
x=(x-xmean)/xdev;
end
if ydev>eps
y=(y-ymean)/ydev;
end
xoris=x;
yoris=y;
xori=x;
yori=y;
wx=x;
wy=y;
wx1=x;
wy1=y;
%
% If n is large, take a subset
%
nsub=sizesubset;
if n>nsub
ntot=n;
n=nsub;
a=rdraw(n,ntot);
x=wx1(a);
y=wy1(a);
wx=x;
wy=y;
end
%
% Check whether the data is in general position. If not, add a very small random
% number to each of the data points.
%
[x,y, Index, angl, ind1,ind2]=dithering(x,y);
%
numdep(1:n)=0;
kstar=0;
for i=1:n
hdep=halfspacedepth(x(i),y(i),x,y);
if hdep>kstar
kstar=hdep;
end
numdep(hdep)=numdep(hdep)+1;
end
%
% Calculation of the Tukey median
%
tukmed=halfmed(x,y,'kstar',kstar);
%
% Calculation of correct value of k
%
nc=floor(n/2);
j=kstar+1;
stay8=1;
while stay8==1
stay8=0;
j=j-1;
if numdep(kstar)<=nc
numdep(kstar)=numdep(kstar)+numdep(j-1);
stay8=1;
else
k=j+1;
numk1=numdep(kstar);
numk=numk1-numdep(k-1);
lambdanc=(nc-numk)/(numk1-numk);
end
end
%
% Calculation of the vertices of Dk
%
[kount1,Dk1,empty1]=isodepth(x,y,k);
wx1=Dk1(:,1);
if Dk1==0
wy1=0;
else
wy1=Dk1(:,2);
end
%
% Calculation of the vertices of Dk-1
%
[kount2,Dk2,empty2]=isodepth(x,y,k-1);
wx2=Dk2(:,1);
if Dk2==0
wy2=0;
else
wy2=Dk2(:,2);
end
if corrcoef(wx2,wy2)>0.98
error('Dk-1 is a line segment. Too many data points coincide.')
end
%
if n>=10
wx1=wx1-tukmed(1);
wy1=wy1-tukmed(2);
wx2=wx2-tukmed(1);
wy2=wy2-tukmed(2);
x(1:n)=x(1:n)-tukmed(1);
y(1:n)=y(1:n)-tukmed(2);
xori(1:ntot)=xori(1:ntot)-tukmed(1);
yori(1:ntot)=yori(1:ntot)-tukmed(2);
%
% Compute the angles of the data points.
%
numdattm=0;
in=1:n;
for i=1:n
if (abs(x(i))<eps)&&(abs(y(i))<eps)
numdattm=numdattm+1;
dattm(numdattm)=i;
angz(i)=1000;
else
dist=sqrt(x(i)^2+y(i)^2);
xcord=x(i)/dist;
ycord=y(i)/dist;
if abs(xcord)>abs(ycord)
if xcord>=0
angz(i)=asin(ycord);
if angz(i)<0
angz(i)=angz(i)+pi*2;
end
else
angz(i)=pi-asin(ycord);
end
else
if ycord>=0
angz(i)=acos(xcord);
else
angz(i)=pi*2-acos(xcord);
end
end
if angz(i)>=(2*pi-eps)
angz(i)=0;
end
end
end
%
% numdattm is the number of datapoints that are equal to the Tukey median
%
[angz,optel]=sort(angz);
in=in(optel);
n=n-numdattm;
%
% Compute the angles of the vertices of the depth region Dk.
%
if not(kount1==1)
ind=1:kount1;
for i=1:kount1
dist=sqrt(wx1(i)^2+wy1(i)^2);
if dist>eps
xcord=wx1(i)/dist;
ycord=wy1(i)/dist;
if abs(xcord)>abs(ycord)
if xcord>=0
angy1(i)=asin(ycord);
if angy1(i)<0
angy1(i)=angy1(i)+pi*2;
end
else
angy1(i)=pi-asin(ycord);
end
else
if ycord>=0
angy1(i)=acos(xcord);
else
angy1(i)=pi*2-acos(xcord);
end
end
if angy1(i)>=(pi*2-eps)
angy1(i)=0;
end
else
nointer=1;
end
end
[angy1,optel]=sort(angy1);
ind=ind(optel);
end
%
% Compute the angles of the vertices of the depth region Dk-1.
%
jnd=1:kount2;
for i=1:kount2
dist=sqrt(wx2(i)^2+wy2(i)^2);
if dist>eps
xcord=wx2(i)/dist;
ycord=wy2(i)/dist;
if abs(xcord)>abs(ycord)
if xcord>=0
angy2(i)=asin(ycord);
if angy2(i)<0
angy2(i)=angy2(i)+pi*2;
end
else
angy2(i)=pi-asin(ycord);
end
else
if ycord>=0
angy2(i)=acos(xcord);
else
angy2(i)=pi*2-acos(xcord);
end
end
if angy2(i)>=(pi*2-eps)
angy2(i)=0;
end
else
nointer=1;
end
end
[angy2,optel]=sort(angy2);
jnd=jnd(optel);
%
% Calculation of arrays px and py for Dk
%
if not(kount1==1)
jk=0;
wx1(kount1+1)=wx1(1);
wy1(kount1+1)=wy1(1);
angy1(kount1+1)=angy1(1);
ind(kount1+1)=ind(1);
if angz(1)<angy1(1)
j=kount1;
end
if angz(1)>=(angy1(1)-eps)
j=1;
end
for i=1:n
while (angz(i)>=(angy1(j+1)-eps))&&(jk==0)
j=j+1;
if j==kount1
jk=1;
end
if j==(kount1+1)
j=1;
end
end
if (abs(wx1(ind(j+1))-wx1(ind(j)))>eps)&&(abs(x(in(i)))>eps)
dist=y(in(i))/x(in(i))-(wy1(ind(j+1))-wy1(ind(j)))/(wx1(ind(j+1))-wx1(ind(j)));
px(i,1)=(wy1(ind(j))*wx1(ind(j+1))-wx1(ind(j))*wy1(ind(j+1)))/(wx1(ind(j+1))-wx1(ind(j)));
px(i,1)=px(i,1)/dist;
py(i,1)=px(i,1)*y(in(i))/x(in(i));
else
if abs(wx1(ind(j+1))-wx1(ind(j)))<=eps
px(i,1)=wx1(ind(j));
py(i,1)=px(i,1)*y(in(i))/x(in(i));
else
px(i,1)=0;
py(i,1)=((wy1(ind(j))-wy1(ind(j+1)))*wx1(ind(j))/(wx1(ind(j+1))-wx1(ind(j))))+wy1(ind(j));
end
end
end
end
%
% Calculation of arrays px and py for Dk-1
%
jk=0;
wx2(kount2+1)=wx2(1);
wy2(kount2+1)=wy2(1);
angy2(kount2+1)=angy2(1);
jnd(kount2+1)=jnd(1);
if angz(1)<angy2(1)
j=kount2;
end
if angz(1)>=(angy2(1)-eps)
j=1;
end
for i=1:n
while (angz(i)>=(angy2(j+1)-eps))&&(jk==0)
j=j+1;
if j==kount2
jk=1;
end
if j==(kount2+1)
j=1;
end
end
if (abs(wx2(jnd(j+1))-wx2(jnd(j)))>eps)&&(abs(x(in(i)))>eps)
dist=(y(in(i))/x(in(i)))-((wy2(jnd(j+1))-wy2(jnd(j)))/(wx2(jnd(j+1))-wx2(jnd(j))));
px(i,2)=(wy2(jnd(j))*wx2(jnd(j+1))-wx2(jnd(j))*wy2(jnd(j+1)))/(wx2(jnd(j+1))-wx2(jnd(j)));
px(i,2)=px(i,2)/dist;
py(i,2)=px(i,2)*y(in(i))/x(in(i));
else
if abs(wx2(jnd(j+1))-wx2(jnd(j)))<=eps
px(i,2)=wx2(jnd(j));
py(i,2)=px(i,2)*y(in(i))/x(in(i));
else
px(i,2)=0;
py(i,2)=((wy2(jnd(j))-wy2(jnd(j+1)))*wx2(jnd(j))/(wx2(jnd(j+1))-wx2(jnd(j))))+wy2(jnd(j));
end
end
end
if kount1==1
px(:,1)=wx1(1);
py(:,1)=wy1(1);
end
%
% Mergesort of angy1 and angy2 to obtain gamma and calculation of arrays px
% and py
%
if kount1==1
gamma=angy2;
px(:,3)=wx1;
py(:,3)=wy1;
px(:,4)=wx2;
py(:,4)=wy2;
else
ja=1;
jb=1;
i=1;
stay9=1;
while stay9==1
if angy1(ja)<=angy2(jb)
if ja<=kount1
gamma(i)=angy1(ja);
px(i,3)=wx1(ind(ja));
py(i,3)=wy1(ind(ja));
if jb==1
wxjb1=wx2(jnd(kount2));
wyjb1=wy2(jnd(kount2));
else
wxjb1=wx2(jnd(jb-1));
wyjb1=wy2(jnd(jb-1));
end
if (abs(wx2(jnd(jb))-wxjb1)>eps)&&(abs(wx1(ind(ja)))>eps)
dist=(wy1(ind(ja))/wx1(ind(ja)))-((wy2(jnd(jb))-wyjb1)/(wx2(jnd(jb))-wxjb1));
px(i,4)=(wyjb1*wx2(jnd(jb))-wxjb1*wy2(jnd(jb)))/(wx2(jnd(jb))-wxjb1);
px(i,4)=px(i,4)/dist;
py(i,4)=px(i,4)*wy1(ind(ja))/wx1(ind(ja));
else
if abs(wx2(jnd(jb))-wxjb1)<=eps
px(i,4)=wxjb1;
py(i,4)=px(i,4)*wy1(ind(ja))/wx1(ind(ja));
else
px(i,4)=0;
py(i,4)=((wyjb1-wy2(jnd(jb)))*wxjb1/(wx2(jnd(jb))-wxjb1))+wyjb1;
end
end
ja=ja+1;
i=i+1;
else
angy1(ja)=angy1(ja)+10;
end
else
if jb<=kount2
gamma(i)=angy2(jb);
px(i,4)=wx2(jnd(jb));
py(i,4)=wy2(jnd(jb));
if ja==1
wxja1=wx1(ind(kount1));
wyja1=wy1(ind(kount1));
else
wxja1=wx1(ind(ja-1));
wyja1=wy1(ind(ja-1));
end
if (abs(wx1(ind(ja))-wxja1)>eps)&&(abs(wx2(jnd(jb)))>eps)
dist=(wy2(jnd(jb))/wx2(jnd(jb)))-((wy1(ind(ja))-wyja1)/(wx1(ind(ja))-wxja1));
px(i,3)=(wyja1*wx1(ind(ja))-wxja1*wy1(ind(ja)))/(wx1(ind(ja))-wxja1);
px(i,3)=px(i,3)/dist;
py(i,3)=px(i,3)*wy2(jnd(jb))/wx2(jnd(jb));
else
if abs(wx1(ind(ja))-wxja1)<=eps
px(i,3)=wxja1;
py(i,3)=px(i,3)*wy2(jnd(jb))/wx2(jnd(jb));
else
px(i,3)=0;
py(i,3)=((wyja1-wy1(ind(ja)))*wxja1/(wx1(ind(ja))-wxja1))+wyja1;
end
end
jb=jb+1;
i=i+1;
else
angy2(jb)=angy2(jb)+10;
end
end
if i<=(kount1+kount2)
stay9=1;
else
stay9=0;
end
end
end
%
% Interpolation of two arrays px and py
%
c=3;
if nointer==1
kount1=0;
end
if nointer==0
wx2(1:(kount1+kount2))=lambdanc*px(1:(kount1+kount2),4)+(1-lambdanc)*px(1:(kount1+kount2),3);
wy2(1:(kount1+kount2))=lambdanc*py(1:(kount1+kount2),4)+(1-lambdanc)*py(1:(kount1+kount2),3);
end
if xdev>eps
xcord=(wx2(1:(kount1+kount2))+tukmed(1))*xdev+xmean;
else
xcord=(wx2(1:(kount1+kount2))+tukmed(1));
end
if ydev>eps
ycord=(wy2(1:(kount1+kount2))+tukmed(2))*ydev+ymean;
else
ycord=wy2(1:(kount1+kount2))+tukmed(2);
end
interpol(1:(kount1+kount2),1)=xcord;
interpol(1:(kount1+kount2),2)=ycord;
%
if nointer==1
if xdev>eps
xcord=(x(in(1:n))+tukmed(1))*xdev+xmean;
else
xcord=x(in(1:n))+tukmed(1);
end
if ydev>eps
ycord=(y(in(1:n))+tukmed(2))*ydev+ymean;
else
ycord=y(in(1:n))+tukmed(2);
end
datatyp(1:n,1)=xcord;
datatyp(1:n,2)=ycord;
datatyp(1:n,3)=1;
if xdev>eps
xcord=(x(dattm(1:numdattm))+tukmed(1))*xdev+xmean;
else
xcord=x(dattm(1:numdattm))+tukmed(1);
end
if ydev>eps
ycord=(y(dattm(1:numdattm))+tukmed(2))*ydev+ymean;
else
ycord=y(dattm(1:numdattm))+tukmed(2);
end
datatyp((n+1):(n+numdattm),1)=xcord;
datatyp((n+1):(n+numdattm),2)=ycord;
datatyp((n+1):(n+numdattm),3)=0;
return
end
%
% Repeat some calculations for the whole data set
%
if ntot>nsub
n=ntot;
x=xoris;
y=yoris;
%
% Angles of the data points
%
x=x-tukmed(1);
y=y-tukmed(2);
numdattm=0;
in=1:n;
for i=1:n
if (abs(x(i))<eps)&&(abs(y(i))<eps)
angz(i)=1000;
numdattm=numdattm+1;
dattm(numdattm)=i;
else
dist=sqrt(x(i)^2+y(i)^2);
xcord=x(i)/dist;
ycord=y(i)/dist;
if abs(xcord)>abs(ycord)
if xcord>=0
angz(i)=asin(ycord);
if angz(i)<0
angz(i)=angz(i)+pi*2;
end
else
angz(i)=pi-asin(ycord);
end
else
if ycord>=0
angz(i)=acos(xcord);
else
angz(i)=pi*2-acos(xcord);
end
end
if angz(i)>=(pi*2-eps)
angz(i)=0;
end
end
end
[angz,optel]=sort(angz);
in=in(optel);
n=n-numdattm;
%
% Calculation of arrays px and py for B
%
jk=0;
wx2(kount2+kount1+1)=wx2(1);
wy2(kount2+kount1+1)=wy2(1);
gamma(kount2+kount1+1)=gamma(1);
if angz(1)<gamma(1)
j=kount2+kount1;
end
if angz(1)>=(gamma(1)-eps)
j=1;
end
for i=1:n
while (angz(i)>=gamma(j+1)-eps)&&(jk==0)
j=j+1;
if j==(kount2+kount1)
jk=1;
end
if j==(kount2+kount1+1)
j=1;
end
end
if (abs(wx2(j+1)-wx2(j))>eps)&&(abs(x(in(i)))>eps)
dist=(y(in(i))/x(in(i)))-((wy2(j+1)-wy2(j))/(wx2(j+1)-wx2(j)));
px(i,1)=(wy2(j)*wx2(j+1)-wx2(j)*wy2(j+1))/(wx2(j+1)-wx2(j));
px(i,1)=px(i,1)/dist;
py(i,1)=px(i,1)*y(in(i))/x(in(i));
else
if abs(wx2(j+1)-wx2(j))<=eps
px(i,1)=wx2(j);
py(i,1)=px(i,1)*y(in(i))/x(in(i));
else
px(i,1)=0;
py(i,1)=((wy2(j)-wy2(j+1))*wx2(j)/(wx2(j+1)-wx2(j)))+wy2(j);
end
end
end
end
%
% Decide the type of each data point
%
c=3;
num=0;
num1=0;
num2=0;
num3=0;
lopen=1;
i=1;
while lopen==1
if ntot<=nsub
px(i,1)=lambdanc*px(i,2)+(1-lambdanc)*px(i,1);
py(i,1)=lambdanc*py(i,2)+(1-lambdanc)*py(i,1);
end
if (px(i,1)^2+py(i,1)^2)>eps
lambda(i)=sqrt(x(in(i))^2+y(in(i))^2)/sqrt(px(i,1)^2+py(i,1)^2);
else
if (abs(x(in(i)))<eps) && (abs(y(in(i)))<eps)
lambda(i)=0;
else
lambda(i)=c+1;
end
end
if lambda(i)<eps
num=num+1;
typ(i)=0;
elseif lambda(i)<=(1+eps)
num1=num1+1;
typ(i)=1;
elseif lambda(i)<=(c+eps)
num2=num2+1;
typ(i)=2;
elseif lambda(i)>c
num3=num3+1;
typ(i)=3;
end
i=i+1;
if i==(n+1)
lopen=0;
end
end
if xdev>eps
xcord=(x(in(1:n))+tukmed(1))*xdev+xmean;
else
xcord=x(in(1:n))+tukmed(1);
end
if ydev>eps
ycord=(y(in(1:n))+tukmed(2))*ydev+ymean;
else
ycord=y(in(1:n))+tukmed(2);
end
datatyp(1:n,1)=xcord;
datatyp(1:n,2)=ycord;
typ=typ.';
datatyp(1:n,3)=typ;
%
if not(numdattm==0)
num=num+numdattm;
typ(n+(1:numdattm))=0;
lambda(n+(1:numdattm))=0;
px(n+(1:numdattm),1)=x(dattm);
py(n+(1:numdattm),1)=y(dattm);
if xdev>eps
xcord=(x(dattm)+tukmed(1))*xdev+xmean;
else
xcord=x(dattm)+tukmed(1);
end
if ydev>eps
ycord=(y(dattm)+tukmed(2))*ydev+ymean;
else
ycord=y(dattm)+tukmed(2);
end
datatyp(n+(1:numdattm),1)=xcord;
datatyp(n+(1:numdattm),2)=ycord;
datatyp(n+(1:numdattm),3)=typ(n+(1:numdattm));
end
else
datatyp((1:n),1)=(x+tukmed(1))*xdev+xmean;
datatyp((1:n),2)=(y+tukmed(2))*ydev+ymean;
datatyp((1:n),3)=0;
interpol=0;
end
%
if xdev>eps
tukmed(1)=tukmed(1)*xdev+xmean;
end
if ydev>eps
tukmed(2)=tukmed(2)*ydev+ymean;
end
result=struct('tukmed',tukmed,'datatyp',datatyp,'interpol',interpol,'depth',hdep);
%-----------------------------
function [x,y,Index,angl,ind1,ind2]=dithering(x,y)
%DITHERING is used to check whether the data points are in general position.
% If not, then a very small random number is added to each of
% the data points. Also the angles formed by pairs of data points are
% computed here. This file is used in the functions isodepth.m, halfmed.m
% and bagp.m.
%
% Required input arguments:
% x : vector containing the first coordinates of all the data
% points
% y : vector containing the second coordinates of all the data
% points
if not(length(x)==length(y))
error('The vectors x and y must have the same length.')
end
if sum(isnan(x))>=1 || sum(isnan(y))>=1
error('Missing values are not allowed')
end
n=length(x);
i=0;
blijf3=1;
dith=1;
xorig=x;
yorig=y;
while dith==1
[xs,Index]=sort(x);
ys=y(Index);
dith=0;
while blijf3==1
blijf3=0;
i=i+1;
if (i+1)>n
blijf3=2;
else
j=i+1;
if not(xs(i)==xs(j))
blijf3=1;
else
if ys(i)==ys(j)
dith=1;
blijf3=0;
% two datapoints coincide
else
if ys(i)<ys(j)
even=Index(j);
Index(j)=Index(i);
Index(i)=even;
end
if (j+1)<=n
next=j+1;
if xs(i)==xs(next)
dith=1;
blijf3=0;
% three data points are collinear
else
blijf3=1;
end
end
end
end
end
if dith==1
fac=1000000;
ran=randn(n,2);
x=xorig+ran(:,1)/fac;
wx=x;
y=yorig+ran(:,2)/fac;
wy=y;
else
x=x;
y=y;
end
end
%
if dith==0
%
% Compute all the angles formed by pairs of data points
%
m=0;
for i=1:n
rest=(i+1):n;
spec=intersect(find(x(i)==x),rest);
hoek(spec+m-i)=pi/2;
spec=intersect(find(not(x(i)==x)),rest);
hoek(spec+m-i)=atan((y(i)-y(spec))./(x(i)-x(spec)));
p=m;
m=m+n-i;
if sum(hoek<=0)>=1
spec=find(hoek<=0);
hoek(spec)=hoek(spec)+pi;
end
ind1((p+1):m)=i;
ind2((p+1):m)=rest;
end
%
% Sort all the angles and permute ind1 and ind2 in the same way.
%
[angl,In]=sort(hoek);
%
ind1=ind1(In);
ind2=ind2(In);
%
% Test wether any three datapoints are collinear.
%
ppp=diff(angl);
k=find(ppp==0);
%
if sum(ind1(k)==ind1(k+1))>=1
% There are 3 or more datapoints collineair.
dith=1;
fac=100000000;
ran=randn(n,2);
x=xorig+ran(:,1)/fac;
wx=x;
y=yorig+ran(:,2)/fac;
wy=y;
else
x=x;
y=y;
end
end
end