4
$\begingroup$

I'm sampling from a multivariate normal truncated distribution using rtmvnorm() function from tmvtnorm package in R. Using the rejection algorithm method, I can check the acceptance rate using pmvnorm(), as described in the article by Wilhelm & Manjunath (2010) who explain the package.

Here is an example of my code:

# creating covariance matrix using different random correlations 
n = 8    
c <- c(rnorm((n*(n-1))/2, mean = .1, sd = 1))
sigma <- diag(8)
sigma[lower.tri(sigma, diag = F)] <- c
sigma <- forceSymmetric(sigma, "L") # creating symmetric correlation matrix
sigma <- sigma * (sqrt(1.5) * sqrt(1.5)) #creating covariance matrix from correlation matrix
sigma <- as.matrix(nearPD(sigma)$mat) #make covariance matrix positive definite

#random sampling multivariate truncated distribution using sigma
mean <- rnorm(n, 6.5, 1.5)    
t <- rtmvnorm(10, mean = mean, sigma = sigma, lower = c(rep(1,n)), upper = c(rep(10,n)),
              algorithm = "rejection")
alpha <- pmvnorm(lower = c(rep(1,n)), upper = c(rep(10,n)), mean = mean, sigma = sigma)
[1] 0.6900399
attr(,"error")
[1] 0.001841024
attr(,"msg")
[1] "Completion with error > abseps"

Now, my acceptance rate .69 is alright, and increases if I simulate a covariance matrix that has an higher mean (instead of .1). However, the accompanied error message indicates that the completion is with error > abseps. I found that this means my error is higher than the absolute error tolerance of .001. I was wondering if someone knows what this error tolerance indicates exactly?

I've tried looking up the codes but it gets me to mtv() from mvtnorm package in R which refers to a function written in Fortran, MVTDST, http://www.math.wsu.edu/faculty/genz/software/fort77/mvtdstpack.f. As I am unfamiliar with this language. I'm not quite sure what produces the message and how bad it would be to lower the error or to ignore the message.

$\endgroup$
1
  • 3
    $\begingroup$ Ideally, give a reproducible example. At the very least, please explain the inputs that produced the message, and quote the exact message. $\endgroup$
    – Glen_b
    May 20, 2015 at 13:02

1 Answer 1

1
$\begingroup$

For what it is worth, one can get a Normal Completion message in some cases by increasing the maxpts -- which can be read about in ?mvtnorm::GenzBretz. The default for maxpts is 25000. Below is the OP code with library() and set.seed() statements to make it reproducible and then a demonstration of increasing the maxpts from 25000 to 10*25000 in the alpha2 call. The expected answers are commented out below the pmvnorm() calls:

library(tmvtnorm)
library( mvtnorm)

set.seed(424242)

# creating covariance matrix using different random correlations 
n = 8    
c <- c(rnorm((n*(n-1))/2, mean = .1, sd = 1))
sigma <- diag(8)
sigma[lower.tri(sigma, diag = F)] <- c
sigma <- forceSymmetric(sigma, "L") # creating symmetric correlation matrix
sigma <- sigma * (sqrt(1.5) * sqrt(1.5)) #creating covariance matrix from correlation matrix
sigma <- as.matrix(nearPD(sigma)$mat) #make covariance matrix positive definite

#random sampling multivariate truncated distribution using sigma
mean <- rnorm(n, 6.5, 1.5)    
t <- rtmvnorm(10, mean = mean, sigma = sigma, lower = c(rep(1,n)), upper = c(rep(10,n)),
              algorithm = "rejection")
alpha <- pmvnorm(lower = c(rep(1,n)), upper = c(rep(10,n)), mean = mean, sigma = sigma)
alpha
# [1] 0.7572689
# attr(,"error")
# [1] 0.002857951
# attr(,"msg")
# [1] "Completion with error > abseps"


## increase maxpts
alpha2 <- pmvnorm(lower = c(rep(1,n)), upper = c(rep(10,n)), mean = mean, sigma = sigma, 
                  maxpts=10*25000)
alpha2
# [1] 0.7569266
# attr(,"error")
# [1] 0.0007923042
# attr(,"msg")
# [1] "Normal Completion"
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.