# $P(X_1 < \min(X_i,\ldots, X_n))$ across different normal random variables

I have a set of mutually independent normal distributions $X_1$ to $X_5$ (with means and standard deviations) which represent finishing times for swimmers over a certain distance. The actual data is as follows:

$$X_1(60, 3.0)$$ $$X_2(61, 1.0)$$ $$X_3(58, 2.3)$$ $$X_4(63, 2.4)$$ $$X_5(61, 1.7)$$ So swimmer 1 ($X_1$) has a mean finishing time of 60 seconds with a standard deviation of 3.0 seconds.

Question 1: What is the probability of an event where $X_i$ finishes first. e.g.

$$P(X_1 \lt X_i, i=2,\ldots,n)$$

Question 2: If I calculate this for all swimmers, can I simply order the results to determine the most probable finishing order?

This is not homework.

Based on the answers to this Cross Validated question, I have tried to solve this problem based on the first answer. i.e.

$$\Pr(X_1 \le X_i, i=2,\ldots,n) = \int_{-\infty}^{\infty} \phi_1(t) [1 - \Phi_2(t)]\cdots[1 - \Phi_n(t)]dt$$

Where $\phi_i$ is the PDF of $X_i$ and $\Phi_i$ is its CDF.

Based on this formula, the results I obtained were:

$$\Pr(X_1 \le X_i, i=2\ldots5) = 0.259653$$ $$\Pr(X_2 \le X_i, i=1,3\ldots5) = 0.214375$$ $$\Pr(X_3 \le X_i, i=1\ldots2, 4\ldots5) = 0.611999$$ $$\Pr(X_4 \le X_i, i=1\ldots3, 5) = 0.0263479$$ $$\Pr(X_5 \le X_i, i=1\ldots4) = 0.0697597$$ However, the probabilities add to 1.182135 when they should add to 1.0. I’m not sure if the formula is incorrect or my implementation of the integral (I used Excel and the trapezoidal method).

I also attempted to solve the problem using Dillip’s method (from the above mentioned question) as follows:

\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*}

However, the probability results were much greater than 1 in most cases so abandoned this approach. By the way, what exactly does $\max X_i$ mean?

Any assistance in calculating the probability would be appreciated.

• $\max$ is short for 'maximum'; it's the converse of $\min$, and gives the largest, rather than the smallest value. Oct 30 '13 at 12:39

I always find it best in these situations to run a Monte Carlo simulation to check (roughly) what the correct answer should be. Here is some R code for doing that:

doRace <- function()
{
times <- rnorm(mean=c(60,61,58,63,61),sd=c(3,1,2.3,2.4,1.7),n=5)
winner <- which.min(times)
winner
}

winners <- replicate(n=10000,expr=doRace())
table(winners) / length(winners)


Which gives the following output for me (of course you will get slightly different answers depending on the state of your random number generator):

winners
1      2      3      4      5
0.2573 0.0317 0.6108 0.0282 0.0720


This indicates that the issue is with swimmer 2, as these results otherwise agree well with yours. I suspect you just have an incorrect cell reference somewhere. Note that a reasonable resolution to the problem is to use the Monte Carlo simulation not just as a verification method but as the final implementation for calculating the probabilities. After all, numerical integration is itself an approximate and computationally expensive procedure.

In order to be absolutely sure, we can use the integrate() function in R. First define the integral:

integral<- function(x,whichSwimmer)
{
means <- c(60,61,58,63,61)
sds <- c(3,1,2.3,2.4,1.7)

dnorm(x,mean=means[whichSwimmer],sd=sds[whichSwimmer]) *
(1 - pnorm(x,mean=means[-whichSwimmer],sd=sds[-whichSwimmer])) *
(1 - pnorm(x,mean=means[-whichSwimmer],sd=sds[-whichSwimmer])) *
(1 - pnorm(x,mean=means[-whichSwimmer],sd=sds[-whichSwimmer])) *
(1 - pnorm(x,mean=means[-whichSwimmer],sd=sds[-whichSwimmer]))
}


Then we can calculate the probability for each swimmer in turn:

>integrate(integral,whichSwimmer=1,lower=0,upper=100)
0.2596532 with absolute error < 2.5e-05

>integrate(integral,whichSwimmer=2,lower=0,upper=100)
0.03223977 with absolute error < 6.4e-06

>integrate(integral,whichSwimmer=3,lower=0,upper=100)
0.6119995 with absolute error < 1.5e-06

>integrate(integral,whichSwimmer=4,lower=0,upper=100)
0.02634785 with absolute error < 1.4e-06

>integrate(integral,whichSwimmer=5,lower=0,upper=100)
0.06975967 with absolute error < 8.1e-05


Which gives very good agreement with the Monte Carlo simulation.

