# Compute moments of maximum of multivariate normal distribution

I want to compute the first and second order moment of the maximum of a random vector from multivariate normal distribution, i.e., compute

$E[\,\,Y\,],\,\,E[\,Y^2\,],$

where

$Y=\max (X),\,\,\\ X=(x_1,x_2,...,x_k) \sim \mathcal{N}(\mu,\Sigma)\\$

$\mu,\Sigma$ are mean vector (k$\times$1) and covariance matrix (k$\times$k), respectively.

The max operation means: $Y$ equals to the max element of vector $X$.

No independency of $x_1,x_2,...,x_k$ are imposed on $\Sigma$, i.e., $\Sigma$ is a general symmetric covariance matrix.

Any hints on how to solve this problem analytically? Approximate method will also be appreciated, e.g., Monte Carlo or Numerical integration, etc.

• A Monte Carlo resolution is obvious. What is the question? – Xi'an Nov 2 '17 at 11:39
• @Xi'an, Hi, Prof., yes, we can simply draw samples of X then compute empirical moments of Y. I tried this method, however it is slow due to the large sample size of X, so I want to seek some efficient methods (more efficient MC method or numerical integration method) – lynnjohn Nov 2 '17 at 11:53
• The distribution of the maximum may be a skew-normal, as discussed in that thread. – Xi'an Nov 2 '17 at 11:57
• @Xi'an, thank u, I see that post, but is there a way to efficiently estimate the first and second moments of maximum of multivariate Gaussian? – lynnjohn Nov 2 '17 at 12:05

And Arenallo-Valle and Genton (2008) have this representation for the pdf of the Normal maximum, $X_{(N)}$:
which involves $\phi_1$, the marginal pdf of $X_n$ and $\Phi_{n-1}$, the multivariate Normal cumulative distribution function of the (n-1) dimensional Normal. This cdf is available in the R package mvtnorm. The above expression simplifies when the components $X_i$ are exchangeable, i.e. when $\mu$ and $\Sigma$ are invariant by permutation.