# Distribution of the maximum of a multivariate normal random variable

Suppose there is a vector of jointly normally distributed random variables $X \sim \mathcal{N}(\mu_X, \Sigma_X)$. What is the distribution of the maximum among them? In other words, I am interested in this probability $P(max(X_i) < x), \forall i$.

Thank you.

Regards, Ivan

• Have you tried to solve the 2-dimension case? May 15 '12 at 14:36
• this is the CDF of $X$, there is no analytic solution here, but there is a name for it gaussian coupla. Apr 24 '20 at 0:46
• In case someone is interested, here it is the same question but regarding the minimum: stats.stackexchange.com/questions/518032/… Apr 3 at 12:18

For general $(\mu_X,\Sigma_X)$ The problem is quite difficult, even in 2D. Clearly $P(max(X_i)<x) = P(X_1<x \wedge X_2<x \cdots \wedge X_d <x)$, so to get the distribution function of the maximum one must integrate the joint density over that region... but that's not easy for a general gaussian, even in two dimensions. I'd bet there is no simple expression for the general case.

Some references:

http://itc.ktu.lt/itc384/Aksom384.pdf

http://people.emich.edu/aross15/q/papers/bounds_Emax.pdf

• Thanks for the answer, @leonbloy. It is a little bit disappointing to hear, I thought everything was relatively simple with the Gaussian distribution.
– Ivan
May 15 '12 at 18:27
• And how about numerical estimation of the expected value and variance of the maximum of jointly normal random variables?
– Ivan
May 28 '12 at 13:56
• @Ivan: well, the estimation seems quite straightforward to me: produce some sample values and use the common sample estimators of mean and variance. are you thinking of something else? May 28 '12 at 17:21

Multivariate Skew-Normal Distributions and their Extremal Properties

Rolf Waeber February 8, 2008 Abstract In this thesis it is established that the distribution is a skew normal dist.

A paper by Nadarajah and Samuel Kotz gives the expression for the max of any bivariate normal F(x,y). IEEE TRANSACTIONS ON VERY LARGE SCALE INTEGRATION (VLSI) SYSTEMS, VOL. 16, NO. 2, FEBRUARY 2008 Exact Distribution of the Max/Min of Two Gaussian Random Variables Saralees Nadarajah and Samuel Kotz If F(x,y) is a standard normal (means=0 and variances=1, r>0) the dist of the maximum is a skew normal.

• While this link may answer the question, it is better to include the essential parts of the answer here and provide the link for reference. Link-only answers can become invalid if the linked page changes. Oct 26 '12 at 12:05
• @rschwieb: What link? There is no link. Oct 26 '12 at 12:26
• @AsafKaragila Explain to me why a naked citation with no explanation is different from a link with no explanation. You know the options are limited in the menu. Oct 26 '12 at 12:39
• @rschwieb: Because the idea behind a citation is to give a reference that can be found later on; whereas links can indeed rot. Regardless to the status of a link to fundamenta mathematicae archives, the reference will not rot at least as long as there are libraries around the globe. You can always add a comment on your own and tell the software not to use the default comments at all. Oct 26 '12 at 12:41
• @AsafKaragila OK, I agree, that is one difference. However something else in the menu description (which I cannot find right now without any review tasks available) was my basis. I thought it said something to the effect of "while the (link) itself may answer the question, you really should add more content." I will have to check when it is available to me again. I think it is pretty clear that bare citations are poor quality answers, especially when the OP was not a reference request. Yes there are libraries, but not everyone can access the material. Oct 26 '12 at 12:48

Please see paper by Reinaldo B. Arellano-Vallea and Marc G. Genton:

On the exact distribution of the maximum of absolutely continuous dependent random variables

Corollary 4 (page 31) gives the general form for the distribution of the maximal of a multivariate Gaussian. Discussion following the corollary says that the distribution of maximal is skew-normal when $$X$$ is bivariate normal.

