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