Rank-sum test and conditional probability

Question

Consider three samples, i.e. three different lists of numbers, A, B and C. The sample sizes of the respective samples might be different. How can I calculate the probability that a draw from sample A exceeds a draw from sample B given that a draw from sample A exceeds sample C?

I am looking for an elegant an resource saving method as proposed here. With the reply to this question, which was also posted by myself, I hoped to figure out the rest by myself, but I was wrong in that. Any helpful comments are appreciated!

Answer 1

Pr(A > B | A > C) = Pr(A > max(B,C)) / Pr(A > C). Here is Mathematica code for it. The counts are obtained using the same technique used to answer the previous question.

na = Length@a; 
(Tr@Ordering[Ordering@Join[a, Max/@Tuples@{b,c}], na] - na(na+1)/2) /
(Length@b * (Tr@Ordering[Ordering@Join[a,  c   ], na] - na(na+1)/2))

EDIT - Here's some test data:
a = {30,20,1,24,27}
b = {18,9,21,3,12,26,14,13,10,6}
c = {22,2,5,29,15,11,25,28,4,16,23,19,17,8,7}
#{a > Max[b,c]} = 468 = Tr@Ordering[Ordering@Join[a,Max/@Tuples@{b,c}],na]-na(na+1)/2
#{a > c} = 50 = Tr@Ordering[Ordering@Join[a,c],na]-na(na+1)/2

EDIT 2 - Here's a much faster algorithm:
r = Ordering@Ordering@a is a list of the ranks of a. (r is a permutation of 1,..,na.)
s = Ordering[Ordering@Join[a,b],na] is a list of the ranks of a in the combined {a,b} data.
t = Ordering[Ordering@Join[a,c],na] is a list of the ranks of a in the combined {a,c} data.
Then #{a > Max[b,c]} = (s-r).(t-r), and #{a > c} = Tr[t-r].

Edit (chameau13):

This is the corresponding R code:

    prob <- function(a,ie1,b,a1,ie2,b2,...){
    ipf <- function(a,b,...){
    m <- length(a) 
    n <- length(b)  
    if (m < n) {
    r <- rank(c(a,b), ...)[1:m] - 1:m
    } else {
    r <- rank(c(a,b), ...)[(m+1):(m+n)] - 1:n
    }
    s <- ifelse ((n+m)^2 > 2^31, sum(as.double(r)), sum(r)) / (as.double(m)*n)
    return (ifelse(m < n, s, 1-s))
    }

    expand.grid.alt <- function(seq1,seq2){
    cbind(rep.int(seq1, length(seq2)),
    c(t(matrix(rep.int(seq2, length(seq1)), nrow=length(seq2)))))}

    if(missing(a1) | missing(b2) | missing(ie2) ){
    if(ie1==">"){
    return(ipf(a,b))
    } else {
    return(ipf(b,a))  
    }  
    } else {
    if(ie1==">"){
    if(ie2==">"){  
    return(ipf(a,apply(expand.grid.alt(b,b2),1,max))/ipf(a1,b2))   
    } else {
    return(1-ipf(apply(expand.grid.alt(b,b2),1,min),a)/(1-ipf(a1,b2)))     
    }  
    } else {
    if(ie2==">"){  
    return(1-ipf(a,apply(expand.grid.alt(b,b2),1,max))/ipf(a1,b2))  
    } else {
    return(ipf(apply(expand.grid.alt(b,b2),1,min),a)/(1-ipf(a1,b2)))     
    }  
    }
    }
    }

Example:

    df  <-  
    data.frame(A=rnorm(200,1,4),B=rnorm(200,1.4,3),C=rnorm(200,0.3,5))

    #the brute force method
    df1 <- expand.grid(df$A,df$B,df$C)
    names(df1) <- c("A","B","C")

    #check if the results are correct
    all.equal(sum(df1$A>df1$B &              df1$A>df1$C)/sum(df1$A>df1$C),prob(df$A,">",df$B,df$A,">",df$C))
    all.equal(sum(df1$A<df1$B & df1$A>df1$C)/sum(df1$A>df1$C),prob(df$A,"<",df$B,df$A,">",df$C))
    all.equal(sum(df1$A>df1$B & df1$A<df1$C)/sum(df1$A<df1$C),prob(df$A,">",df$B,df$A,"<",df$C))
    all.equal(sum(df1$A<df1$B & df1$A<df1$C)/sum(df1$A<df1$C),prob(df$A,"<",df$B,df$A,"<",df$C))

    #compare execution time
    #brutforce
    ptm <- proc.time()
    df1 <- expand.grid(df$A,df$B,df$C)
        names(df1) <- c("A","B","C")
        pct <- sum(df1$A>df1$B & df1$A>df1$C)/sum(df1$A>df1$C)
    proc.time() - ptm

    user  system elapsed 
    0.930   0.214   1.145 

    #rank-sum
    system.time(prob(df$A,">",df$B,df$A,">",df$C))

    user  system elapsed 
    0.108   0.000   0.108