# Unstable symbolic expectation of gaussian distribution on absolute values?

Context

I use mathematica 11.3 and 12.

11.3.0 for Mac OS X x86 (64-bit) (March 7, 2018)

12.0.0 for Mac OS X x86 (64-bit) (April 7, 2019)

I am interested in evaluating the following integral

  exp = Expectation[
Abs[x] Abs[y], {x, y} \[Distributed]MultinormalDistribution[
{0, 0}, {{a, c}, {c, b}}]]// FullSimplify[#, c > 0] &;


(* (Sqrt[a b-c^2] (2 a b-c^2)+a b c cot^-1(Sqrt[a b-c^2]/c))/([Pi] a b) *)

If I then evaluate the same integral but while giving a,b,c values first I get

  expN =  Expectation[Abs[x] Abs[y], {x, y} \[Distributed]
MultinormalDistribution[{0, 0}, {{1, 1/4}, {1/4, 2}}]] // N


(* 0.914421 *)

But if I take the symbolic result and evaluate it with those same values I get

 exp /. {a -> 1, b -> 2, c -> 1/4} // N


(* 0.886433 *)

Question

Can anyone reproduce what seems to be a bug?

Can you suggest any workaround?

Update

In view of the comments below I would like to label this issue as a bug. The fact that mathematica yields different answers to a given analytical integral as a function of evaluation time qualifies as a bug IMHO.

Any objection?

• I cannot reproduce the problem. On my machine, exp /. {a -> 1, b -> 2, c -> 1/4} // N also produces 0.914421. Maybe because exp is actually (2 (Sqrt[a b - c^2] + c ArcTan[c/Sqrt[a b - c^2]]))/\[Pi] instead of what you write. Using $Version == "12.0.0 for Mac OS X x86 (64-bit) (April 7, 2019)". May 22, 2019 at 18:18 • Windows 10, Mathematica 11.3: I also get 0.886433. – JimB May 22, 2019 at 18:26 • @Roman I don't understand why you get a different expression. May 22, 2019 at 18:35 • But if one uses exp = Expectation[ Abs[x] Abs[y], {x, y} \[Distributed] MultinormalDistribution[{0, 0}, {{a, c}, {c, b}}], Assumptions -> {a > 0, b > 0, c > 0}] // FullSimplify, then the substitution works fine and one gets a different symbolic result:$\frac{4 \sqrt{a b-c^2}+2 c \tan ^{-1}\left(\frac{c}{\sqrt{a b-c^2}}\right)-2 c \cot ^{-1}\left(\frac{c}{\sqrt{a b-c^2}}\right)+\pi c}{2 \pi }$– JimB May 22, 2019 at 18:37 • I revise my post and now agree with @MichaelE2. Delayed-defining exp := Expectation[Abs[x] Abs[y], {x, y} \[Distributed] MultinormalDistribution[{0, 0}, {{a, c}, {c, b}}]] and then executing A = Table[exp, {10}]; I get seven different answers out of ten attempts (Length[DeleteDuplicates[A]] gives 7), some of which give a numerical 0.914421 and some a 0.886433. May 22, 2019 at 19:43 ## 1 Answer $Version

(* "12.0.0 for Mac OS X x86 (64-bit) (April 7, 2019)" *)


Checking by comparing the results for the distribution defined with MultinormalDistribution and BinormalDistribution

Clear["Global*"]

distMN = MultinormalDistribution[{0, 0}, {{a, c}, {c, b}}];


The implied assumptions (i.e., required for a valid distribution) are

assumeMN = DistributionParameterAssumptions[distMN]

(* a > 0 && b > 0 && c ∈ Reals && a b - c^2 > 0 *)


Note that this is different than assuming c > 0.

Similarly,

distBN = BinormalDistribution[{0, 0}, {Sqrt[a], Sqrt[b]}, c/Sqrt[a*b]];

assumeBN = DistributionParameterAssumptions[distBN]

(* Sqrt[a] > 0 && Sqrt[b] > 0 && -1 < c/Sqrt[a b] < 1 *)


Verifying that the distributions are equivalent under either form of the assumptions:

Simplify[
PDF[distMN, {x, y}] == PDF[distBN, {x, y}],
assumeMN] &&
Simplify[
PDF[distMN, {x, y}] == PDF[distBN, {x, y}],
assumeBN]

(* True *)

expMNgen = Assuming[assumeMN,
Expectation[Abs[x] Abs[y], {x, y} \[Distributed] distMN] // FullSimplify]


Assuming c > 0 (i.e., positive correlation coefficient) with MultinormalDistribution

expMN = expMNgen // FullSimplify[#, c > 0] &


params = {a -> 1, b -> 2, c -> 1/4};

(expMN /. params) // N

(* 0.914421 *)

expMNN = Expectation[Abs[x] Abs[y],
{x, y} \[Distributed] (distMN /. params)] // N

(* 0.914421 *)

expBNgen = Assuming[assumeBN,
Expectation[Abs[x] Abs[y], {x, y} \[Distributed] distBN] // FullSimplify]


Despite the different forms, expMNgen and expBNgen are equal

expMNgen == expBNgen // FullSimplify

(* True *)


Assuming c > 0 with BinormalDistribution

expBN = expBNgen // Simplify[#, c > 0] &


As expected, expMN and expBN are equal for c > 0

FullSimplify[expBN == expMN, c > 0]

(* True *)

expBN /. params // N

(* 0.914421 *)

expBNN = Expectation[Abs[x] Abs[y],
{x, y} \[Distributed] (distBN /. params)] // N

(* 0.914421 *)
`
• Thank you for your answer. Have you noticed any instability in the answer when you request it more than once? It seems this is the most serious concern here. May 23, 2019 at 5:47