"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])
}