[R] mnormt and mvtnorm

[Shh]

Hi R user,

I came across a question when I played around with mnormt package (mvtnorm as well), could any of you let me know where I made the mistake?  Thanks!

For instance, with a multivariate normal distribution with known mean and variance (to be simple, I tried the i.i.d case first)

 u1 <- u2 <- u3 <- 1
 test <- function(x) {2*x*dnorm(x,u1)* pmnorm(rep(x,1),c(u2),diag(1)) }
 value1 <- integrate(test , lower=-Inf, upper=Inf)$value
 testA <- function(x) {2*x*dnorm(x,u1,1)*pnorm(x,u2)}
 value1A <- integrate(testA , lower=-Inf, upper=Inf)$value
 value1
#[1] 1.564190
 value1A
#[1] 1.564190

# when pmnorm handle the univariate case, it is the same as pnorm. However, when I extended to multivariate, I got a confusing result:
 
 test2 <- function(x) {x*dnorm(x,u1)*pmnorm(rep(x,2),c(u2,u3),diag(2)) }
 value2 <- integrate(test2 , lower=-Inf, upper=Inf)$value
 test2A <- function(x) {x*dnorm(x,u1)*pnorm(x,u2)*pnorm(x,u3)}
 value2A <- integrate(test2A , lower=-Inf, upper=Inf)$value
 value2
# [1] 0.6997612
 value2A
# [1] 0.6154281

# I assume value2=value2a since they are independent, am I right? Also when I used some acutual values, I got:
 
 pmnorm(rep(0,2),c(u2,u3),diag(2)) 
# [1] 0.02517149
 pnorm(0,u2)*pnorm(0,u3)
# [1] 0.02517149

 test3  <- function(x) {x*dnorm(x,u1)*dmnorm(rep(x,2),c(u2,u3),diag(2)) }
 value3 <- integrate(test3 , lower=-Inf, upper=Inf)$value
 test3A <- function(x) {x*dnorm(x,u1)*dnorm(x,u2)*dnorm(x,u3)}
 value3A <- integrate(test3A , lower=-Inf, upper=Inf)$value
 value3
#[1] 0
 value3A
# [1]  0.09188815

# I got they same with an actual value:
 
 dmnorm(rep(0,2),c(u2,u3),diag(2)) 
#[1] 0.05854983
 dnorm(0,u2)*dnorm(0,u3)
# [1] 0.05854983


Similar thing occurred in mvtnorm I got three different results with mvtnorm, mnormt and pnorm.