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
)
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 |
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 |
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 |
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.
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;
Matthew Fidler, Zdravko Botev and some from Matteo Fasiolo
## 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