What is $P(X_1>X_2 , X_1>X_3,... , X_1>X_n)$? All $X$ are mutually independent and from normal distributions, each with its own mean and variance. If it's easier, $P(X_1 \geq X_i \forall i \in \{1, ..., n\})$ is fine although I suspect it's the same. If it matters, $n$ is between 5 and 20.
I found three similar questions:


*

*The answer to this one is for only three random variables.

*The answer to this one is for only mean 0 and variance 1.

*I'm unsure if this one applies. If it does, I don't know how to apply it. Its top answer is for three random variables.


(This is not homework.)
 A: My answer to the second question you list 
has a simple form of the more general result given by @whuber, but is readily adapted
to the general case.  Instead of
$$P(X_1 > \max X_i  \mid X_1 = \alpha) = \prod_{i=2}^n P\{X_i < \alpha \mid X_1 = \alpha\} 
= \left[\Phi(\alpha)\right]^{n-1}$$
which applies when the $X_i$ are independent $N(0,1)$ random variables, we have
$$P(X_1 > \max X_i  \mid X_1 = \alpha) = \prod_{i=2}^n P\{X_i < \alpha \mid X_1 = \alpha\} 
= \prod_{i=2}^n \Phi\left(\frac{\alpha-\mu_i}{\sigma_i}\right)$$
since the $X_i$ are independent $N(\mu_i, \sigma_i^2)$ random variables, and instead of 
$$P(X_1 > \max X_i)  
= \int_{-\infty}^{\infty}\left[\Phi(\alpha)\right]^{n-1}
\phi(\alpha-\mu)\,\mathrm d\alpha$$
we have
$$P(X_1 > \max X_i)  
= \int_{-\infty}^{\infty}\prod_{i=2}^n \Phi\left(\frac{\alpha-\mu_i}{\sigma_i}\right)
\frac{1}{\sigma}\phi\left(\frac{\alpha-\mu_1}{\sigma_1}\right)\,\mathrm d\alpha$$
where $\Phi(\cdot)$ and $\phi(\cdot)$ are the cumulative distribution function
and probability density function of the standard normal random variable.
This is just whuber's answer expressed in different notation.  
The complementary probability 
$P(X_1 < \max X_i) = P\{(X_1 < X_2) \cup \cdots \cup (X_1 < X_n)$ can also
be bounded above by the union bound discussed in my answer to the other question.
We have that
$$\begin{align*}
P(X_1 < \max X_i) &=  P\{(X_1 < X_2) \cup \cdots \cup (X_1 < X_n)\\
&\leq \sum_{i=2}^n P(X_1 < X_i)\\
&= \sum_{i=2}^n Q\left(\frac{\mu_1 - \mu_i}{\sqrt{\sigma_1^2 + \sigma_i^2}}\right)
\end{align*}$$
since $X_i-X_1 \sim N(\mu_i-\mu_1,\sigma_i^2+\sigma_1^2)$.  Note that
$Q(x) = 1-\Phi(x)$ is the complementary  standard normal
distribution function. The union bound is very tight when
 $\mu_1 \gg \max \mu_i$ and the variances are roughly
comparable even for large $n$, but for small $n$, the bound
can exceed $1$ and thus be useless.
A: For $n \gt 2$ this needs numeric integration, as indicated in several of the links.
To be explicit, let $\phi_i$ be the PDF of $X_i$ and $\Phi_i$ be its CDF.  Conditional on $X_1 = t$, the chance that $X_1 \gt X_i$ for the remaining $i$ is the product of the individual chances (by independence):
$$\Pr(t \ge X_i, i=2,\ldots,n) = \Phi_2(t)\Phi_3(t)\cdots\Phi_n(t).$$
Integrating over all values of $t$, using the distribution function $\phi_1(t)dt$ for $X_1$, gives the answer
$$= \int_{-\infty}^{\infty} \phi_1(t) \Phi_2(t)\cdots\Phi_n(t)dt.$$
For $n=20$, the integral takes 5 seconds with Mathematica, given vectors $\mu$ and $\sigma$ of the means and SDs of the variables:
\[CapitalPhi] = MapThread[CDF[NormalDistribution[#1, #2]] &, {\[Mu], \[Sigma]}];
\[Phi] = PDF[NormalDistribution[First[\[Mu]], First[\[Sigma]]]];
f[t] := \[Phi][t] Product[i[t], {i, Rest[\[CapitalPhi]]}]
NIntegrate[f[t], {t, -Infinity, Infinity}]

The value can be checked (or even estimated) with a simulation. In the same five seconds it takes to do the integral, Mathematica can do over 2.5 million iterations and summarize their results:
m = 2500000;
x = MapThread[RandomReal[NormalDistribution[#1, #2], m] &, {\[Mu],\[Sigma]}]\[Transpose];
{1, 1./m} # & /@ SortBy[Tally[Flatten[Ordering[#, -1] & /@ x]], First[#] &]

For instance, we can generate some variable specifications at random:
{\[Mu], \[Sigma]} = RandomReal[{0, 1}, {2, n}];

In one case the integral evaluated to $0.152078$; a simulation returned 

{{1, 0.152387}, ... }

indicating that the first variable was greatest $0.152387$ of the time, closely agreeing with the integral.  (With this many iterations we expect agreement to within a few digits in the fourth decimal place.)
