IsoMeans <- function(y, w=rep(1,times=length(y)), display=FALSE) # IsoMeans(x,y) # fits an isotonic vector mhat to a data vector # y such that # sum((y - mhat).^2) is minimal. # (Pool-adjacent-violators algorithm). # In case of a second input vector w with positive # entries (and the same size as y), IsoRegr(y,w) # produces an isotonic vector mhat minimizing # sum(w*(y - mhat)^2) # # Lutz Duembgen, October 3, 2008 { n <- length(y) 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 j <- 1 # ci is the number of the interval considered currently, # where we start with {1}. The index j is the right endpoint # of the current interval. weight[ci] <- w[1] mhat[ci] <- y[1] # mhat[ci] is the weighted mean of y-values within # the current interval. while (j < n) { j <- j + 1 # a new index intervall {j} is created: ci <- ci + 1 index[ci] <- j weight[ci] <- w[j] mhat[ci] <- y[j] 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 } if (display) { plot((1:n),y, col="blue", xlab="i", ylab="y_i, mhat_i", main="Isotonic weighted least squares") lines((1:n),mhat,type="l") } return(mhat) }