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 > \max\{B,C,D\}$$
 A: 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) &= P(X_1 \ge X_2, X_1 \ge X_3, \ldots, X_1 \ge X_N)\\
&= \int_{-\infty}^\infty P(X_1 \ge X_2, X_1 \ge X_3, \ldots, X_1 \ge X_N, X_1 = x) dx\\
&= \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\\
&= \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\\
&= \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\\
&= \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 >= 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.
A: If you want to do it with numerical integration, that's straightforward enough. The joint pdf is 
$$f(x_1,\dots,x_N)=\prod_{i=1}^N \frac{1}{\sqrt{2 \pi \sigma_i}} e^{-\frac{(x_i-\mu_i)^2}{2\sigma_i^2}}.$$
Now supposing you want $X_1$ to be the biggest, then you just integrate over the region $x_1>x_2,x_1>x_3,\dots,x_1>x_N$. This is 
$$\int_{-\infty}^\infty \int_{-\infty}^\infty \dots \int_{\max \{ x_2,\dots,x_N \}}^\infty f(x_1,\dots,x_N) dx_1 dx_2 \dots dx_N.$$
However, I will warn that if $N$ is even slightly large (say, 10), this integral will probably be faster to compute by Monte Carlo integration than by direct integration. This is a manifestation of the "curse of dimensionality". If you really want to use direct integration, one approach would be Gauss-Hermite quadrature: given the $m$ nodes $x_j$ and $m$ weights $w_j$ for one dimensional Gauss-Hermite quadrature (meaning that $\sum_{j=1}^m w_j f(y_j) \approx \int_{-\infty}^\infty \frac{1}{\sqrt{2 \pi}} f(x) e^{-x^2/2} dx$), your integral can be approximated as
$$\sum_{i_1=1}^m w_{i_1} \sum_{i_2=1}^m w_{i_2} \dots \sum_{i_N=1}^m w_{i_N} f(\mu_1+\sigma_1 x_{i_1},\mu_2+\sigma_2 x_{i_2},\dots,\mu_N+\sigma_N x_{i_N})$$
where $f(x_1,\dots,x_N)=1$ if $x_1=\max x_i$ and $0$ otherwise. Since this function is discontinuous, the convergence rate will probably be just decent, not great. Moreover, as you can see you will need $m^N$ evaluations to use this rule at all, which is large when $N$ is large even if $m$ is just $2$.
A: To ask for $p_{A} = P(A > max(B,~C,~D))$ is equivalent to ask for $P(A > B~ \land ~A > C~ \land ~A > D)$.

Edit: I think my initial answer was wrong, so I reworked on it, without being able to show a final solution yet.
$P(A > B~ \land ~A > C~ \land ~A > D) = P(A > B) \cdot P( A > C~ \land ~A > D|A > B) = P(A > B) \cdot P(A > C|A > B) \cdot P(A > D|A > B, A>C)$. 
A: 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.
