Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It only takes a minute to sign up.
Sign up to join this community
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
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:
Ker, Alan P., On the maximum of bivariate normal random variables, Extremes 4, No. 2, 185–190 (2001). ZBL1003.60017
Aksomaitis, A.; Burauskaitė-Harju, A., The moments of the maximum of normally distributed dependent values, Information Technology and Control 38, No. 4, 301–302 (2009)
Ross, Andrew M., Useful bounds on the expected maximum of correlated normal variables, ISE Working Paper 03W-004 (2003)
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.
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
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
Link: https://stsda.kaust.edu.sa/Documents/2008.AG.SPL.pdf
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.
