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