Why is the “pmvnorm” result different when the input matrix are covariance and correlation matrix?

Question

The result is different when using a variance matrix and a correlation matrix. Why is this happening?

I will write down the results directly for convenience.

Variance matrix - naming as co
0.1234 0.125
0.1250 0.245

Correlation matrix - naming as coo (made by cov2cor function)
1.0000 0.7189
0.7189 1.0000

Result
pmvnorm(mean=c(1,1),sigma=co, lower=rep(-Inf,2), upper=c(0.7,4)
0.1965493

pmvnorm(mean=c(1,1),corr=coo, lower=rep(-Inf,2), upper=c(0.7,4)
0.3820885

I made a covariance matrix, and we got a correlation matrix using covariance matrix. And these two values were implemented, and the result was different.

It is code.

install.packages("mvtnorm")
library(mvtnorm)
co <- matrix(c(0.1234,0.125,0.125,0.245),2,2)
coo <- cov2cor(co)
pmvnorm(mean=c(1,1),sigma=co, lower=rep(-Inf,2), upper=c(0.7,4)
pmvnorm(mean=c(1,1),corr=coo, lower=rep(-Inf,2), upper=c(0.7,4)

Please let me know why.

Maurits Evers Maurits Evers · Accepted Answer · 2017-11-21T01:35:41

As per ?pmvnorm (emphasis mine)

sigma: the covariance matrix of dimension n. Either ‘corr’ or ‘sigma’ can be specified. If ‘sigma’ is given, the problem is standardized. If neither ‘corr’ nor ‘sigma’ is given, the identity matrix is used for ‘sigma’.

So to make both calculations consistent you need to give standardised upper limits when giving a correlation matrix.

# Using covariance matrix sigma
cov <- matrix(c(0.1234,0.125,0.125,0.245), 2, 2);
x1 <- pmvnorm(mean = c(1, 1), sigma = cov, lower = -Inf, upper = c(0.7, 4));
x1;
#[1] 0.1965493
#attr(,"error")
#[1] 1e-15
#attr(,"msg")
#[1] "Normal Completion"

# Using correlation matrix corr
# Note: Need to scale the upper limits
cor <- cov2cor(cov);
x2 <- pmvnorm(mean = c(0, 0), corr = cor, lower = -Inf, upper = (c(0.7, 4) - c(1, 1)) / sqrt(diag(cov)));
x2;
#[1] 0.1965493
#attr(,"error")
#[1] 1e-15
#attr(,"msg")
#[1] "Normal Completion"

PS. It's a bit hidden, but ?pmvnorm includes a simpler example on the complementarity of both approaches.

 # Correlation and Covariance

 a <- pmvnorm(lower=-Inf, upper=c(2,2), sigma = diag(2)*2)
 b <- pmvnorm(lower=-Inf, upper=c(2,2)/sqrt(2), corr=diag(2))
 stopifnot(all.equal(round(a,5) , round(b, 5)))

Why is the “pmvnorm” result different when the input matrix are covariance and correlation matrix?

1 Answers