# How to use Mathematica to find the maximum expected value of multiple normal distributions?

Is it possible to calculate the maximum expected value of $$n (n\geq 2)$$ independent (identical distributed) random variables that follow a normal distribution? If possible, how to find it? If it cannot be solved, is there an approximate solution method?

Related post:

Below is my code.

But it takes too long time to calculate meanValues, are there some mistakes in my code? Any help would be appreciated.

Clear["Global*"];

MultiIIDGaussianMaximumDistribution[n_Integer?Positive, \[Mu]_, \[Sigma]_] :=
Module[{xSymbols, condition, dist},
If[n < 2, Message[MultiIIDGaussianMaximumDistribution::arg, n];
Return[\$Failed]];
xSymbols = Table[Subscript[x, i], {i, 1, n}];
condition =
Table[xSymbol \[Distributed]
NormalDistribution[\[Mu], \[Sigma]], {xSymbol, xSymbols}];
dist = TransformedDistribution[Max @@ xSymbols, condition];
Return[dist];
]

pdfs = Table[
PDF[MultiIIDGaussianMaximumDistribution[n, \[Mu], \[Sigma]], y],
{n, 2, 7}
]

Plot[
Evaluate[
Table[pdfs[[n - 2]] /. {\[Mu] -> 0, \[Sigma] -> 1}, {n, 2, 7}]
],
{y, -10, 10},
PlotRange -> Full,
PlotLegends -> Table["n = " <> ToString[n], {n, 2, 7}],
PlotLabel -> "PDF of Maximum of n IID Gaussian Distributions",
PlotStyle -> Table[ColorData[1][n - 2], {n, 2, 7}]
]

meanValues = Table[
Mean[MultiIIDGaussianMaximumDistribution[n, 0, 1]],
{n, 2, 7}
]


$$\frac{e^{-\frac{(y-\mu )^2}{2 \sigma ^2}} \text{erfc}\left(\frac{\mu -y}{\sqrt{2} \sigma }\right)}{\sqrt{2 \pi } \sigma },\frac{3 e^{-\frac{(y-\mu )^2}{2 \sigma ^2}} \text{erfc}\left(\frac{\mu -y}{\sqrt{2} \sigma }\right)^2}{4 \sqrt{2 \pi } \sigma },\frac{e^{-\frac{(y-\mu )^2}{2 \sigma ^2}} \text{erfc}\left(\frac{\mu -y}{\sqrt{2} \sigma }\right)^3}{2 \sqrt{2 \pi } \sigma },\frac{5 e^{-\frac{(y-\mu )^2}{2 \sigma ^2}} \text{erfc}\left(\frac{\mu -y}{\sqrt{2} \sigma }\right)^4}{16 \sqrt{2 \pi } \sigma },\frac{3 e^{-\frac{(y-\mu )^2}{2 \sigma ^2}} \text{erfc}\left(\frac{\mu -y}{\sqrt{2} \sigma }\right)^5}{16 \sqrt{2 \pi } \sigma },\frac{7 e^{-\frac{(y-\mu )^2}{2 \sigma ^2}} \text{erfc}\left(\frac{\mu -y}{\sqrt{2} \sigma }\right)^6}{64 \sqrt{2 \pi } \sigma }$$

(* Monte Carlo simulation for mean values *)
meanValuesMonteCarlo[n_, \[Mu]_, \[Sigma]_, numSamples_] :=
Module[{samples, maxSamples},
samples = RandomVariate[NormalDistribution[\[Mu], \[Sigma]], {numSamples, n}];
maxSamples = Max /@ samples;
Mean[maxSamples]
]

(* Parameters for Monte Carlo simulation *)
numSamples = 10^6;
meanValues = Table[
meanValuesMonteCarlo[n, 0, 1, numSamples],
{n, 2, 7}
]


Monte Carlo simulation results are:

$$\{0.563841,0.84589,1.02998,1.16336,1.26683,1.35187\}$$

PDFs of these synthesised distributions:

• Is there a reason for you to suspect that there is an explicit solution? For example, the answers in your linked MSE post only give the approximations, which you can validate against a numerical integration: Table[{n, MultiIIDGaussianMaximumDistribution[n, 0, 1]]}, {n, 1, 3}] gives {{1, 0.}, {2, 0.56419}, {3, 0.846284}}. Commented May 16 at 15:03

This is only a partial answer in that it only provides the means for $$n=2$$ through $$n=5$$. With $$n>6$$ Mathematica only returns the input.

means = TableForm[
Table[{n, Mean[OrderDistribution[{NormalDistribution[μ, σ], n}, n]] // FullSimplify},
{n, 2, 5}], TableHeadings -> {None, {"n", "Mean"}}]


For values of $$n>5$$ (as @Domen suggests) use numerical integration:

TableForm[
Table[{n, μ + σ  NIntegrate[x
PDF[OrderDistribution[{NormalDistribution[0, 1], n}, n], x] // N,
{x, -∞, ∞}]}, {n, 6, 20}], TableHeadings -> {None, {"n", "Mean"}}]


• Oooops! Somehow I deleted PDF from the code. Will fix.
– JimB
Commented May 16 at 15:40
• Sorry. Now fixed.
– JimB
Commented May 16 at 15:44
• μ + σ NIntegrate[x PDF[OrderDistribution[{NormalDistribution[0, 1], n}, n], x] // N uses the properties of expectation operator. Thank you for your answer! Another question: Can we use mma for demonstration of asymptotic analysis? Commented May 16 at 15:48
• ExtremeValueDistribution Applications section would help? Commented May 16 at 16:00
• I don't see how the ExtremeValueDistribution helps with your specific question. You could plot the density of the ExtremeValueDistribution with the same mean and variance but I'm not sure what that tells you for small sample sizes. But certainly whatever "asymptotic analysis" means, Mathematica has many related functions: Limit, AsymptoticSum`, etc.
– JimB
Commented May 17 at 0:20