IsoQuantRegr <- function(x=(1:length(y)), y, probs=c(0.25,0.5,0.75)) # IsoQuantRegr(x,y,probs) # fits for each entry probs[k] (a number strictly between 0 and 0) # an isotonic function mhat(.,k) to data vectors # x and y such that # sum(rho(mhat(x) - y, probs[k])) is minimal, # where # rho(t,p) := (1 - 2p)t + |t| . # (Pool-adjacent-violators algorithm with suitable sample # quantiles). # # The output is a list of the following objects: # - ord: a permutation of the vector (1:length(x)) # such that x[ord[i]] is non-decreasing in i. # - xs: a vector of the order statistics of x # - ys: a vector with the components of y rearranged # accordingly # - mhat: a matrix containing the fitted values # mhat[i,k] = mhat(xs[i],k) # # Lutz Duembgen, 01.10.2008 { ord <- order(x) xs <- x[ord] ys <- y[ord] n <- length(ys) kmax <- length(probs) mhat <- array(0, dim=c(n,kmax)) for (k in (1:kmax)) { zs <- ys # zs contains the entries of ys, partially ordered. index <- rep(0,times=n) # An interval of indices is represented by its left # endpoint ("index"). ci <- 1 index[ci] <- 1 # ci is the number of the interval considered currently. j <- 1 while (j < n && xs[j+1] == xs[j]) { j <- j + 1 } # j is the right endpoint of the current index interval. zs[1:j] <- sort(ys[1:j]) mhat[ci,k] <- LocalQuantile(zs[1:j],probs[k]) # mhat[ci,k] is a sample-probs[k]-quantile of the y-values # within the current interval. while (j < n) { j <- j + 1 # a new index intervall with left endpoint j is created: ci <- ci + 1 index[ci] <- j while (j < n && xs[j+1] == xs[j]) { j <- j + 1 } zs[index[ci]:j] <- sort(ys[index[ci]:j]) mhat[ci,k] <- LocalQuantile(zs[index[ci]:j], probs[k]) while (ci >= 2 && mhat[ci,k] <= mhat[ci-1,k]) { # "pool adjacent violators": zs[index[ci-1]:j] <- LocalMerge( zs[index[ci-1]:(index[ci]-1)], zs[index[ci]:j]) ci <- ci-1 mhat[ci,k] <- LocalQuantile(zs[index[ci]:j], probs[k]) } } # Now define mhat for all indices: while (j >= 1) { mhat[index[ci]:j,k] <- mhat[ci,k] j <- index[ci] - 1 ci <- ci - 1 } } plot(xs,ys, col="blue", xlab="x", ylab="y, mhat(x)", main="Isotonic regression quantiles") for (k in (1:kmax)) { lines(xs,mhat[,k],type="l") } return(list(ord=ord,xs=xs,ys=ys,mhat=mhat)) } LocalMerge <- function(x,y) # computes a vector z containing the components of x and y # in non-decreasing order. The vectors x and y itself are # supposed to have non-decreasing components. { m <- length(x) n <- length(y) z <- rep(0,times=(m+n)) i <- 1 j <- 1 k <- 1 while (i <= m && j <= n) { if (x[i] <= y[j]) { z[k] <- x[i] i <- i+1 } else { z[k] <- y[j] j <- j+1 } k <- k+1 } if (i <= m) { z[k:(m+n)] <- x[i:m] } if (j <= n) { z[k:(m+n)] <- y[j:n] } return(z) } LocalQuantile <- function(yord, beta) { n <- length(yord) k <- ceiling(beta*n) ell <- floor(beta*n + 1) return((yord[k] + yord[ell])/2) }