wilcox_R: What the package does (short line)

mbmreg <-
function(x,y,iter=20,bend=2*sqrt(ncol(x)+1)/nrow(x),xout=FALSE,outfun=outpro,...){
#
# Compute a bounded M regression estimator using
# Huber Psi and Schweppe weights with
# regression outliers getting a weight of zero.
#
# This is the modified M-regression estimator in Chapter 8
#
# The predictors are assumed to be stored in the n by p matrix x.
#
x<-as.matrix(x)
p1<-ncol(x)+1
p<-ncol(x)
xy<-cbind(x,y)
xy<-elimna(xy)
x<-xy[,1:p]
y<-xy[,p1]
if(xout){
m<-cbind(x,y)
flag<-outfun(x,plotit=FALSE,...)$keep
m<-m[flag,]
x<-m[,1:p]
y<-m[,p1]
}
x<-as.matrix(x)
if(is.matrix(y)){
if(ncol(y)==1)y=as.vector(y)
}
x1<-cbind(1,x)
library(MASS)
reslms<-lmsreg(x,y)$resid
sighat<-sqrt(median(reslms^2))
sighat<-1.4826*(1+(5/(length(y)-ncol(x)-1)))*sighat
if(sighat==0)warning("The estimated measure of scale, based on the residuals using lms regression, is zero")
temp<-ifelse(sighat*reslms>0,abs(reslms)/sighat,0*reslms)
wt<-ifelse(temp<=2.5,1,0)
init<-lsfit(x,y,wt)
resid<-init$residuals
nu<-sqrt(1-hat(x1))
low<-ncol(x)+1
for(it in 1:iter){
ev<-sort(abs(resid))
scale<-median(ev[c(low:length(y))])/qnorm(.75)
rov<-(resid/scale)/nu
psi<-ifelse(abs(rov)<=bend,rov,bend*sign(rov))  # Huber Psi
wt<-nu*psi/(resid/scale)
wt<-ifelse(temp<=2.5,wt,0)
new<-lsfit(x,y,wt)
if(abs(max(new$coef-init$coef)<.0001))break
init$coef<-new$coef
resid<-new$residuals
}
resid<-y-x1%*%new$coef
if(abs(max(new$coef-init$coef)>=.0001))
paste("failed to converge in",iter,"steps")
list(coef=new$coef,residuals=resid,w=wt)
}