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