CODEWITHSUNDEEP
rdistributiongaussiannormal-distribution
FairyOnIce
FairyOnIce

Reputation: 2614

R function pmvnorm: Why do values and errors differ every time I run this function with the same inputs?

In my understanding, pmvnorm in mvtnorm library is a function to compute the CDF over a multivariate normal distribution. So it is a deterministic function. However, I found that the results change every time I run this function with the same inputs. Here is a small example.

 library(mvtnorm)
 lower <- c( -Inf, -0.07,  0.81,  -Inf,  0.89,  -Inf,  1.33,  1.21,  -Inf)
 upper <- c( 1.00,  0.34,  0.98, -0.04,  1.07,  0.01,  1.48,  1.38,  0.09)

 sigma <- matrix(c(0.03, 0.00, -0.01, 0.00, 0.00, 0.00, 0.00, 0.00, 0.02,
                   0.00, 1.00,  0.66, 0.64, 0.64, 0.64, 0.64, 0.64, 0.64,
                  -0.01, 0.66,  0.99, 0.66, 0.64, 0.64, 0.64, 0.64, 0.64,
                   0.00, 0.64,  0.66, 1.00, 0.66, 0.64, 0.64, 0.64, 0.64,
                   0.00, 0.64,  0.64, 0.66, 1.00, 0.66, 0.64, 0.64, 0.64,
                   0.00, 0.64,  0.64, 0.64, 0.66, 1.00, 0.66, 0.64, 0.64,
                   0.00, 0.64,  0.64, 0.64, 0.64, 0.66, 1.00, 0.66, 0.64,
                   0.00, 0.64,  0.64, 0.64, 0.64, 0.64, 0.66, 1.00, 0.66,
                   0.02, 0.64,  0.64, 0.64, 0.64, 0.64, 0.64, 0.66, 1.00),
                   byrow=TRUE,length(lower))

 set.seed(1)
 (try1 <- pmvnorm(lower=lower,upper=upper,sigma=sigma))

This gives me the values:

 [1] 4.42436e-09
 attr(,"error")
 [1] 4.312159e-13
 attr(,"msg")
 [1] "Normal Completion"

Now I re-evaluate the function with a different seed.

 set.seed(2)
 (try2 <- pmvnorm(lower=lower,upper=upper,sigma=sigma))

Then I get:

 [1] 4.424396e-09
 attr(,"error")
 [1] 4.048187e-13
 attr(,"msg")
 [1] "Normal Completion"

And

 try1 == try2

gives me:

 [1] FALSE

Can anyone explain why this is happening?

Upvotes: 2

Views: 1588

Answers (1)

shadow
shadow

Reputation: 22333

Check out the first reference given in ?pmvnorm. For instance on http://www.math.wsu.edu/faculty/genz/papers/mvn.pdf. In short, pmvnorm usea a Monte Carlo sampling algorithm to calculate the distribution function of the multivariate normal distribution.

Upvotes: 4

Related Questions