Note that although you can technically give lower and upper bounds of negative/positive infinity to integrate() I found that this caused the procedure to break down, giving clearly incorrect results.

EDIT: I've just noticed you had a second question regarding the most likely ordering of swimmers. Again, we can easily check whether the intuition about just ordering the probability of each swimmer winning is correct by running a Monte Carlo simulation. We just need to adapt the sampling function to return the order of the swimmers rather than only the winner:

doRace <- function()
{
times <- rnorm(mean=c(60,61,58,63,61),sd=c(3,1,2.3,2.4,1.7),n=5)
finishOrder <- order(times)

paste(finishOrder,collapse="")
}

finishOrders <- replicate(n=1e6,expr=doRace())
which.max(table(finishOrders) / length(finishOrders))


I get the ouput:

31254
50


In other words, the most likely order is $3,1,2,5,4$ which is not the same as ordering the swimmers by their probability of winning!

For me, this is another reason to prefer the Monte Carlo approach as the final implementation as you can easily answer this and other questions - e.g. what is each swimmer's probability of finishing second or, given that swimmer $1$ finishes first, what is the most likely ordering of the remaining swimmers?

EDIT 2: To be able to answer these other questions, we need to change the sampling function again, this time to return the complete order in which the swimmers finish:

doRace <- function()
{
times <- rnorm(mean=c(60,61,58,63,61),sd=c(3,1,2.3,2.4,1.7),n=5)
finishOrder <- order(times)

finishOrder
}

finishOrders <- replicate(n=1e6,expr=doRace())


finishOrders is a matrix where each column corresponds to a single simulated race, the first row gives the winner of each race, the second row the second placed swimmer of each race and so on. So, to get the probability that each swimmer finishes second we do:

> table(finishOrders[2,]) / ncol(finishOrders)

1        2        3        4        5
0.271749 0.198205 0.235460 0.075165 0.219421


To find the most likely order given that swimmer $1$ wins the race is a little more fiddly. First, extract all races where the first row is equal to $1$:

finishOrdersWhen1WinsRace <- finishOrders[,finishOrders[1,]==1]


Then turn the finish order of each race from a vector of numbers into a character string so we can use the table function to find the most frequent one:

> which.max(table(apply(finishOrdersWhen1WinsRace,2,paste,collapse="")))
13254
8


In other words, given that swimmer $1$ wins the race, the most likely order of the remaining swimmers is $3,2,5,4$, which occurs:

> max(table(apply(finishOrdersWhen1WinsRace,2,paste,collapse="")) / ncol(finishOrdersWhen1WinsRace))
 0.2341227


$23.4\%$ of the time.

I'm not sure whether a Monte Carlo approach is the only way to answer these questions but it seems likely that even if you can obtain closed-form expressions for the probabilities you'll need to rely on numerical integration like you did to find the winning probabilities.

• Thanks very much for your comprehensive answer. I hadn't even thought of doing a monte carlo analysis. You were right about having a problem in my Excel SS. I fixed it and now get the correct results. Oct 31 '13 at 9:19
• You mentioned the calculation of other stats such as each swimmer's probability of finishing second or, given that swimmer 1 finishes first, what is the most likely ordering of the remaining swimmers? Are you able to shed some light on how these might be calculated - either R code and/or a description of the process would be great. I assume these kinds of stats can only be calculated using a monte carlo approach as opposed to implementing a single formula? Oct 31 '13 at 9:33
• Jason, I've added the R code for answering those two additional questions. Once you've changed the sampling function to return the complete order of the swimmers you should be able to answer pretty much any question like this although manipulation of the returned matrix can get fiddly. The other thing to consider is whether you want to return the finish times themselves in order to answer questions such as about the distribution of the winning time or the largest gap between times or so on. Nov 1 '13 at 9:20
• Thanks again for your feedback. I've been able to download R and have run your code successfully. Nov 4 '13 at 22:36

The answers above don't provide a closed-form solution; nor do the ones to the related question here. I will try to give an analytical answer. To fix notation, let:

$X_i$ be the independent random variables, with $X_i \sim N(\mu_i,\sigma_i^2)$

$Y_{i-1}:= X_1 - X_i, i =2,\ldots, n$.

Finally, let $e\in R^{n-1}$ be such that $(e)_i=1$. Then

$P(X_1 \le X_i, i=2,\ldots,n) = P(Y_i \ge 0, i=2,\ldots,n)$

It is $Y_i \sim N(\mu_1-\mu_i, \sigma_1^2+\sigma_i^2)$ and moreover (straightforward calculation) $\text{cov} (Y_i, Y_j)=\sigma_1^2$.

so the covariance matrix $\Sigma$ of $Y$ is given by

