CODEWITHSUNDEEP
r
AAA
AAA

Reputation: 27

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

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.

Upvotes: 1

Views: 535

Answers (1)

Maurits Evers
Maurits Evers

Reputation: 50718

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)))

Upvotes: 1

Related Questions