rxRmvn function

Simulate from a (truncated) multivariate normal

This is simulated with the fast, thread-safe threefry simulator and can use multiple cores to generate the random deviates.

rxRmvn(
  n,
  mu = NULL,
  sigma,
  lower = -Inf,
  upper = Inf,
  ncores = 1,
  isChol = FALSE,
  keepNames = TRUE,
  a = 0.4,
  tol = 2.05,
  nlTol = 1e-10,
  nlMaxiter = 100L
)

Arguments

• n: Number of random row vectors to be simulated OR the matrix to use for simulation (faster).

• mu: mean vector

• sigma: Covariance matrix for multivariate normal or a list of covariance matrices. If a list of covariance matrix, each matrix will simulate n matrices and combine them to a full matrix

• lower: is a vector of the lower bound for the truncated multivariate norm

• upper: is a vector of the upper bound for the truncated multivariate norm

• ncores: Number of cores used in the simulation

• isChol: A boolean indicating if sigma is a cholesky decomposition of the covariance matrix.

• keepNames: Keep the names from either the mean or covariance matrix.

• a: threshold for switching between methods; They can be tuned for maximum speed; There are three cases that are considered:

  case 1: a < l < u

  case 2: l < u < -a

  case 3: otherwise

  where l=lower and u = upper

• tol: When case 3 is used from the above possibilities, the tol value controls the acceptance rejection and inverse-transformation;

  When abs(u-l)>tol, uses accept-reject from randn

• nlTol: Tolerance for newton line-search

• nlMaxiter: Maximum iterations for newton line-search

Returns

If n==integer (default) the output is an (n x d) matrix where the i-th row is the i-th simulated vector.

If is.matrix(n) then the random vector are store in n, which is provided by the user, and the function returns NULL invisibly.

Examples

## From mvnfast
## Unlike mvnfast, uses threefry simulation

d <- 5
mu <- 1:d

# Creating covariance matrix
tmp <- matrix(rnorm(d^2), d, d)
mcov <- tcrossprod(tmp, tmp)

set.seed(414)
rxRmvn(4, 1:d, mcov)

set.seed(414)
rxRmvn(4, 1:d, mcov)

set.seed(414)
rxRmvn(4, 1:d, mcov, ncores = 2) # r.v. generated on the second core are different

###### Here we create the matrix that will hold the simulated
#  random variables upfront.
A <- matrix(NA, 4, d)
class(A) <- "numeric" # This is important. We need the elements of A to be of class "numeric".

set.seed(414)
rxRmvn(A, 1:d, mcov, ncores = 2) # This returns NULL ...
A # ... but the result is here

## You can also simulate from a truncated normal:

rxRmvn(10, 1:d, mcov, lower = 1:d - 1, upper = 1:d + 1)

# You can also simulate from different matrices (if they match
# dimensions) by using a list of matrices.

matL <- lapply(1:4, function(...) {
  tmp <- matrix(rnorm(d^2), d, d)
  tcrossprod(tmp, tmp)
})

rxRmvn(4, setNames(1:d, paste0("a", 1:d)), matL)

References

John K. Salmon, Mark A. Moraes, Ron O. Dror, and David E. Shaw (2011). Parallel Random Numbers: As Easy as 1, 2, 3. D. E. Shaw Research, New York, NY 10036, USA.

The thread safe multivariate normal was inspired from the mvnfast package by Matteo Fasiolo https://CRAN.R-project.org/package=mvnfast

The concept of the truncated multivariate normal was taken from Zdravko Botev Botev (2017) tools:::Rd_expr_doi("10.1111/rssb.12162")

and Botev and L'Ecuyer (2015) tools:::Rd_expr_doi("10.1109/WSC.2015.7408180") and converted to thread safe simulation;

Author(s)

Matthew Fidler, Zdravko Botev and some from Matteo Fasiolo