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}}];
NProbability[-Infinity < x < 0 && -Infinity < y < 0, {x, 
   y} \[Distributed] mnd]

yields: 0.957547

| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.