Let $X = \left( {{X_1},...,{X_n}} \right) \sim \mathcal{N}\left( {{\mathbf{\mu }},{\mathbf{\Sigma }}} \right)$ be a Gaussian random vector and $I = \mathop {\arg \max }\limits_{i = 1,n} {X_i}$.

$I$ has probability mass function

$\mathbb{P}\left( {I = i} \right) = \mathbb{P}\left( {{X_i} = \mathop {\max {X_j}}\limits_{j = 1,n} } \right) = \mathbb{P}\left( {{X_i} - \mathop {\max {X_j}}\limits_{j \ne i} > 0} \right)$

and mathematical expectation

$\mathbb{E}I = \sum\limits_{i = 1}^n {i\mathbb{P}\left( {I = i} \right)} $

Generally speaking, for large $n$ and arbitrary covariance matrix ${\mathbf{\Sigma }}$ , computing $\mathbb{E}I$ is very difficult because it requires the numerical evaluation of high-dimensional normal orthant integrals. So, apart from the IID and INID cases with a diagonal covariance matrix ${\mathbf{\Sigma }}$, banded covariance matrices and degenerate cases such as ${\mu _j} \gg {\mu _{i \ne j}}$ , under which conditions on ${\mathbf{\Sigma }}$ (e.g. correlation decay) can we get simple, easy-to-evaluate numerical approximations to $\mathbb{E}I$ (and $\mathbb{V}I$ as well)?

The covariance matrices ${\mathbf{\Sigma }}$ I'm interested in look like this:

Covariance matrix

Until now, I’ve not been able to find anything about this problem.

Related questions:

https://mathoverflow.net/questions/153039/maximal-component-of-a-multivariate-gaussian-distribution Expectation of the softmax transform for Gaussian multivariate variables

  • 1
    $\begingroup$ By "argmax" are you referring to what we call the "maximum order statistic" if we view the random vector as a sample, or to something else? $\endgroup$ Jul 20, 2018 at 14:28
  • $\begingroup$ @AlecosPapadopoulos Yes the argmax is the maximum order statistic also noted ${X_{(n)}}$ in extreme value theory $\endgroup$
    – user193726
    Jul 20, 2018 at 14:31
  • 1
    $\begingroup$ It is not at all certain that we know of conditions that will allow for "easy numerical approximations" of the $n$-dimensional normal integral with correlation. $\endgroup$ Jul 20, 2018 at 14:46
  • $\begingroup$ @AlecosPapadopoulos ...Hence my question. For instance, I know from some papers in extreme value theory that if the correlations/covariances decay sufficiently fast, we can get simplified formulae for the maximum ${X_{(n)}}$. Therefore, we might get simplified formulae for the argmax as well, but I've not been able to find a single paper about this problem. Besides, sorry, ${X_{(n)}}$ is the maximum, not the argmax which is the index of the maximum. $\endgroup$
    – user193726
    Jul 20, 2018 at 14:56
  • 1
    $\begingroup$ @UthsavChitra Unfortunately, not. $\endgroup$
    – user193726
    Jul 1, 2019 at 9:06

1 Answer 1


You can use the law of large numbers to approximate your expectation pretty easily.

Edit: Analyticaly you can multiply a bunch of normal cdf evaluations together. For $i > 0$

\begin{align*} \mathbb{P}\left( {I = i} \right) &= \mathbb{P}\left( {{X_i} = \mathop {\max {X_j}}\limits_{j = 1,n} } \right) \\ &= \prod_{j \ne i} \mathbb{P}\left( {{X_i} > X_j } \right) \\ &= \prod_{j \ne i} \mathbb{P}\left( Y_j \le 0 \right) \end{align*}

where $Y_j \sim \text{Normal}(\mu_j - \mu_i, \Sigma_j + \Sigma_i)$.

If we code $i=0$ as the event that no data point is clearly a maximum, then $\mathbb{P}\left( I = 0 \right) = 1 - \sum_{k=1}^n \mathbb{P}\left( I = k \right)$

  • $\begingroup$ Yes, sampling from $X$ is the easy way. But what if $n=81$ (as above), $n=1000$ or $n=100000$ (my own $n$'s are virtually unbounded)? Can we efficiently sample (e.g. in polynomial time) an arbitrary large multivariate normal distribution (I should know, but I don't) and how many samples do we need in order to properly estimate $\mathbb{E}I$? If there is an efficient Monte Carlo method for arbitrary $n$ and ${\mathbf{\Sigma }}$ then the problem is solved. Otherwise, we'd like to get efficient, analytic approximations depending on the particular structure of ${\mathbf{\Sigma }}$ $\endgroup$
    – user193726
    Jul 20, 2018 at 19:09
  • 1
    $\begingroup$ Both of those sound like separate questions. If you want to simulate quickly for higher dimensional data, you might have to write more bespoke code and introduce parallelization. Required sample size calculations can be done using the CLT. I can add that in to the answer shortly if you’re interested. $\endgroup$
    – Taylor
    Jul 20, 2018 at 19:17
  • $\begingroup$ ...sure I am. As stated in the question (and other related questions such as stats.stackexchange.com/questions/326320/…, I'm considering the asymptotic case. Is there a polynomial time Monte Carlo algo for arbitary ${\mathbf{\Sigma }}$? If not, is there a polytime algo is some special cases? I'll add this in the question if you don't mind $\endgroup$
    – user193726
    Jul 20, 2018 at 19:28
  • 1
    $\begingroup$ @FabricePautot "How many samples do we need in order to properly estimate $\mathbb{E}(I)$?" We could say that the sample distribution (for the mean of a sample from $I$) can be well approximated by a normal distribution. The standard deviation of the sample distribution will depend on the standard deviation of the distribution of $I$. This standard deviation will be in the worst case equal to $0.5 \cdot n$ (when the distribution of $I$ concentrates on the ends $P(I=1)=(I=n)=0.5$). You can use this to compute the desired necessary sample size for the Monte Carlo simulation. $\endgroup$ Sep 2, 2019 at 11:38
  • 1
    $\begingroup$ @FabricePautot Yes, the distribution of the maximum does not follow a Gaussian distribution. However, the mean of a sample from that distribution does asymptotically approach a Gaussian distribution which will become more narrow as we increase the sample size. $\endgroup$ Sep 3, 2019 at 8:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.