R/lsqs2.R

lsqs2 <-
function(x,y,MD=FALSE,tr=.05,plotit=TRUE){
#  cf Liu and Singh, JASA 1993, 252-260
#
if(is.list(x))x<-matl(x)
if(is.list(y))y<-matl(y)
disyx<-NA # depth of y in x
disxy<-NA # depth of x in y
if(!is.matrix(x) && !is.matrix(y)){
x<-x[!is.na(x)]
y<-y[!is.na(y)]
#
tempxx<-NA
for(i in 1:length(x)){
tempxx[i]<-sum(x[i]<=x)/length(x)
if(tempxx[i]>.5)tempxx[i]<-1-tempxx[i]
}
for(i in 1:length(x)){
temp<-sum(x[i]<=y)/length(y)
if(temp>.5)temp<-1-temp
disxy[i]<-mean(temp>tempxx)
}
tempyy<-NA
for(i in 1:length(y)){
tempyy[i]<-sum(y[i]<=y)/length(y)
if(tempyy[i]>.5)tempyy[i]<-1-tempyy[i]
}
for(i in 1:length(y)){
temp<-sum(y[i]<=x)/length(x)
if(temp>.5)temp<-1-temp # depth of y_i in x
disyx[i]<-mean(temp>tempyy)
}
qhatxy<-mean(disyx)
qhatyx<-mean(disxy)
qhat<-(qhatxy+qhatyx)/2
}
if(is.matrix(x) && is.matrix(x)){
if(!MD){
if(ncol(x)!=2 || ncol(y)!=2){
# Use approximate depth
tempyy<-fdepth(y)
temp<-fdepth(y,x)
for(i in 1:nrow(x)){
disxy[i]<-mean(temp[i]>tempyy)
}
tempxx<-NA
tempxx<-fdepth(x)
temp<-fdepth(x,pts=y)
for(i in 1:nrow(y)){
disyx[i]<-mean(temp[i]>tempxx)
}}
if(ncol(x)==2 && ncol(y)==2){
if(plotit){
plot(rbind(x,y),type="n",xlab="Var 1",ylab="VAR 2")
points(x)
points(y,pch="o")
temp<-NA
for(i in 1:nrow(x)){
temp[i]<-depth(x[i,1],x[i,2],x)
}
flag<-(temp>=median(temp))
xx<-x[flag,]
xord<-order(xx[,1])
xx<-xx[xord,]
temp<-chull(xx)
xord<-order(xx[,1])
xx<-xx[xord,]
temp<-chull(xx)
lines(xx[temp,])
lines(xx[c(temp[1],temp[length(temp)]),])
temp<-NA
for(i in 1:nrow(y)){
temp[i]<-depth(y[i,1],y[i,2],y)
}
flag<-(temp>=median(temp))
xx<-y[flag,]
xord<-order(xx[,1])
xx<-xx[xord,]
temp<-chull(xx)
flag<-(temp>=median(temp))
xord<-order(xx[,1])
xx<-xx[xord,]
temp<-chull(xx)
lines(xx[temp,],lty=2)
lines(xx[c(temp[1],temp[length(temp)]),],lty=2)
}
tempyy<-NA
for(i in 1:nrow(y))tempyy[i]<-depth(y[i,1],y[i,2],y)
for(i in 1:nrow(x)){
temp<-depth(x[i,1],x[i,2],y)
disxy[i]<-mean(temp>tempyy)
}
tempxx<-NA
for(i in 1:nrow(x))tempxx[i]<-depth(x[i,1],x[i,2],x)
for(i in 1:nrow(y)){
temp<-depth(y[i,1],y[i,2],x)
disyx[i]<-mean(temp>tempxx)
}
}}
if(MD){
mx<-apply(x,2,median)
my<-apply(y,2,median)
vx<-apply(x,2,winval,tr=tr)-apply(x,2,mean,trim=tr)+mx
vx<-var(vx)
vy<-apply(y,2,winval,tr=tr)-apply(y,2,mean,trim=tr)+my
vy<-var(vy)
tempxx<-1/(1+mahalanobis(x,mx,vx))
tempyx<-1/(1+mahalanobis(y,mx,vx))
for(i in 1:nrow(y)){
disyx[i]<-mean(tempyx[i]>tempxx)
}
tempyy<-1/(1+mahalanobis(y,my,vy))
tempxy<-1/(1+mahalanobis(x,my,vy))
for(i in 1:nrow(x)){
disxy[i]<-mean(tempxy[i]>tempyy)
}
}
qhatxy<-sum(disxy)
qhatyx<-sum(disyx)
qhat<-(qhatxy+qhatyx)/(length(disxy)+length(disyx))
}
qhatyx<-mean(disyx)
qhatxy<-mean(disxy)
list(qhatxy,qhatyx,qhat)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.