# Probability that a draw from one sample is greater than a draw from a another sample

I have a sample of b for id=1 and I want to compare this sample with the sample of id={2,3} of a i.e. I want to compute the probability that a draw from sample b (id=1) is greater than a draw from sample a (id={2,3}).

    library(data.table)
df <- data.table(id=c(rep(1,500),rep(2,500),rep(3,500)),a=rnorm(1500),b=rnorm(1500,1))
ot <- df[id!=1,a]
ow <- df[id==1,b]
sa <- df[id!=1,sample(a,500)]
#all permutations of eig and sam
fall <- expand.grid(ow,sa)
#all permutations of eig and other
all  <- expand.grid(ow,ot)

#1. best possible empirical measure
sum(all[,1]>all[,2])/nrow(all)
#2.
sum(fall[,1]>fall[,2])/nrow(fall)
#3.
sum(ow>sa)/length(sa)


Nb. 1 would obviously be the best solution as in my real data set the distributional properties of a and b are unknown. Nb. 1 is however not feasible because I would have to repeat this procedure for every ID. My data set is huge.

The solution I used until now was Nb. 3: I a draw a random sample from a (id={2,3}) in the size of the sample for b (id=1). But this is of course imprecise. Solution Nb. 2 would be a compromise but is also not feasible computation-wise.

I would like to find a data table (this is a faster implementation of data frames in the form of the R package "data.table") solution because those are normally much faster.

I didn't post this question on Stackoverflow, because I also hoped that one of you could maybe hint me to a mathematical trick to overcome this problem, but of course also coding advises are welcome. Thanks in advance!

• Can you explain what you mean when you say you always get a bias when comparing samples from group 1 with samples from groups 3 and 10 combined? Oct 1, 2013 at 9:39
• @M.Berk I completely rewrote the question. I hope it is clearer now. Oct 1, 2013 at 12:24
• Can you write in plain math instead of R code? Oct 1, 2013 at 12:30
• This isn't even working R code! -- What is data.table? There are no variables a or b to make sense of the second through fourth lines, either.
– whuber
Oct 1, 2013 at 12:59
• @Penguin_Knight Thanks, I know how to search for packages (at rseek.org and other places). The point is that people should not have to search the Web in order to be able to read questions here! If they have to, the question needs clarification.
– whuber
Oct 1, 2013 at 13:26

There is a simple efficient solution. It uses the ideas common to all rank-sum tests, such as the Wilcoxon tests. This answer derives the solution and provides an R implementation.

The code in the question simulates data that have a vanishingly small chance of exhibiting any ties at all between a sample of group $a$ and a sample of group $b$, so let's assume there exist no ties.

Let there be $m$ elements in $a$ (and they can have ties among each other) and $n$ elements in $b$ (which also may have ties). Let $A$ be the random variable modeled by drawing one element randomly and uniformly from group $A$ and similarly let $B$ be the random variable for one draw from group $B$. The desired value (as I interpret the question) is the chance that $A$ exceeds $B$.

Notice that the test of whether a realization of $A$ exceeds one of $B$ is a simple comparison. Thus, the problem is unchanged if we replace all elements of $a$ and $b$ by their indexes when the two sets are sorted in increasing order. These indexes are their ranks, provided that ties are resolved in some arbitrary manner (that is, do not average the ranks of any groups of ties).

For example, let $a = (0,2)$ and $b = (1,1,3)$. Sorting the combined two (multi)sets gives the sequence $(0,1,1,2,3)$. The indexes of the values coming from $a$ are $1$ and $4$ while the indexes of the values from $b$ are $2, 3,$ and $5$.

Compute the chance $\Pr(A\gt B)$ by summing over the possible values of $A$, each of which has the probability $1/m$. Let the ranks of these values be $r_1 \lt r_2 \lt \cdots \lt r_m$. Suppose $r_i$ is picked. Then the chance, conditional on this selection, that $B$ has a smaller value equals the number of smaller values in $b$ divided by $n$. The number of smaller values altogether in both $a$ and $b$ is, by definition, $r_i-1$, but we know exactly $i-1$ of them (namely, $r_1, r_2, \ldots, r_{i-1}$) are in $a$. Thus

