# R Fortran using rejection sampling in rtmvnorm() gives error

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. • Ideally, give a reproducible example. At the very least, please explain the inputs that produced the message, and quote the exact message. May 20, 2015 at 13:02 ## 1 Answer 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"