IsoWLSRegr <- function(x=(1:length(y)), y, w=rep(1,times=length(y))) # IsoWLSRegr(x,y) # fits an isotonic function mhat to data vectors # x and y such that # sum((y - mhat(x)).^2) is minimal. # (Pool-adjacent-violators algorithm). # In case of a third input vector w with positive # entries (and the same size as x, y), IsoRegr(x,y,w) # produces an isotonic function mhat minimizing # sum(w*(y - mhat(x))^2) # # 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 # - ws: a vector with the components of w rearranged # accordingly # - mhat: a vector of the fitted values # mhat[i] = mhat(xs[i]) # # Lutz Duembgen, 30.09.2008 { ord <- order(x) xs <- x[ord] ys <- y[ord] ws <- w[ord] n <- length(ys) index <- rep(0,times=n) weight <- rep(0,times=n) # An interval of indices is represented by its left # endpoint ("index") and its "weight" (sum of w's). mhat <- rep(0,times=n) 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. weight[ci] <- sum(ws[1:j]) mhat[ci] <- sum(ys[1:j]*ws[1:j])/weight[ci] # mhat[ci] is the weighted mean of 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 } weight[ci] <- sum(ws[index[ci]:j]) mhat[ci] <- sum(ys[index[ci]:j]*ws[index[ci]:j])/weight[ci] while (ci >= 2 && mhat[ci] <= mhat[ci-1]) { # "pool adjacent violators": nw <- weight[ci-1] + weight[ci] mhat[ci-1] <- mhat[ci-1] + (weight[ci] / nw) * (mhat[ci] - mhat[ci-1]) weight[ci-1] <- nw ci <- ci-1 } } # Now define mhat for all indices: while (j >= 1) { mhat[index[ci]:j] <- mhat[ci] j <- index[ci] - 1 ci <- ci - 1 } plot(xs,ys, col="blue", xlab="x", ylab="y, mhat(x)", main="Isotonic weighted least squares") lines(xs,mhat,type="l") return(list(ord=ord,xs=xs,ys=ys,ws=ws,mhat=mhat)) }