Let us take $$d=2$$ and assume for simplicity that $$\mu_X=(0,0)$$. Denote $$\Sigma_X:=((\sigma_1^2,\rho),(\rho,\sigma_2^2))$$ where $$\sigma_{1,2} \ge 0$$ and $$-\sigma_1 \sigma_2 \le \rho \le \sigma_1 \sigma_2$$. Furthermore denote $$\delta:=-\rho^2 + \sigma_1^2 \sigma_2^2$$. Now the result is as follows: $$\begin{eqnarray} &&\sqrt{\frac{2 \pi}{\delta}}{\mathbb P}\left( max(X_1,X_2) \le X \right)=\\ && \frac{\sqrt{2 \pi } \left(T\left[\frac{X}{\sigma_1},\frac{\rho -\sigma_1^2}{\sqrt{\delta}}\right]+T\left[\frac{X}{\sigma_2},\frac{\rho -\sigma_2^2}{\sqrt{\delta}}\right]\right)}{\sqrt{\delta}}+\frac{\sqrt{\frac{\pi }{2}} \left(\text{erf}\left(\frac{X}{\sqrt{2} \sigma_1}\right)+\text{erf}\left(\frac{X}{\sqrt{2} \sigma_2}\right)+2\right)}{2 \sqrt{\delta}} \end{eqnarray}$$

Here $$T[\cdot,\cdot]$$ is the Owen's T function. In deriving the result we used Generalized Owen's T function.

In[894]:= (*The generalized Owen T function.*)
(*NIntegrate[Exp[-x^2/2]/Sqrt[2Pi]1/2 Erf[(a \
x+b)/Sqrt[2]],{x,h,Infinity},WorkingPrecision\[Rule]20]*)
Clear[T];
T[h_, a_, b_] :=
1/(2 Pi) (ArcTan[a] - ArcTan[a + b/h] - ArcTan[a + h/b + a^2 h/b]) +
1/4 Erf[b/Sqrt[2 (1 + a^2)]] + OwenT[h, a + b/h] +
OwenT[b/Sqrt[1 + a^2], a + h/b + a^2 h/b];

{s1, s2} = RandomReal[{0, 2}, 2, WorkingPrecision -> 50];
rho = RandomReal[{-s1 s2, s1 s2}, WorkingPrecision -> 50];
CC = {{s1^2, rho}, {rho, s2^2}};
dd = -rho^2 + s1^2 s2^2;
CCI = 1/dd {{s2^2, -rho}, {-rho, s1^2}}; d = 2; x =.; xx =
Table[x[p], {p, 1, d}];
X = RandomReal[{-1, 1}, WorkingPrecision -> 50];

Sqrt[2 Pi]/Sqrt[Det[CC]] NIntegrate[
Exp[-1/2 xx.(CCI.xx)]/Sqrt[(2 Pi)^2 Det[CC]], {x[1], -Infinity,
X}, {x[2], -Infinity, X}]
(Sqrt[\[Pi]/2] (1 + Erf[X/(Sqrt[2] s2)]))/(
2 Sqrt[-rho^2 + s1^2 s2^2]) +
Sqrt[2 Pi]/Sqrt[-rho^2 + s1^2 s2^2]
NIntegrate[(
E^(-(x[2]^2/
2 )) (Erf[(s2 X + rho x[2])/(
Sqrt[2] Sqrt[- rho^2 + s1^2 s2^2])]))/(
Sqrt[2 Pi] 2 ), {x[2], -X/s2, Infinity}]
(Sqrt[\[Pi]/2] (1 + Erf[X/(Sqrt[2] s2)]))/(
2 Sqrt[-rho^2 + s1^2 s2^2]) +
Sqrt[2 Pi]/Sqrt[-rho^2 + s1^2 s2^2]
T[-X/s2, rho /Sqrt[- rho^2 + s1^2 s2^2], (s2 X )/
Sqrt[- rho^2 + s1^2 s2^2]]
Sqrt[\[Pi]/2]  (2 + Erf[X/(Sqrt[2] s1)] + Erf[X/(Sqrt[2] s2)])/(
2 Sqrt[dd]) +
Sqrt[2 Pi]/Sqrt[
dd] (OwenT[X/s1, (rho - s1^2)/Sqrt[dd]] +
OwenT[X/s2, (rho - s2^2)/Sqrt[dd]])

Out[901]= 0.0356087

Out[902]= 0.0356087

Out[903]= 0.03560872573422025123907864650438155988048895417

Out[904]= 0.03560872573422025123907864650438155988048895417