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! 
 A: 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))
[1] 0.92124
   user  system elapsed 
   0.51    0.00    0.51 
> system.time(print(1 - prob(b,a), digits=5))
[1] 0.92124
   user  system elapsed 
   0.51    0.00    0.52 
> print(pnorm((mu.a - mu.b)/sqrt(2)), digits=5)
[1] 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.
