Difference in output from mvtnorm and Mathematica/Java [closed]

Question

Let $a \sim\mathcal{N}(6.532056,0.06532056)$,$b \sim\mathcal{N}(8.390961,0.08390961)$ and $c \sim\mathcal{N}(8.736566,0.08736566)$. We use $\mathcal{N}(\mu,\sigma^2)$ notation unless specified otherwise.

We construct two normal variables $x = a-b$ and $y = a-c$. So, $x \sim\mathcal{N}(-1.858905,0.14923017)$ and $y \sim\mathcal{N}(-2.20451,0.15268622)$

Correlation between $x$ and $y$ i.e. $\textrm{cor}(x,y)= \frac{\textrm{var}(a)}{\sigma_x \times \sigma_y}$ by using the basic properties of covariance.

Solving, $\textrm{cor}(x,y)=0.4327346392418512 \approx 0.433$ which can be written in matrix (let it be called $ mat $)for as \begin{bmatrix} 1.0 & 0.433 \\[0.3em] 0.433 & 1.0 \\[0.3em] \end{bmatrix}

Now I want to find the $\Pr[-\infty<x<0 \textrm{ and } -\infty<y<0]$.

I used the mvtnorm package's pmvnorm method by invoking

pr<-pmvnorm(mean=c(-1.858905,-2.20451), corr=mat, lower=rep(-Inf, 2), upper=rep(0,2))

The result was 0.9575448.

The same when I want to compute in Mathematica/Java(an implementation of the Genz algorithm found on internet) I am getting the result as 0.9999992445813132.

I am including the Mathematica code below. Please note here, the second argument in this case, is standard deviation.

px = NormalDistribution[-1.858905000000001, 0.38630320992712447]

py = NormalDistribution[-2.20451, 0.3907508413298685]

pxy = ProductDistribution[px, py]

Probability[-Infinity < x < 0 && -Infinity < y < 0, {x, 
   y} \[Distributed] pxy]

What is the mistake I am doing? For my application domain (visualization) this difference is turning out to be very costly,especially because I run it over more than 700,000 data points.

Answer 1

Note that the covariance matrix of $(x, y)$ doesn't have unit diagonals, so specifying corr as an argument to pmvnorm provides incomplete information. You need to specify sigma, which is the covariance matrix.