huber2 <- function (y, k = 1.345, initmu = median(y), 
					s = mad(y), iters = FALSE, tol = 1e-08) 
{
require(robustbase)
y <- y[!is.na(y)]
n <- length(y)
mu <- as.numeric(initmu)
s <- as.numeric(s)
Niter <- 300
Converged <- FALSE
if (s == 0) 
stop("cannot estimate scale: s is zero for this sample")
for (i in 0:Niter) {
if (iters) s <- mad(y,mu)  #s ignored, MAD of current residuals
yy <- Mpsi((y - mu)/s, cc = k, psi="huber")*s
mu1 <- sum(yy)/n
mu <- mu + mu1
if (abs(mu1) < tol * s ) {Converged <- TRUE; Niter <- i; break}
}
list(mu = mu, s = s, Niter = Niter, Converged = Converged)
}


cauchy <- function (y, initmu = median(y), 
					s = mad(y), iters = FALSE, tol = 1e-08) 
{
y <- y[!is.na(y)]
n <- length(y)
mu <- as.numeric(initmu)
s <- as.numeric(s)
Niter <- 300
Converged <- FALSE
if (s == 0) 
stop("cannot estimate scale: s is zero for this sample")
for (i in 0:Niter) {
if (iters) s <- mad(y,mu)  #s ignored, MAD of current residuals
x <- (y - mu)/s
yy <- 2 * x/(1 + x^2) * s
mu1 <- sum(yy)/n
mu <- mu + mu1
if (abs(mu1) < tol * s ) {Converged <- TRUE; Niter <- i; break}
}
list(mu = mu, s = s, Niter = Niter, Converged = Converged)
}


biweight <- function (y, k = 4.685061, initmu = median(y), 
					  s = mad(y), iters = FALSE, tol = 1e-08) 
{
require(robustbase)
y <- y[!is.na(y)]
n <- length(y)
mu <- as.numeric(initmu)
s <- as.numeric(s)
Niter <- 300
Converged <- FALSE
if (s == 0) 
stop("cannot estimate scale: s is zero for this sample")
for (i in 0:Niter) {
if (iters) s <- mad(y,mu)  #s ignored, MAD of current residuals
yy <- Mpsi((y - mu)/s, cc = k, psi="bisquare") * s
mu1 <- sum(yy)/n
mu <- mu + mu1
if (abs(mu1) < tol * s ) {Converged <- TRUE; Niter <- i; break}
}
list(mu = mu, s = s, Niter = Niter, Converged = Converged)
}


hampel <- function (y, a = 2, b = 4, c = 8, initmu = median(y),
s = mad(y), iters = FALSE, tol = 1e-08) 
{
require(robustbase)
y  <- y[!is.na(y)]
n  <- length(y)
mu <- as.numeric(initmu)
s  <- as.numeric(s)
Niter <- 300
Converged <- FALSE
if (s == 0) 
stop("cannot estimate scale: s is zero for this sample")
for (i in 0:Niter) {
if (iters) s <- mad(y,mu) #s ignored, MAD of current residuals
yy  <- Mpsi((y - mu)/s, cc = c(a, b, c), psi="hampel") * s
mu1 <- sum(yy)/n
mu  <- mu + mu1
if (abs(mu1) < tol * s ) {Converged <- TRUE; Niter <- i; break}
}
list(mu = mu, s = s, Niter = Niter, Converged = Converged)
}


skipped.mean <- function (y, k = 2, initmu = median(y), 
						  s = mad(y), iters = FALSE, tol = 1e-08) 
{
y <- y[!is.na(y)]
n <- length(y)
mu <- as.numeric(initmu)
s <- as.numeric(s)
Niter <- 300
Converged <- FALSE
if (s == 0) 
stop("cannot estimate scale: s is zero for this sample")
for (i in 0:Niter) {
if (iters) s <- mad(y,mu)  #s ignored, MAD of current residuals
x <- (y - mu)/s
yy <- ifelse(abs(x) > k, 0, x) * s
mu1 <- sum(yy)/n
mu <- mu + mu1
if (abs(mu1) < tol * s ) {Converged <- TRUE; Niter <- i; break}
}
list(mu = mu, s = s, Niter = Niter, Converged = Converged)
}


winsorized.mean<-function(x, trim=0)
{
x <-sort(x[!is.na(x)])
n<-length(x)
if ((trim < 0) | (trim>0.5) ) 
stop("cannot estimate: alpha<0 or alpha>0.5")
ntr<-trunc(trim*n)
(ntr*x[ntr+1]+sum(x[(ntr+1):(n-ntr)])+ntr*x[n-ntr])/n
}


hodges.lehmann <- function(x)
{
x <- x[!is.na(x)]
diffs <- outer(x, x, "+")
diffs <- diffs[!upper.tri(diffs)] 
median(diffs)/2
}


tls.KB<-function(formula, data,  alpha=0.05)
{
require(quantreg)
resid1 <-residuals(rq(formula,data,tau=alpha))
resid2 <-residuals(rq(formula,data,tau=1-alpha))
c1 <- c(resid1 >= 0)
c2 <- c(resid2 <= 0)
coefficients(lm(formula,data[c(c1 & c2),]))
} 


rfit0 <- function(formula, data=NULL, tau=0.5) 
{
require(Rfit)
n <- dim(model.frame(formula,data))[1]
my.phi0 <- function(u, param = tau)
{
s1 = param[1]
c1 = ceiling(s1*n)
ifelse(u < c1/(n+1), s1 - 1, ifelse(u > c1/(n+1), s1,
	   c1-1-s1*(n-1)))
}
my.Dphi0 <- function(u, param = tau)
{
s1 = param[1]
ifelse(u < s1, 0, ifelse(u > s1, 0, 1))
}
myscores0 <- new("scores", phi = my.phi0, Dphi = my.Dphi0,
				 param = tau)
rfit(formula=formula, data=data,scores=myscores0)
}


