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

Question

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.

Answer 1

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][1],sd=sds[-whichSwimmer][1])) *
    (1 - pnorm(x,mean=means[-whichSwimmer][2],sd=sds[-whichSwimmer][2])) *
    (1 - pnorm(x,mean=means[-whichSwimmer][3],sd=sds[-whichSwimmer][3])) *
    (1 - pnorm(x,mean=means[-whichSwimmer][4],sd=sds[-whichSwimmer][4]))
}

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))
[1] 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.

Answer 2

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

Answer 3

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[1], xlim[2], ...)$value
}

Answer 4

You had two questions

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.

Answer 5

I 'm interested in this, While I doubt there is an analytical solution, it's easy with Monte Carlo. So if you want to apply this to the races (human, animal, snail ...) just use Monte Carlo. But for some purposes such as if you are trying to work out the standard deviation model of your runners, you may need to do this integration many times and it may slow you down. You can however always use Simpson's method to compute the integrals.

I used Laplace distributions instead of normal and worked out the integral for two runners, in analytical form. Further than that I don't know. In (5) above what is Q ?