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

Value

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.

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) doi: 10.1111/rssb.12162 and Botev and L'Ecuyer (2015) doi: 10.1109/WSC.2015.7408180 and converted to thread safe simulation;

Author

Matthew Fidler, Zdravko Botev and some from Matteo Fasiolo

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)
#>           [,1]        [,2]      [,3]     [,4]     [,5]
#> [1,] 0.6136552  1.25812013  5.303832 4.525316 5.365598
#> [2,] 0.2516852 -0.92219464  1.078185 1.306427 2.255198
#> [3,] 3.1149151  0.09902607 -1.774739 4.105226 6.016358
#> [4,] 0.2015307  2.66830901  2.041887 2.379933 3.566044

set.seed(414)
rxRmvn(4, 1:d, mcov)
#>           [,1]        [,2]      [,3]     [,4]     [,5]
#> [1,] 0.6136552  1.25812013  5.303832 4.525316 5.365598
#> [2,] 0.2516852 -0.92219464  1.078185 1.306427 2.255198
#> [3,] 3.1149151  0.09902607 -1.774739 4.105226 6.016358
#> [4,] 0.2015307  2.66830901  2.041887 2.379933 3.566044

set.seed(414)
rxRmvn(4, 1:d, mcov, ncores = 2) # r.v. generated on the second core are different
#>           [,1]      [,2]       [,3]     [,4]     [,5]
#> [1,] 0.6136552 4.5619767  3.7840009 3.820684 6.181438
#> [2,] 0.3652661 2.3724185  5.0871725 4.070487 5.749454
#> [3,] 0.2516852 0.5206738 -0.4249075 1.026492 1.182587
#> [4,] 1.3164949 2.0964747  1.4628823 3.371607 5.266312

###### 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
#>           [,1]      [,2]       [,3]     [,4]     [,5]
#> [1,] 0.6136552 4.5619767  3.7840009 3.820684 6.181438
#> [2,] 0.3652661 2.3724185  5.0871725 4.070487 5.749454
#> [3,] 0.2516852 0.5206738 -0.4249075 1.026492 1.182587
#> [4,] 1.3164949 2.0964747  1.4628823 3.371607 5.266312

## You can also simulate from a truncated normal:

rxRmvn(10, 1:d, mcov, lower = 1:d - 1, upper = 1:d + 1)
#>             [,1]     [,2]     [,3]     [,4]     [,5]
#>  [1,] 1.40194307 2.188182 2.829579 4.998496 4.804913
#>  [2,] 0.02971623 2.555220 2.694816 4.462161 5.022664
#>  [3,] 0.73953497 2.082153 2.918758 4.984464 5.057942
#>  [4,] 0.18577166 2.453988 2.695229 4.462135 4.990005
#>  [5,] 1.00918396 2.113297 3.013401 3.182227 5.226960
#>  [6,] 0.02691318 1.336583 3.166658 4.139035 5.284742
#>  [7,] 0.75784792 1.776195 3.369402 4.195169 4.996142
#>  [8,] 1.12122571 1.662802 2.679519 4.951974 5.025238
#>  [9,] 0.65771584 1.992188 2.946981 4.750599 5.124848
#> [10,] 0.08909610 1.817789 2.978310 3.115197 5.301765


# 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)
#>               a1         a2         a3       a4       a5
#>  [1,]  2.1027689 -3.1277802 -1.7058664 3.877449 6.160238
#>  [2,] -0.1528749  0.3372425  1.4079033 3.204818 6.771887
#>  [3,] -1.6178920  2.4185331  4.2854089 4.204923 6.115151
#>  [4,]  5.2850177  6.0162640  1.1522397 6.009218 8.423744
#>  [5,] -0.8389005  1.8108360  3.9390837 2.635257 3.078044
#>  [6,]  2.4044119  0.2584573  2.4396218 4.687467 7.747550
#>  [7,]  7.2635180  5.5533803 -0.2744354 9.627871 4.788654
#>  [8,] -0.3358839 -1.7501935  4.1309667 2.961912 6.318101
#>  [9,]  3.4028706  1.1584222  4.5627803 3.408984 3.627014
#> [10,]  3.4927361 -1.5055971  1.0890217 2.867575 4.036955
#> [11,]  2.2089477 -1.0108105 -0.6234824 1.745858 4.178149
#> [12,]  1.9045759 -3.5353019  3.6738287 3.085417 2.755276
#> [13,] -1.0910176  1.6896066 -0.1344867 5.810620 3.303876
#> [14,]  0.7032309  2.5654064  2.7642640 3.706315 5.057585
#> [15,] -3.4067706  2.3051387  4.4260234 6.692458 8.170364
#> [16,] -3.7727314  4.2737997 -0.8379328 6.990844 1.255809