Imagine we have $k$ Gaussian RVs $$ X_i \sim N(\mu_i, \sigma_i^2) \text{ for } i=1, \ldots, k $$ and we sample from each of them independently to produce a vector, $\vec{x} = (x_1, \ldots, x_k)$.

For one of the Gaussian RVs, say $X_j$, I am interested in computing the probability that it exceeds all others, i.e. $$ \Pr\left\{ \cap_{i\not= j} \, X_j > X_i \right\}. $$

I know I can use Monte Carlo sampling to estimate this probability. But are there any closed-form analytical methods or approximations?

  • $\begingroup$ You can condition on the value of $X_j$ and then use the total probability formula. $\endgroup$ – Ian Jan 24 at 4:20
  • $\begingroup$ @Math1000 What is written there is both deleted and incorrect. $\endgroup$ – Ian Jan 24 at 4:21
  • $\begingroup$ @Ian - I see, but then I need to integrate over all possible conditional values for $X_j$, which doesn't seem nice. Is there an easy way to approximate this integral? $\endgroup$ – ted Jan 24 at 4:30
  • $\begingroup$ A similar question has been asked before. It's not an exact duplicate -- in that case all the $X_i (i \neq j)$ are equi-distributed, which makes it potentially easier, but the question asked for the rank of $X_j$, which makes it potentially harder. Regardless, despite two years and a bounty, nobody was able to come up with anything beyond the obvious. :( In limit cases (e.g. $\sigma_j$ very large) some qualitative stmts can be made, but not in general... $\endgroup$ – antkam Jan 24 at 7:24

The distribution of the maximum of several normal random variables is approximately known

One method you could use to approximate this would be to use some results from extreme value theory. As is mentioned in this answer to Is there an analytical expression for the distribution of the max of a normal k sample? the distribution of the maximum of several normal random variables is approximated by a Gumbel distribution. Similar answers can also be found for Asymptotic distribution of maximum order statistic of IID random normals. While I think these may be tangible for uncorrelated random normals, you get several awkward issues when you introduce correlation which I am not sure can be trivially overcome (Distribution of the maximum of two correlated normal variables).

This reduces the problem down to a single integral

Given you have an expression for the distribution of the maximum, it is then relatively easy to write the probability as a two dimensional integral with the joint distribution (using appropriate integration limits corresponding to $X_j > \max(\{X_i\}_{i\neq j})$). Given the normal should be independent to the Gumbel the joint distribution should be separable and this should collapse nicely to just a single integral. This integral is likely quite ugly, but given the functions are all fairly nice it might be feasible to find where the integrand is maximal and just approximate this region to approximate the probability. While this doesn't solve the problem, it should hopefully reduce it to a simpler problem.

Bounding the error is your biggest challenge if you wish to compare it to Monte Carlo

The biggest issue will be bounding the error from using the Gumbel distribution, and then on top of that having to likely approximate the resultant integral.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.