[R] Question regarding the mvtnorm R package

[shlomi lifshits]

Hi,

I am trying to compare two methods for the calculation of a trivariate normal cdf:

1) "direct" - using the original covariance matrix
2) "indirect" -  transforming the integration interval using the Mahalanobis transformation and then integrating as a standard random variable.

The following code demonstrates the two methods.

Unfortunately, I get different results.

Am I missing something?

Thanks in advance.

Shlomi Lifshits
Tel-Aviv University

####################

library(mvtnorm)

#create a symmetric positive definite matrix for the covariance matrix
my_mat<-matrix(data=c(1,2,3,0,0.1,0,0,0,0.1), nrow=3, ncol=3, byrow=TRUE)
covmat<-my_mat%*%t(my_mat)


#calculate the Mahalanobis transformation 
spectral_dec<-eigen(covmat)
new_eigenval<-c(sqrt(spectral_dec$values[1]),sqrt(spectral_dec$values[2]),sqrt(spectral_dec$values[3]))
covmat_trans<-solve(spectral_dec$vectors%*%diag(new_eigenval)%*%t(spectral_dec$vectors))

# integrate with the original limit
orig_limit<-c(-0.5,0.5,0.5)
a<-pmvnorm(lower=-Inf, upper=orig_limit, mean=rep(0, 3), sigma=covmat)

#transform the limit and integrate
new_limit_raw<-covmat_trans%*%orig_limit
new_limit<-c(new_limit_raw[1,1],new_limit_raw[2,1],new_limit_raw[3,1])
new_sigma<-diag(3)
b<-pmvnorm(lower=-Inf, upper=new_limit, mean=rep(0, 3), sigma=new_sigma)

# a,b should be the same but they are not!!!?