2
$\begingroup$

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

Mathematica graphics

(* (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?

$\endgroup$
  • 1
    $\begingroup$ 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)". $\endgroup$ – Roman May 22 at 18:18
  • 1
    $\begingroup$ Windows 10, Mathematica 11.3: I also get 0.886433. $\endgroup$ – JimB May 22 at 18:26
  • $\begingroup$ @Roman I don't understand why you get a different expression. $\endgroup$ – chris May 22 at 18:35
  • $\begingroup$ 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 }$ $\endgroup$ – JimB May 22 at 18:37
  • 2
    $\begingroup$ 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. $\endgroup$ – Roman May 22 at 19:43
2
$\begingroup$
$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]

enter image description here

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

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

enter image description here

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]

enter image description here

Despite the different forms, expMNgen and expBNgen are equal

expMNgen == expBNgen // FullSimplify

(* True *)

Assuming c > 0 with BinormalDistribution

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

enter image description here

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 *)
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – chris May 23 at 5:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.