$\Sigma = \text{diag}(\sigma_2^2,\ldots, \sigma_n^2) + \sigma_1^2 ee'$

Define $\nu_{i-1} := \mu_i-\mu_1, i =2,\ldots, n$. Now we use the standard affine transformation of a multivariate standard normal to obtain an arbitrarily distributed multivariate normal. Let $\xi \sim N(0, I)$ be a $(n-1)$-dimensional standard normal. We have $Y =_\text{dist} \nu + Q^{1/2} \xi$.

The square root $Q^{1/2}$ is uniquely identified.

The probability we want to estimate becomes

$P(Y \ge 0) = P(\xi \ge -Q^{-1/2} \nu) = \prod_{i=1}^{n-1} \bar \Phi( (-Q^{-1/2} \nu)_i)$

where the last equality follows from the independence of the $\xi_i$ an d $\bar \Phi$ is the complement of the cumulative distribution of the standard normal.

Last observation: $Q^{1/2}$ is trivial in the case where $\kappa:=\sigma_2=\ldots=\sigma_n$. I am omitting the derivation, but point out that there is one eigenvector $e$ with eigenvalue $\kappa^2+\sigma_1^2$. The other eigenvector lie in $[e]^\perp$ and have all eigenvalues equal to $\kappa^2$.

Question 1: What is the probability of an event where $X_i$ finishes first.

Your proposed solution look correct but, as you say, you clearly have an error in implelentation as the probablities do not add up to $1$. M. Berk has shown that it is likely to be an issue with Swimmer 2.

Question 2: If I calculate this for all swimmers, can I simply order the results to determine the most probable finishing order?

Not quite - if you have two swimmers with the same mean time in the middle of the group then the chance of one beating the other is $\frac12$, but the one with the higher standard deviation is more likely to be the overall winner: in your particular example, the chance of $X_2$ beating $X_5$ is $0.5$ but $X_5$ is more likely than $X_2$ to be the overall winner, so $X_5$ is less likely than $X_2$ to be third because of its higher standard deviation.

Using simulation, the most likely finishing order is $X_3,X_1,X_2,X_5,X_4$ (with a probability about 9.0%) above $X_3,X_1,X_5,X_2,X_4$ (about 8.0%), $X_1,X_3,X_2,X_5,X_4$ (about 6.1%), $X_3,X_2,X_5,X_1,X_4$ (about 5.8%), $X_1,X_3,X_5,X_2,X_4$ (also about 5.8%), $X_3,X_5,X_1,X_2,X_4$ (about 4.5%) and other less likely outcomes. Your idea would have predicted 31524, the second most likely outcome.

• Hi Henry, thanks for your response. The finishing order analysis is interesting. Did you do a monte carlo analysis such as that mentioned by M.Berk above to determine the the probabilities? Is this what is meant by simulation? How did you calculate the probability of a particular finishing order eg. 9.0% for order: 31254. Did you use R or some other app? Oct 31 '13 at 9:25
• I used R in for the simulation in a similar way to M Berk. The previous paragraphs were logical argument Oct 31 '13 at 15:47

The formula at https://stats.stackexchange.com/a/44142/919 is correct.

Below is R code to implement it (as a function p.max). The code to input the data of the question and apply p.max is short:

X <- data.frame(Mean = c(60, 61, 58, 63, 61),
SD = c(3, 1, 2.3, 2.4, 1.7))
X$p <- sapply(seq_len(nrow(X)), p.max, mu=-X$Mean, sigma=X$SD)  (Note how supplying the negatives of the times enables p.max to compute the chances of being the smallest.) This plot illustrates the five densities, colored according to the winning chances: The sum of the computed chances is within $$10^{-10}$$ of equaling $$1.$$ # # Given X[i] ~ Normal(mu[i], sigma[i]), # compute the chance that the jth component is largest. # ... is passed to integrate. # p.max <- function(j, mu, sigma, Z=7, ...) { if (is.infinite(mu[j])) { if (mu[j] < 0) { if (all(is.infinite(mu) & mu < 0)) return(1 / length(mu)) return(0) } return(1 / sum(is.infinite(mu) & mu > 0)) } # # Choose units that make X[j] standard Normal. # mu <- (mu[-j] - mu[j]) / sigma[j] sigma <- sigma[-j] / sigma[j] # # Set integration limits. # xlim <- range(c(-Z, Z, mu-Z*sigma, mu+Z*sigma)) # # Integrate. # h <- Vectorize(function(x) exp(dnorm(x, log=TRUE) + sum(pnorm(x, mu, sigma, log.p=TRUE)))) integrate(h, xlim, xlim, ...)$value
}

• Any chances of getting this in python? Adding suitable alternatives reference, albertotb.com/… Jan 15 at 22:10