R/lts1reg.R

In musto101/wilcox_R: What the package does (short line)

lts1reg <-
function(x,y,tr=.2,h=NA){
#
# Compute the least trimmed squares regression estimator.
# Only a single predictor is allowed in this version
#
if(is.na(h))h<-length(x)-floor(tr * length(x))
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]
val<-NA
inter<-v3/2-slope*v4/2
for(i in 1:length(slope)){
#risk<-(y[vec2>0]-slope[i]*x[vec2>0]-inter[i])^2
risk<-(y-slope[i]*x-inter[i])^2
risk<-sort(risk)
val[i]<-sum(risk[1:h])
}
best<-min(val)
coef<-NA
coef[2]<-mean(slope[val==best])
coef[1]<-mean(inter[val==best])
res<-y-coef[2]*x-coef[1]
list(coef=coef,residuals=res)
}