disp0 <- function (beta, x, y, tau=0.5) 
{
require(Rfit)
my.phi0 <- function(u, param)
{
s1 = param[1]
c1 = ceiling(s1*n)
ifelse(u < c1/(n+1), s1 - 1, ifelse(u > c1/(n+1), s1, 
c1-1-s1*(n-1)))
}
x <- as.matrix(x)
e <- y - x %*% beta
r <- rank(e, ties.method = "first")/(length(e) + 1)
my.phi0(r,param=tau) %*% e
}

"LADLS"<-function(ydata, xdata, resultls, resultl1, int=T, kern=F)
{
 n <- length(ydata)
 xdata <- as.matrix(xdata)
 p <-length(xdata[1,])
 if(!int)
   p<-p-1
 if(missing(resultl1))
    resultl1 <- rq(ydata~xdata,tau=0.5,ci=F)
 if(missing(resultls))
    resultls <- lm(ydata~xdata)
 if (kern) sn<-s.kernel(xdata,ydata,sqrt(3),int=int)
 else sn<-s.histog(xdata,ydata,int=int)
 absresl1 <- abs(resultl1[[2]])
 e10 <- 1/(sn * (n - p -1)) * sum(absresl1)
 sigma02 <- (1/((sn^2) * (n - p - 1))) * sum(resultl1[[2]] * resultl1[[2]])
 if(int)
   xdata<-cbind(rep(1,n),xdata)
 if(e10 >= 0.5) {
   delta0 <- 1
   tn <- resultl1[[1]]
 }
 if(e10 < 0.5 && e10 >= 2 * sigma02) {
   delta0 <- 0
   tn <- resultls[[1]]
 }
 if(e10 < 0.5 && e10 < 2 * sigma02) {
    delta0 <- (4 * sigma02 - 2 * e10)/(4 * sigma02 - 4 * e10 + 1)
    tn0 <- delta0 * resultl1[[1]] + (1 - delta0) * resultls[[1]]
    tn<-nlminb(start=tn0, objective=my2.fcn, matx=xdata,
                      vecy=ydata,delta=delta0)$parameters
 }
 tn <- as.vector(tn)
 if(int)
    xlabel <- c("Intercept", paste("X", 1:p, sep= ""))
 else xlabel <- paste("X", 1:(p+1), sep = "")
 names(tn)<-xlabel
 return(coef=tn)
}


"my.fcn" <- function(t,matx,vecy,delta,gama,sn,k)
         { (1-delta)*2/gama* sum(((((vecy-matx%*%t)/sn)^2)/2)[c(compare(rep(k,length(vecy)),
           (vecy-matx%*%t)/sn))>=0]) +
         (1-delta)*2/gama* sum((k*abs((vecy-matx%*%t)/sn)^(-k*k/2) )[c(compare(rep(k,length(vecy)),
           (vecy-matx%*%t)/sn))<0]) +
         delta*sum(abs(vecy-matx%*%t)/sn)
         }


"my2.fcn" <- function(t,matx,vecy,delta)
          { (1-delta)* sum(((vecy-matx%*%t)^2)) + delta*sum(abs(vecy-matx%*%t))
          }

 "s.histog"<-function(xdata,ydata,int=T)
{
 nun <- (((1/8)^(2/5)) * (9/2)^(1/5) * (length(ydata))^(-1/5))
 cxdata <- scale(xdata, center = T, scale = F)
 quantplus <- rq(ydata~cxdata,  tau=(0.5 + nun),ci=F)$coef[1]
 quantminus <- rq(ydata~cxdata, tau=(0.5 - nun),ci=F)$coef[1]
 sn <- max((quantplus-quantminus)/(2 * nun),0.0001)
 return(sn)
}


"s.kernel"<-function(xdata,ydata,b,int=T)
{
 nun <- (((1/8)^(2/5)) * (9/2)^(1/5) * (length(ydata))^(-1/5))
 sol<-rq(ydata~xdata,ci=F)$sol[c(1,3),]
 sol<-sol[,abs((0.5-sol[1,])/nun)<b]
 if (sol[1,1]>0){
              xl<-0.5-b*nun
              xl<-c(xl,rq(ydata~xdata,tau=xl,ci=F)$coef[1])
              sol<-cbind(xl,sol)
 }
 if (sol[1,length(sol[1,])]<1){
              xu<-0.5+b*nun
              xu<-c(xu,rq(ydata~xdata,int=int,tau=xu,ci=F)$coef[1])
              sol<-cbind(sol,xu)
 }
 nsol<-length(sol[1,])
 kernel<--1.5/b^3*(0.5*sol[1,]/nun-sol[1,]^2/(2*nun))
 snker<-0
 for(i in 1:(nsol-1))
  {snker<-snker + sol[2,i]*(kernel[i+1]-kernel[i])}
 snker<-max(snker/nun^2,0.0001)
 return(snker)
}


rank.test.meas<-
function(x, y, score = "wilcoxon")
{
n<-length(x[,1]) 
xp<-apply(x,2,mean)
x2<-x-matrix(rep(xp,n),nrow=n,byrow=T)
Q=t(x2)%*%x2/n
por<-rank(y)
if(score == "wilcoxon") {
skor<-por/(n+1)-0.5
A2<-12
}
 else if(score == "normal") {
  skor<-qnorm(por/(n+1))
  A2<-1
}
else if(score == "sign") {
  skor<-0.5*sign(por/(n+1)-0.5)
  A2<-1
}
sn<-n^(-0.5)*apply(x2*skor,2,sum)
t(sn)%*%solve(Q)%*%sn*A2
}