What is $P(X_1>X_2 , X_1>X_3,... , X_1>X_n)$?

Question

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:

1. The answer to this one is for only three random variables.
2. The answer to this one is for only mean 0 and variance 1.
3. 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.)

Answer 1

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.)

Answer 2

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.