# Difference in output from R's mvtnorm and Mathematica/Java

Let a=$\mathcal{N}(6.532056,0.06532056)$,b~$\mathcal{N}(8.390961,0.08390961)$ and c~$\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=$\mathcal{N}(-1.858905,0.14923017)$ and y=$\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.

I believe you are after MultinormalDistribution:
mnd=MultinormalDistribution[{-1.858905000000001, -2.20451}, {{1,
0.4327346392418512}, {0.0 .4327346392418512, 1}}];
`