IsoCondDistr <- function(X=(1:length(Y)), Y,
	probs=c(0.25,0.5,0.75))
# IsoDistr(x,y,probs):
# Let (X,Y) be a pair of real random variables such that
#    F(y | x) := P(Y <= y | X = x)
# is non-increasing in x. The goal is to estimate
# F(. | x) for each entry X[i] of an input vector X,
# based on X and an additional random input vector Y
# with independent components Y[i] ~ F(. | X[i]).
# The estimator, proposed by ElBarmi and Mukerjee (2005),
# uses the Pool-adjacent-violators algorithm up to
# n = length(Y) times.
# 
# If one is only interested in isotonic regression quantiles
# for a few probability values, the program IsoRegrQuants.R
# is typically much faster and yields the same results.
# The latter fact is not obvious but can be proved rigorously;
# see Duembgen and Jongbloed (2009).
# 
# The output is a list of the following objects:
# - xs: a vector containing the different elements
#          xs[1] < xs[2] < ... < xs[mx]
#       of X.
# - ys: a vector containing the different elements
#          ys[1] < ys[2] < ... < ys[my]
#       of Y.
# - FhM: The resulting estimators Fhat(. | xs[i])
#          are discrete with support ys. Thus it suffices
#          to put out a matrix FhatM with entries
#               FhM[i,j] = Fhat(ys[j] | xs[i])
#          for 1 <= i <= mx, 1 <= j <= my.
# - QhM: For 1 <= i <= mx and 1 <= k <= length(probs),
#        QhM[i,k] is a probs[k]-quantile (type 2) of the
#        distribution Fhat(. | xs[i]).
# 
# Auxiliary program:
# - IsoMeans.R
# 
# Lutz Duembgen, October 2, 2008
# 
# Sources:
# - Hammou El Barmi and Hari Mukerjee (2005):
#   Inferences under stochastic ordering: the k-sample case.
#   Journal of the American Statistical Association 100, 252-261.
# - Lutz Duembgen and Geurt Jongbloed (2009):
#   A new link between isotonic means and quantiles
#    and some implications.
#   Preprint in preparation.
{
	n <- length(X)

	xs <- sort(X)
	xs <- xs[ xs > c(-Inf,xs[1:(n-1)]) ]
	mx <- length(xs)

	ys <- sort(Y)
	ys <- ys[ ys > c(-Inf,ys[1:(n-1)]) ]
	my <- length(ys)

	Xind <- rep(0,times=n)
	Yind <- rep(0,times=n)
	for (i in (1:n))
	{
		Xind[i] <- min((1:mx)[xs >= X[i]])
		Yind[i] <- min((1:my)[ys >= Y[i]])
	}
	# X[i] = xs[Xind[i]] and Y[i] = ys[Yind[i]]
	
	W <- rep(0,times=mx)
	FhM0 <- array(0,dim=c(mx,my))
	for (i in (1:n))
	{
		W[Xind[i]] <- W[Xind[i]] + 1
		FhM0[Xind[i],(Yind[i]:my)] <-
			FhM0[Xind[i],(Yind[i]:my)] + 1
	}

	for (i in (1:mx))
	{
		FhM0[i,] <- FhM0[i,]/W[i]
	}

	FhM <- array(0,dim=c(mx,my))
	for (j in (1:my))
	{
		FhM[,j] <- - IsoMeans(-FhM0[,j],W)
	}

	kmax <- length(probs)
	
	QhM <- array(0, dim=c(mx,kmax))
	
	for (i in (1:mx))
	{
		for (k in (1:kmax))
		{
			l1 <- min((1:my)[FhM[i,] >= probs[k]])
			l2 <- min((1:my)[FhM[i,] >  probs[k]])
			QhM[i,k] <- (ys[l1] + ys[l2])/2
		}
	}
	plot(X,Y, col="blue", xlab="x", ylab="y, Qhat(x,beta)",
		main="Isotonic regression quantiles")
	for (k in (1:kmax))
	{
		lines(xs,QhM[,k],type="l")
	}
	return(list(xs=xs,ys=ys,FhM=FhM,QhM=QhM))
}