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