R/depth.R In WRS2: A Collection of Robust Statistical Methods

```depth<-function(U,V,m){
#
#  Compute the halfspace depth of the point (u,v) for the pairs of points
#  in the n by 2 matrix m.
#
X<-m[,1]
Y<-m[,2]
FV<-NA
NUMS<-0
NUMH<-0
SDEP<-0.0
HDEP<-0.0
N<-length(X)
P<-acos(-1)
P2<-P*2.0
EPS<-0.000001
ALPHA<-NA
NT<-0
for(i in 1:nrow(m)){
DV<-sqrt(((X[i]-U)*(X[i]-U)+(Y[i]-V)*(Y[i]-V)))
if (DV <= EPS){
NT<-NT+1
}
else{
XU<-(X[i]-U)/DV
YU<-(Y[i]-V)/DV
if (abs(XU) > abs(YU)){
if (X[i] >= U){
ALPHA[i-NT]<-asin(YU)
if(ALPHA[i-NT] < 0.0)
ALPHA[i-NT]<-P2+ALPHA[i-NT]
}
else{
ALPHA[i-NT]<-P-asin(YU)
}
}
else{
if (Y[i] >= V)
ALPHA[i-NT]<-acos(XU)
else
ALPHA[i-NT]<-P2-acos(XU)
}
if (ALPHA[i-NT] >= P2-EPS) ALPHA[i-NT]<-0.0
}
}
NN<-N-NT
if(NN<=1){
NUMS<-NUMS+depths1(NT,1)*depths1(NN,2)+depths1(NT,2)*depths1(NN,1)+
depths1(NT,3)
if(N >= 3)SDEP<-(NUMS+0.0)/(depths1(N,3)+0.0)
NUMH<-NUMH+NT
HDEP<-(NUMH+0.0)/(N+0.0)
return(HDEP)
}
ALPHA<-sort(ALPHA[1:NN])
ANGLE<-ALPHA[1]-ALPHA[NN]+P2
for(i in 2:NN){
ANGLE<-max(c(ANGLE,ALPHA[i]-ALPHA[i-1]))
}
if(ANGLE > (P+EPS)){
NUMS<-NUMS+depths1(NT,1)*depths1(NN,2)+depths1(NT,2)*depths1(NN,1)+
depths1(NT,3)
if(N >= 3)SDEP<-(NUMS+0.0)/(depths1(N,3)+0.0)
NUMH<-NUMH+NT
HDEP<-(NUMH+0.0)/(N+0.0)
return(HDEP)
}
ANGLE<-ALPHA[1]
NU<-0
for (i in 1:NN){
ALPHA[i]<-ALPHA[i]-ANGLE
if(ALPHA[i]<(P-EPS))NU<-NU+1
}
if(NU >= NN){
NUMS<-NUMS+depths1(NT,1)*depths1(NN,2)+depths1(NT,2)*depths1(NN,1)+
depths1(NT,3)
if(N >= 3)SDEP<-(NUMS+0.0)/(depths1(N,3)+0.0)
NUMH<-NUMH+NT
HDEP<-(NUMH+0.0)/(N+0.0)
return(HDEP)
}
#
#  Mergesort the alpha with their antipodal angles beta,
#  and at the same time update I, F(I), and NBAD.
#
JA<-1
JB<-1
ALPHK<-ALPHA[1]
BETAK<-ALPHA[NU+1]-P
NN2<-NN*2
I<-NU
NF<-NN
for(J in 1:NN2){
NF<-NF+1
if(JA < NN){
JA<-JA+1
ALPHK<-ALPHA[JA]
}
else
ALPHK<-P2+1.0
}
else{
I<-I+1
NN1<-NN+1
if(I==NN1){
I<-1
NF<-NF-NN
}
FV[I]<-NF
NFI<-NF-I
if(JB < NN){
JB<-JB+1
if(JB+NU <= NN)
BETAK<-ALPHA[JB+NU]-P
else
BETAK<-ALPHA[JB+NU-NN]+P
}
else
BETAK<-P2+1.0
}
}
#
#  Computation of NUMH for halfspace depth.
#
GI<-0
JA<-1
ANGLE<-ALPHA[1]
dif<-NN-FV[1]
NUMH<-min(FV[1],dif)
for(I in 2:NN){
AEPS<-ANGLE+EPS
if(ALPHA[I] <= AEPS){
JA<-JA+1
}
else{
GI<-GI+JA
JA<-1
ANGLE<-ALPHA[I]
}
KI<-FV[I]-GI
NNKI<-NN-KI
NUMH<-min(c(NUMH,min(c(KI,NNKI))))
}
NUMS<-NUMS+depths1(NT,1)*depths1(NN,2)+depths1(NT,2)*depths1(NN,1)+
depths1(NT,3)
if(N >= 3)SDEP<-(NUMS+0.0)/(depths1(N,3)+0.0)
NUMH<-NUMH+NT
HDEP<-(NUMH+0.0)/(N+0.0)
HDEP
}
```

Try the WRS2 package in your browser

Any scripts or data that you put into this service are public.

WRS2 documentation built on May 2, 2019, 4:46 p.m.