Incorrect results of Mean and Variance of TransformedDistribution

Question

I am solving a problem from 17th Kolmogorov competition. Let $(x,y)$ be a Gaussian random vector with mean $(0,0)$ and a covariance matrix $$\left(\begin{array}{cc} 1&r\\r&1\end{array} \right),$$ where $|r|\le 1$. Find mean and variance of $\max(x,y)$. My try in 13.2 on Windows 10 is as follows.

Mean[TransformedDistribution[Max[x, y], {x, y} \[Distributed] 
BinormalDistribution[{0, 0}, {1, 1}, r]]]

(2 Sqrt[1 - r] + Sqrt[2] r - Sqrt[2] r^2 + Sqrt[2] (-1 + r) Abs[r])/(2 Sqrt[\[Pi]])

Variance[TransformedDistribution[ Max[x, y], {x, y} \[Distributed] 
BinormalDistribution[{0, 0}, {1, 1}, r]]]

Piecewise[{{(-1 + Pi + r)/Pi, Inequality[0, LessEqual, r, Less, 1]}}, (-1 + Pi + r + 2*r^2 - 4*r^3 + 2*r^4)/Pi]

If $r=-1$, then $\max(x,y)=|x|$ and the above results are not in accordance with

Mean[TransformedDistribution[RealAbs[x], x \[Distributed] NormalDistribution[0, 1]]]

Sqrt[2/\[Pi]]

and

Variance[TransformedDistribution[RealAbs[x], x \[Distributed] NormalDistribution[0, 1]]]

(-2 + \[Pi])/\[Pi]

I think it should be r instead of Abs[r] in the result for the mean.

How to obtain the correct results with Mathematica?

PS.Rewriting $\max(x,y)$ in other formulas does not help:

Mean[TransformedDistribution[(RealAbs[x - y] + RealAbs[x + y])/2, 
{x, y} \[Distributed] BinormalDistribution[{0, 0}, {1, 1}, r]]]

(Sqrt[1 - r] - Sqrt[2] r^2 + Sqrt[1 + r] + Sqrt[2] r Abs[r])/Sqrt[\[Pi]]

Mean[TransformedDistribution[ RealAbs[x - y]/2, {x, y} \[Distributed] 
BinormalDistribution[{0, 0}, {1, 1}, r]]] + 
Mean[TransformedDistribution[RealAbs[x + y]/2, {x, y} \[Distributed]BinormalDistribution[{0, 0}, {1, 1}, r]]]

Sqrt[1 - r]/Sqrt[\[Pi]] + Sqrt[1 + r]/Sqrt[\[Pi]]

Answer 1

Without adding Assumptions (on Windows 10, Mathematica 13.2.0.0) I get the same "wrong" mean (at least wrong for $-1 \leq r < 0$). From @BobHanlon 's answer, the answer depends on the operating system. However, on Windows 10, Mathematica 13.2.0.0 one can get the right formulas for the mean and variance for $-1<r<0$ by adding in that assumption in TransformedDistribution:

Mean[TransformedDistribution[Max[x, y], 
  {x, y} \[Distributed] BinormalDistribution[{0, 0}, {1, 1}, r], 
  Assumptions -> -1 < r < 0]]
(* Sqrt[1 - r]/Sqrt[π] *)

Variance[TransformedDistribution[Max[x, y], 
  {x, y} \[Distributed] BinormalDistribution[{0, 0}, {1, 1}, r], 
  Assumptions -> -1 < r < 0]]
(* (-1 + π + r)/π *)

The rest of this is just an extended comment providing support (rather than "verification") that for the formula for the mean, Abs[r] should be r and the desired function for the variance.

When one doesn't have an exact formula or is unsure of a formula, simulations can be helpful. Below a million samples are taken for several values of r and the sample mean determined. That is plotted against the two formulas for the mean.