$$\Pr(A \gt B | A = r_i) = \frac{1}{n}\left(r_i-1 - (i-1)\right) = \frac{1}{n}\left(r_i-i\right)$$

entailing

$$\Pr(A\gt B) = \sum_{i=1}^m \Pr(A\gt B | A=r_i)\Pr(A=r_i) = \frac{1}{mn}\sum_i \left(r_i-i\right).$$

In the example, $\Pr(A\gt B) = \frac{1}{2\times 3}((1-1) + (4-2)) = \frac{2}{6}$ and (reversing the roles of $a$ and $b$ as a quick check) $\Pr(B\gt A) = \frac{1}{3\times 2}((2-1) + (3-2) + (5-3)) = \frac{4}{6} = 1 - \frac{2}{6}$ as one would expect.

This calculation (when implemented as a general-purpose algorithm) requires sorting all $m+n$ values to find their ranks and then summing either $m$ or $n$ values (for efficiency, one would pick whichever is smaller). Therefore the computational burden is $O((m+n)\log(m+n)),$ and can be reduced to $O(\min(m\log(n), n\log(m)))$ when the larger of $a$ and $b$ is already sorted. That's pretty efficient.

When there are ties between elements of $a$ and $b$, the idea to reduce the question to rank sums and then compute a sum over conditional probabilities still works, but the calculations of the conditional probabilities get more complicated.

### R Code

R will have trouble with calculations that overflow its integer data type. The following solution handles that possibility.

prob <- function(a, b, ...) {
# Returns chance that a random sample of a will exceed one of b
# (assuming no ties between elements of a and b)
# Optional args are passed to rank to control handling of NAs and
# how to resolve any ties.
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))
}


To emulate the data in the question, let's simulate sets of normally distributed values. When the elements of $a$ come from a normal distribution with mean $\mu$ and standard deviation $\sigma$ and those of $b$ come from a normal distribution with mean $\nu$ and SD $\sigma$, we may analytically compute that $\Pr(A\gt B)$ (prior to simulating the elements) equals $\Phi(\frac{\mu-\nu}{\sigma\sqrt{2}})$ where $\Phi$ is the cumulative standard normal distribution function. This enables us to test prob, as in the following:

set.seed(17)
m <- 10^6; n <- 10^4
mu.a <- 0; mu.b <- -2
a <- rnorm(m, mu.a)
b <- rnorm(n, mu.b)
system.time(print(prob(a,b), digits=5))
system.time(print(1 - prob(b,a), digits=5))
print(pnorm((mu.a - mu.b)/sqrt(2)), digits=5)


The output is

> system.time(print(prob(a,b), digits=5))
 0.92124
user  system elapsed
0.51    0.00    0.51
> system.time(print(1 - prob(b,a), digits=5))
 0.92124
user  system elapsed
0.51    0.00    0.52
> print(pnorm((mu.a - mu.b)/sqrt(2)), digits=5)
 0.92135


It shows that the computation time does not depend on the order in which a and b are provided to prob. The computation time of $1/2$ second is reasonably quick (for over one million numbers). The close agreement of $0.92124$ and $0.92135$ is evidence in favor of the correctness of this solution.

This solution can easily be iterated over groups using the usual R idioms for looping.

• awesome! how were you able to write such a reponse so fast after my last edit!? thank you so much! Oct 1, 2013 at 18:03
• I would have one more question: can you also extend this method to conditional probabilities like P(A>B|A>C)? Oct 3, 2013 at 17:40
• I think so, although obviously it would be more complicated.
– whuber
Oct 3, 2013 at 17:49
• i don't want to bother you further with my questions but could you maybe hint me to literature which could help me solving this question myself? Oct 3, 2013 at 18:07
• The literature is vast, because underlying these calculations is the entire field of combinatorics. You might consider asking for recommendations suitable for your interests and mathematical background. The Mathematics site could be a good place for such advice.
– whuber
Oct 3, 2013 at 18:25