[R] Question regarding the mvtnorm R package
shlomi lifshits
lifshits_il at yahoo.com
Tue Apr 10 11:49:21 CEST 2007
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!!!?
More information about the R-help
mailing list