# 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:

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

• Indeed, $P(X_1>X_2 , X_1>X_3,... , X_1>X_n)=P(X_1 \geq X_i\, \forall i \in \{2, ..., n\})$. This can be reduced to $P(X_1 \geq \max_{i=2}^n X_i)$ which is analysed in one of your links.
– user10525
Nov 21, 2012 at 19:49
• @Procreastinator, are you referring to the first post, suggesting I can use a joint multivariate normal distribution? Nov 21, 2012 at 20:13
• @Proc You may be overlooking that the $X_i$ can have different parameters.
– whuber
Nov 21, 2012 at 20:25
• @whuber Yes, in my second comment (deleted), you are right. I agree with your answer (+1).
– user10525
Nov 21, 2012 at 20:26
• @whuber As I discuss in a separate answer below, the different parameters for the $X_i$ are not much of a hindrance to the analysis presented in my answer to the second question listed above. Nov 21, 2012 at 21:48

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

• For future reference: this answer generalizes to the case where the $X_i$ may have entirely different distributions altogether, provided they are continuous. (Care is needed in dealing with ties for distributions that are not continuous.) The independence assumption is still essential. (Without this assumption, solutions can still be obtained with numerical integration or simulation, but both are more difficult to carry out and require the dependence to be explicitly represented.)
– whuber
Nov 21, 2012 at 20:58
• I started with a simulation, written in Python, and it didn't perform very well for $n = 12$. My preference is to use your numerical integration solution. Perhaps I can use numpy since most of my other work is in Python. Otherwise, I'll probably use Octave or R (or perhaps even C). Nov 21, 2012 at 22:01
• It is also worth noting that if all the $X_i$'s are iid continuous random variables, then the probability is $\frac{1}{n}$ without the need for numerical integration or simulation Nov 21, 2012 at 22:36
• Chris, The integration can go very quickly, because the integrand is so nicely behaved. For instance, applying Simpson's Rule in my example over the range $[-1,5]$ with intervals of $1/10$ (just $61$ evaluations, 0.006 seconds) agreed with the correct result to more than five significant figures (and using just $13$ evaluations with a spacing of $1/2$, equal to a typical SD of the $X_i$, still got the answer correct almost to four sig figs--as good as the simulation). @Dilip That's right, but even slight deviations from iid can create enormous variation among the answers when $n$ is sizable.
– whuber
Nov 21, 2012 at 22:41
• Thank you, @whuber, I appreciate your attention. I'm surprised it takes so few intervals to get such a good result. I suppose that's because the tails of the Gaussian are sufficiently small. Nov 21, 2012 at 23:50

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.

• (1) How good is the union bound? I would expect it to be poor for $n \ge 3$, because there is so much overlap among the subsets, but I haven't checked. (2) You need to adjust $n$ to $n-1$ in your first integral. (3) I think your notation gets in the way. By writing out the location-scale transformations, you actually make your answer less general than it could be: it is unnecessary to restrict the analysis to such a situation. (4) Presumably, "$Q$" should be the complementary CDF $1-\Phi$ at the end.
– whuber
Nov 21, 2012 at 22:02
• @whuber I have incorporated your comments in my revised answer. I wrote out the location-scale transformations because the OP seemed to be not able to recognize that the different means and variances are easily handled, and I agree with your comment on your own answer that the idea applies not just to normal distributions. Nov 21, 2012 at 22:32
• OK, great. Thanks for answering my question about the union bound. (+1)
– whuber
Nov 21, 2012 at 22:43