# Rank-sum test and conditional probability

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!

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

• Thank you very much for your reply! Would you please be so nice to explain your Mathematica code a bit. When I just apply the R-code from the last example along Pr(A > B | A > C) = Pr(A > max(B,C)) / Pr(A > C) it yields the wrong result. What does for example Max/@Tuples@{b,c} do? Thank you! Oct 5, 2013 at 17:45
• ok, I got it: Max/@Tuples@{b,c} corresponds to apply(expand.grid(x,y),1,max); I will write the code in R and add it to your response. Oct 5, 2013 at 19:58