Probability that two correlated Normal random variables have the same sign?

Questions
1

When computing the probability for two correlated (zero mean) Normal random variables, I’m getting conflicting results. One of them seems nonsensical and potentially indicative of a bug in Mathematica.

In the attached notebook, I define a MultiNormalDistribution and ask Mathematica to compute the probability that x*y > 0. I believe this should lead to the same result as the probability that (x > 0 && y > 0) || (x < 0 && y < 0). The latter result is quite sensible whereas the probability that the product is positive is very odd.

Does this make sense?

SignAgreementBug.nb (12.5 KB)

2

It’s more than odd. The two approaches you give should give the same result but only the second one below does give the right answer. (You should inform Wolfram, Inc.)

FullSimplify[Probability[x y >= 0, {x, y} \[Distributed] MultinormalDistribution[{{1, \[Rho]}, {\[Rho], 1}}]],
   Assumptions -> -1 <= \[Rho] <= 1]
(* 1/2 *)

2 Probability[x >= 0 && y >= 0, {x, y} \[Distributed] MultinormalDistribution[{{1, \[Rho]}, {\[Rho], 1}}]] // FullSimplify
(* 1/2 + ArcSin[\[Rho]]/\[Pi]  *)

An alternative way that does work with x*y is the following:

distxy =  TransformedDistribution[x y, {x, y} \[Distributed]  MultinormalDistribution[{{1, \[Rho]}, {\[Rho], 1}}]]]
1 - CDF[distxy, 0] // FullSimplify
(* 1/2 + ArcSin[\[Rho]]/\[Pi] *)
3

Thanks Jim! Will do.

4

Thank you for catching the discrepancy - it looks like it is an issue with Probability, I will investigate.

5

Hi Gosia Konwerska,

I just purchased a new license for Mathematica 14.0.0 and observe the same issue is still present. Has there been any progress on this?

Thanks,
Greg Sternberg

6

Unfortunately has not been resolved yet…

7

It seems to me that the problem lies with Integrate with parameters, rather than with Probability. In the following examples the result depends on how we specify the integration region and on the assumptions on the parameter:

dist = MultinormalDistribution[{{1, \[Rho]}, {\[Rho], 1}}];
pdf = PDF[dist, {x, y}];
reg = ImplicitRegion[a >= 0 && b >= 0, {a, b}];
int1 = Integrate[pdf,
  {x, y} \[Element] reg,
  Assumptions -> -1 < \[Rho] < 1]
int2 = Integrate[pdf,
  {x, y} \[Element] reg,
  Assumptions -> -1 < \[Rho] < 0]
int3 = Integrate[pdf,
  {x, 0, Infinity}, {y, 0, Infinity},
  Assumptions -> -1 < \[Rho] < 1]
int4 = Integrate[pdf,
  {x, 0, Infinity}, {y, 0, Infinity},
  Assumptions -> -1 < \[Rho] < 0]
{int1, int2, int3, int4} /. \[Rho] -> -1/2