# Incorrect results of Mean and Variance of TransformedDistribution

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]]

• For Abs[r]==1 the covariance matrix {{1, r}, {r, 1}} is not invertible. How is the Gaussian defined then? Apr 18 at 8:08
• @HenrikSchumacher: Thank you for your interest to the question. Sorry for the delay in my reply: force majeure. First, the case r==-1 should be considered as the limit when r tends to - 1 from above. BTW, OK with r==1. Second, there are problems with other negative values of r. For example, (2 Sqrt[1 - r] + Sqrt[2] r - Sqrt[2] r^2 + Sqrt[2] (-1 + r) Abs[r])/( 2 Sqrt[\[Pi]]) /. r -> -0.5 produces 0.0925749, whereas Mean[TransformedDistribution[ Max[x, y], {x, y} \[Distributed] BinormalDistribution[{0, 0}, {1, 1}, -0.5]]] results in 0.690988. Apr 18 at 13:32

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 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]]


• Thank you for your work and +1. Indeed, the true answers are Sqrt[1 - r]/Sqrt[π] and (-1 + π + r)/π. The ones may be derived by making use of the identity Max[x,t]===RealAbs[x - y] + RealAbs[x + y])/2. Apr 19 at 1:00
• @user64494 Just added in the use of Assumptions but that doesn't fix the somewhat pathological case of $r=-1$ as TransformedDistribution complains if one sets r to -1.
– JimB
Apr 19 at 1:03
• In my comment I have had in mind the proof by hand. Apr 19 at 9:01
• JimB (@ does not work): It is not a good practice to substantially edit a post without any indication of changes. Hope I am clear. Apr 19 at 23:41
• @user64494 "substantial" is in the eye of the beholder. I placed adequate text about the edits in the "Edit summary" (which I think were better than some that I've seen before such as "added 586 characters in body").
– JimB
Apr 20 at 0:07
\$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} *)

• Thank you for your work and +1. This is basically that I do, but on Mac. On Windows 10 I still obtain (2 Sqrt[1 - r] - Sqrt[2] r^2 + Sqrt[2] r Abs[r])/(2 Sqrt[\[Pi]]) for the mean. Apr 18 at 14:16
• @user64494, I can confirm that "12.3.1 for Microsoft Windows (64-bit) (June 24, 2021)" , "13.0.0 for Microsoft Windows (64-bit) (December 3, 2021)" and 13.2.0 for Linux x86 (64-bit) (December 12, 2022) (Wolfram Cloud) all return the mean with Abs[r]. Apr 18 at 14:25
• @Domen: Thank you. Apr 18 at 14:34
• On an older laptop with 13.2.1 for Mac OS X x86 (64-bit) (January 27, 2023) I get an incorrect result initially, but repeated evaluations are somewhat quicker and get to the correct result. I would guess that the initial evaluation is timing out before it comes to the correct result and an unverified intermediate result is returned. However, subsequent calculations have access to cached intermediate results and get to the correct solution in time. Apr 18 at 14:56
• @BobHanlon, I can confirm your guess on "13.2.1 for Microsoft Windows (64-bit) (January 27, 2023)". After the 3rd evaluation of the Mean, I get the same result as you did. Apr 18 at 15:08

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.

• The rest is Variance[ TransformedDistribution[(x + y)/2, {x, y} \[Distributed] BinormalDistribution[{0, 0}, {1, 1}, r]]]+Mean[TransformedDistribution[((x - y)/2)^2, {x, y} \[Distributed] BinormalDistribution[{0, 0}, {1, 1}, r]]]-(Sqrt[1 - r]/Sqrt[\[Pi]])^2. Apr 19 at 16:31
• Do you think it is a bug ultimately? Apr 20 at 2:23

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.

• Can you kindly give a link to the mathStatica package? I don't find it at this forum. In other case this is empty talk. Apr 24 at 15:05
• There are more than 60 articles on this forum on mathStatica. Just do a search. Apr 24 at 15:39
• Sorry, I repeat I don't find a link to download mathStatica at this forum as well as a link with its description. Apr 24 at 15:41
• Sorry, I did not understand. Is there a link to download Mathematica on this forum? If you can provide that, many people might find that helpful. Apr 24 at 16:05
• Can you give any link with the description of mathStatica and any link to download the one? In other case this is empty talk. Apr 24 at 16:45