# 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)". – Roman May 22 at 18:18 • Windows 10, Mathematica 11.3: I also get 0.886433. – JimB May 22 at 18:26 • @Roman I don't understand why you get a different expression. – chris May 22 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 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. – Roman May 22 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. – chris May 23 at 5:47