Bug in the exact 2d/3d orthant probabilities

[vittorioerba]

Hi! First of all, thanks for this awesome package. It has been extremely crucial for a research project of mine.

I think I may have identified a bug in the implementation. In the file functions.jl, at line 159, there is a code that handles 2d and 3d special cases, namely the cases when the integration range is [0,∞] for all variables. In my understanding, these special cases require also the mean to be zero, which is not checked in line 165 or before.

One can check that the current implementation is incorrect using the following code (check 1 and check 2), showing that mvnormcdf is not depedent on the value of the mean if a=0 and b=Inf, and that it does not agree with HCubature.

using MvNormalCDF, HCubature

gausspdf(x, mu, sigma) = exp(-0.5 * dot(x-mu, inv(sigma), x-mu)) / sqrt((2pi)^length(x) * det(sigma))

mu = [1, 1]
sigma = [1.0 0.5; 0.5 1.0]
a = [0.0, 0.0]
b = [Inf, Inf]

# check 1: agains hcubature
mvnormcdf(mu, sigma, a, b)[1] |> println
hcubature( x -> gausspdf(x, mu, sigma), [0., 0.], [1000,1000] )[1] |> println

# check 2: the integral is not dependent on mu
mvnormcdf([0,0], sigma, a, b)[1] |> println
mvnormcdf(mu, sigma, a, b)[1] |> println

# check 3: axis inversion not taken into account
mvnormcdf([0,0], [1.0 -0.5; -0.5 1.0], [0., -Inf], [Inf, 0]) |> println
mvnormcdf([0,0], [1.0 0.5; 0.5 1.0], [0., 0.], [Inf, Inf]) |> println

A possible bug fix is to set a = a - mu and b = b - mu before checking for the special cases, which currently happens only after at line 182/185.

Also, I think the current implementation will not consider the integral with mu=[0,0], a=[0, -Inf] and b=[Inf, 0] as a special case, as check 3 in the code above shows (see the difference in the error estimate). But a coordinate transformation which inverts the second axis suffices to recognise this as a special case again. A possible improvement would be to invert all axis whose boundaries contain exactly one infinity such that the infinity is the upper bound (i.e. (a,b) = (1, 3) would not get inverted, (-Inf, Inf) = (1, 3) would not get inverted, (a,b) = (1, Inf) would not get inverted and (a,b) = (-Inf, 1) would get inverted). I would do this after shifting the boundaries to eliminate the mean. But I had not time to carefully check this proposal, so take this just as a random thought.

Thanks again for publishing this package, I hope my comment is useful :)

[PharmCat]

A possible bug fix is to set a = a - mu and b = b - mu before checking for the special cases, which currently happens only after at line 182/185.

Hi! Thank you! It's done on dev brunch.

So IMHO simplest way to add directly checking for special cases a=[0, -Inf], b=[Inf, 0]; a=[-Inf, -Inf], b=[0, 0]; a=[Inf, 0], b=[0, -Inf],

[CarloLucibello]

can this be closed?