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