# Probability that randomly chosen value from one distribution is greater than randomly chosen value from another distribution

Say I have $n$ values sampled from two distributions, $A$ and $B$ . That is, I have a sample $A_1, A_2, \dots, A_n$ and a sample $B_1, B_2, \dots, B_n$. How would I go about finding $P\left(A_i>B_j\right)$ for any given value $i\in(1,\dots,n)$ and $j\in(i,\dots,n)$?

I know I could get a bootstrapped solution fairly easily using the following code:

bootstrapProcedure <- function(A, B, sample.size = 100) {
# Calculates the fraction of times a sample of size (sample.size) from A is
# greater than a sample of the same size from B (both drawn with replacement).
#
# Args:
#   A: vector of values for sample 1
#   B: vector of values for sample 2
#   sample.size: integer of the size of the bootstrapped sample to draw
#
# Returns:
#   The fraction of times the sample from A is greater than the sample from B
mean(sample(A, sample.size, replace = T) > sample(B, sample.size, replace = T))
}

# Draw 2 populations
A <- rnorm(1000, mean = 1, sd = 2)
B <- rnorm(1000, mean = 2, sd = 4)

# Get the bootstrapped probability 1,000 times
replicate(1000, bootstrapProcedure(A, B))


but it seems like there should be a simple, analytical solution to this. Any ideas how I should go about finding it?

• Is $A_i$ independent of $B_j$? – TrynnaDoStat Mar 4 '15 at 21:07
• The way you phrase this question suggests your data are not independent, for otherwise $P(A_i\gt B_j)$ should not depend on $i$ or $j$ at all. How do they depend on each other, then? It's also unclear whether you are asking how to estimate such probabilities from data or derive them from assumptions about $A$ and $B$. Could you please edit the question to address those issues? – whuber Mar 4 '15 at 21:07
• It's not clear to me what relation the 'nonparametric' tag has to the question. I suggest you replace it with 'independence' (if, as it seems, independence actually applies, which you should state explicitly in your question). – Glen_b Mar 5 '15 at 6:09

Your bootstrap simulation suggests, the A's and B's are independent; I will assume so.

Note that $P(A>B)=P(A-B>0)$.

In the case where A and B are also normal, $D=A-B \sim N(\mu_A-\mu_B,\sigma^2_A+\sigma^2_B)$

In that case

\begin{eqnarray*} P(D>0) &=& P(\frac{D-\mu_D}{\sigma_D}>\frac{0-\mu_D}{\sigma_D}) \\ &=& P(Z>-\frac{\mu_A+\mu_B}{\sqrt{\sigma^2_A+\sigma^2_B}})=P(Z<\frac{\mu_A+\mu_B}{\sqrt{\sigma^2_A+\sigma^2_B}})\\ &=&\Phi\left(\frac{\mu_A+\mu_B}{\sqrt{\sigma^2_A+\sigma^2_B}}\right) \end{eqnarray*}

In the case of other distributions there may be no simple "closed" form (you'll get one for a few distributions but you can't expect in in general). For specific instances the value of the probability can be calculated via numerical convolution.

(If you have bivariate normality but not independence, you can do a similar calculation.)

• +1. Others reading this thread may also be interested in reading about the Mann-Whitney U test / Wilcoxon rank-sum test, which could be considered basically the hypothesis-test counterpart to the estimation problem that the OP asked about. – Jake Westfall Oct 6 '17 at 21:04
• @Jake that's a great point. Indeed, the statistic may also be turned into an estimate of the probability (divide the Mann-Whitney statistic -- the number of $i,j$ pairs where $A_i$ exceeded $B_j$ - by the total number of such pairs of observations, to get the corresponding "sample proportion"). – Glen_b Oct 6 '17 at 21:21

To use standard packages, organize your data with "labels" as 0 for sample a and 1s for sample b then use the measurements as the "score". AUC is then an estimate for $P\left(A_i<B_j\right)$