N different normal distributions - probability A is largest Following on from this question, I now have a harder version I need to solve. If we have N numbers selected from N normal distributions, each with its own mean and variance, what is the probability that a particular number is the largest. So for example we choose $A,B,C$ and $D$ from distributions with means $M_a,M_b,M_c,$ and $M_d$ and variances $V_a,V_b,V_c$ and $V_d$, then what is the probability that:

$$A &gt; \max\{B,C,D\}$$
Thank you Ian for pointing out my mistake (independence issues). I've got ANOTHER answer and some code for you. I tested it out properly with R via Monte Carlo simulation. I'll provide the R code and the Mathematica code. I used Mathematica to numerically evaluate the integral I derive below.

So let's say we have $N$ distributions, and each $X_i$ is distributed according to the $i^{th}$ distribution. Let's assume without loss of generality that we are interested in $P(X_1 \ge \max_i X_i)$. I'm going to be a little loose with the notation, but I think it'll be clear.

$$\begin{align}
P(X_1 \ge \max_i X_i) &amp;= P(X_1 \ge X_2, X_1 \ge X_3, \ldots, X_1 \ge X_N)\\
&amp;= \int_{-\infty}^\infty P(X_1 \ge X_2, X_1 \ge X_3, \ldots, X_1 \ge X_N, X_1 = x) dx\\
&amp;= \int_{-\infty}^\infty P(X_1 \ge X_2, X_1 \ge X_3, \ldots, X_1 \ge X_N| X_1 = x) f_{X_1}(x) dx\\
&amp;= \int_{-\infty}^\infty P(X_1 \ge X_2, X_1 \ge X_3, \ldots, X_1 \ge X_N| X_1 = x) f_{X_1}(x) dx\\
&amp;= \int_{-\infty}^\infty P(X_1 \ge X_2|X_1)P(X_1 \ge X_3|X_1)\ldots P(X_1 \ge X_N| X_1 = x) f_{X_1}(x) dx\\
&amp;= \int_{-\infty}^\infty f_{X_1}(x) \prod_{i=2}^N F_{X_i}(x) dx\\
\end{align}$$
where $f$ denotes the PDF and $F$ is the CDF. At this point, it's reasonable to evaluate the integral numerically.

Here is some Mathematica code to show an example of how to evaluate the above expression. Say


$X_1$ ~ N(1,2)
$X_2$ ~ N(0,1)
$X_3$ ~ N(2,1)
$X_4$ ~ N(1.5,3)


Then to get $P(X_1 \ge \max_i X_i)$ we just run the following in Mathematica

NIntegrate[
 CDF[NormalDistribution[0, 1], x]*CDF[NormalDistribution[2, 1], x]*
  CDF[NormalDistribution[1.5, 3], x]*
  PDF[NormalDistribution[1, 2], x], {x, -Infinity, Infinity}]


We can get answer of 0.216147.

Now let's check our answer with R. Here's code for running a Monte Carlo simulation. Note that it took me around 2 minutes to run.

x1 = rnorm(10^8, 1, 2)
x2 = rnorm(10^8, 0, 1)
x3 = rnorm(10^8, 2, 1)
x4 = rnorm(10^8, 1.5, 3)
indicators = ifelse(x1 &gt;= pmax(x2, x3, x4), 1, 0)
mean(indicators)


I just ran it and I got 0.2161365, which is pretty close to what we expect. I tried different numbers of variables and different distributions, and it always worked. Hopefully I understood your problem correctly and this helps you.
Let us say the $A$ we are looking for is $A_k$ and we want to see if it is largest in the set $\{A_1,\cdots,A_n\}$. My initial thought is to fix $A_k$ and then calculate the probability that none of the other $A_i$ is pairwise larger and then use complementary probability:

$$P\left(\max_i A_i\leq A_k| A_k\right) =1-\prod_{\forall i \neq k} P\left(A_i\leq A_k\right)$$

Then what would be needed is to sum or integrate over all $A_k$:

$$P\left(\max_i A_i \leq A_k\right) = \sum_{\forall A_k} P\left(\max_i A_i \leq A_k|A_k\right)\cdot P(A_k)$$

If you have a table for normal distribution function, this becomes easy, just the product of a series of table lookups ( maybe you want interpolation in between table entries ). You will need to calculate which entry of the table to check for each $A_i$ but it will be a a really nice function of its mean and variance. This will be nice to express with vectors too, in case you want to try it in a computational language like R or Matlab or Python or Octave.