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.
