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

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

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