R/depreg.R

depreg <-
function(x,y,xout=FALSE,outfun=out,...){
#
# Compute the depth regression estimator.
# Only a single predictor is allowed in this version
#
if(is.matrix(x)){
if(ncol(x)>=2)stop("Only a single predicor is allowed")
x<-as.vector(x)
}
xy=cbind(x,y)
xy=elimna(xy)
if(xout){
flag<-outfun(xy[,1],plotit=FALSE,...)$keep
xy<-xy[flag,]
}
x=xy[,1]
y=xy[,2]
ord<-order(x)
xs<-x[ord]
ys<-y[ord]
vec1<-outer(ys,ys,"-")
vec2<-outer(xs,xs,"-")
v1<-vec1[vec2>0]
v2<-vec2[vec2>0]
slope<-v1/v2
vec3<-outer(ys,ys,"+")
vec4<-outer(xs,xs,"+")
v3<-vec3[vec2>0]
v4<-vec4[vec2>0]
deep<-NA
inter<-v3/2-slope*v4/2
temp<-matrix(c(inter,slope),ncol=2)
deep<-apply(temp,1,rdepth,x,y)
best<-max(deep)
coef<-NA
coef[2]<-mean(slope[deep==best])
coef[1]<-mean(inter[deep==best])
res<-y-coef[2]*x-coef[1]
list(coef=coef,residuals=res)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.