t = Table[{r, Mean[Max[#] & /@ RandomVariate[BinormalDistribution[{0, 0}, {1, 1}, r], 1000000]]},
  {r, -0.99, 0.99, 0.02}];

Show[Plot[{(2*Sqrt[1 - r] - Sqrt[2]*r^2 + Sqrt[2]*r*Abs[r])/(2*Sqrt[Pi]), 
  Sqrt[1 - r]/Sqrt[Pi]}, {\[Rho], -1, 1}, 
  PlotRange -> {All, {0, 0.8}}, PlotStyle -> {{LightGray, Thickness[0.015]}, {Cyan, Thickness[0.015]}}, 
  PlotLegends -> LineLegend[{Red, LightGray, Cyan}, {"Simulations", 
  "\!\(\*FractionBox[\(2\\\ \*SqrtBox[\(1 - r\)] - \*SqrtBox[\(2\)]\\\ \
\*SuperscriptBox[\(r\), \(2\)] + \*SqrtBox[\(2\)]\\\ r\\\ Abs[r]\), \
\(2\\\ \*SqrtBox[\(\[Pi]\)]\)]\)", 
    "\!\(\*FractionBox[\(2\\\ \*SqrtBox[\(1 - r\)] - \*SqrtBox[\(2\)]\
\\\ \*SuperscriptBox[\(r\), \(2\)] + \*SqrtBox[\(2\)]\\\ \
\*SuperscriptBox[\(r\), \(2\)]\), \(2\\\ \*SqrtBox[\(\[Pi]\)]\)]\)=\!\
\(\*FractionBox[SqrtBox[\(1 - r\)], SqrtBox[\(\[Pi]\)]]\)"}]],
 ListPlot[t, PlotStyle -> Red, Joined -> True]]

For the variance the formula Mathematica produces is also wrong for $-1\leq r <0$:

v = Table[{r, Variance[Max[#] & /@ RandomVariate[BinormalDistribution[{0, 0}, {1, 1}, r], 1000000]]},
   {r, -0.99, 0.99, 0.02}];

Show[Plot[{Piecewise[{{(-1 + Pi + r)/Pi, Inequality[0, LessEqual, r, Less, 1]}}, 
  (-1 + Pi + r + 2*r^2 - 4*r^3 + 2*r^4)/Pi], (-1 + Pi + r)/Pi}, {r, -1, 1}, 
  PlotRange -> {{-1, 1}, {0, 1.8}}, 
  PlotStyle -> {{Thickness[0.015], LightGray}, {Thickness[0.015], Cyan}},
  PlotLegends -> LineLegend[{Red, LightGray, Cyan}, 
  {"Simulations", "Piecewise function", "(-1+Pi+r)/Pi"}]], 
 ListPlot[v, PlotStyle -> Red, Joined -> True]]

Answer 2

$Version

(* "13.2.1 for Mac OS X ARM (64-bit) (January 27, 2023)" *)

Clear["Global`*"]

dist = BinormalDistribution[r];

DistributionParameterAssumptions[dist]

(* -1 < r < 1 *)

dist2 = TransformedDistribution[Max[x, y], {x, y} \[Distributed] dist];

μ = Mean[dist2]

(* Sqrt[1 - r]/Sqrt[π] *)

Limit[μ, r -> #] & /@ {-1, 1}

(* {Sqrt[2/π], 0} *)

var = Variance[dist2]

(* (-1 + π + r)/π *)

Limit[var, r -> #] & /@ {-1, 1}

(* {(-2 + π)/π, 1} *)

Answer 3

There is a workaround in 13.2 on Windows 10. Making use of the identity (In my PS this was written incorrectly.) $$ \max(x,y)= \frac {|x-y|+x+y} 2$$ and the property of mean $E(u+v)=E(u)+E(v)$ which is valid for all random variables $u,v$ (of course, s.t. at least two of these means exist), one obtains

Mean[TransformedDistribution[(RealAbs[x - y])/2, {x, y} \[Distributed] 
BinormalDistribution[{0, 0}, {1, 1}, r]]]

Sqrt[1 - r]/Sqrt[\[Pi]]

Mean[TransformedDistribution[(x + y)/2, {x, y} \[Distributed] 
BinormalDistribution[{0, 0}, {1, 1}, r]]]

0

In order to find the variance of $\max(x,y)$, one should use the independence of $x+y$ and $x-y$ which follows from

Mean[TransformedDistribution[(x - y)*(x + y), {x, y} \[Distributed] 
BinormalDistribution[{0, 0}, {1, 1}, r]]]

0

and

Mean[TransformedDistribution[x + y, {x, y} \[Distributed] 
BinormalDistribution[{0, 0}, {1, 1}, r]]]

0

and

Mean[TransformedDistribution[x - y, {x, y} \[Distributed] 
BinormalDistribution[{0, 0}, {1, 1}, r]]]

0

I leave the rest on your own.

Answer 4

We are given $(X,Y) \sim$ standard Bivariate Normal, with correlation coefficient $r$ and pdf $f(x,y)$:

Then $E[\text{max}(X,Y)]$ is given immediately as:

and $\text{Var}[\text{max}(X,Y)]$ is:

where I am using the Expect and Var functions from the mathStatica package for Mathematica. The salient point is that it should not be necessary to first find the transformed distribution, and to then find the moments of the latter: we can find the moments of the transformed variable directly wrt to the parent standard bivariate Normal pdf.