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.