"rqx"<-
function(x, y, tau = 0.5, int = T, max.it = 50)
{
# Barrodale and Roberts quantile regression algorithm - lite
#			author:  R. Koenker, 
#                       date:    June, 1998
# The algorithm of Barrodale and Roberts (1974) for solving discrete l1 regression 
# problems has served as the standard simplex algorithm for quantile regression
# since the mid-70's.  For large problems recent developments of interior point
# methods have proven to be superior, but for small to moderate problems the
# approach of Barrodale and Roberts is still very competitive.  Unfortunately,
# the fortran version of original algorithm is somewhat obscure.  In this pure
# R version I have tried to reduce the algorithm to its essentials -- at each
# iteration we find the Cauchy direction i.e. the direction of steepest descent
# and then we take this direction until it no longer reduces the value of the
# of the objective function.  This "Cauchy Step" involves solving a one dimensional
# "through the origin" rq problem.  This can be formulated as a weighted quantile
# problem that is solved by the function wquantile().  A more efficient version
# of this algorithm should eventually be coded in fortran with an O(n) algorithm
# for solving these subproblems using the approach of Floyd and Rivest (1975), if
# possible.  Meanwhile this is simply a heuristic device.
#
	if(int) x <- cbind(1, x)
	p <- ncol(x)
	n <- nrow(x)	#Phase I -- find a random (!) initial basis
	h <- sample(1:n, size = p)
	it <- 0
	repeat {
		it <- it + 1
		Xhinv <- solve(x[h,  ])
		bh <- Xhinv %*% y[h]
		rh <- y - x %*% bh	
	#find direction of steepest descent along one of the edges
		g <-  - t(Xhinv) %*% t(x[ - h,  ]) %*% c(tau - (rh[ - h] < 0))
		g <- c(g + (1 - tau),  - g + tau)
		ming <- min(g)
		if(ming >= 0 || it > max.it)
			break
		h.out <- seq(along = g)[g == ming]
		sigma <- ifelse(h.out <= p, 1, -1)
		if(sigma < 0)
			h.out <- h.out - p
		d <- sigma * Xhinv[, h.out]	
	#find step length by one-dimensional minimization
		xh <- x %*% d
		step <- wquantile(xh, rh, tau)
		h.in <- step$k
		h <- c(h[ - h.out], h.in)
	}
	if(it > max.it)
		warning("non-optimal solution: max.it exceeded")
	return(bh)
}
"rho.tau"<- function(r, tau = 0.5) tau * pmax(r, 0) + (1 - tau) * pmax( - r, 0)
"wquantile"<- function(x, y, t = 0.5) {
#weighted quantile by brute force (full sorting) for scalar quantile regression
# through the origin model:  Q_{y_i} (t)  = x_i b
#
	ord <- order(y/x)
	b <- (y/x)[ord]
	wabs <- abs(x[ord])
	k <- sum(cumsum(wabs) < ((t - 0.5) * sum(x) + 0.5 * sum(wabs)))
	list(b = b[k + 1], k = ord[k + 1])
}