Created
November 24, 2020 23:21
Star
You must be signed in to star a gist
R generate multivariate normal
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
set.seed(1234) | |
junk <- biospear::simdata(n=500, p=500, q.main = 10, q.inter = 10, | |
prob.tt = .5, m0=1, alpha.tt= -.5, | |
beta.main= -.5, beta.inter= -.5, b.corr = .7, b.corr.by=25, | |
wei.shape = 1, recr=3, fu=2, timefactor=1) | |
## Method 1: MASS::mvrnorm() | |
## This is simdata() has used. It gives different numbers on different OS. | |
## | |
library(MASS) | |
set.seed(1234) | |
n <- 500 | |
p <- 500 | |
b.corr.by <- 25 | |
b.corr <- .7 | |
n.blocks <- p%/%b.corr.by | |
covMat <- diag(n.blocks) %x% | |
matrix(b.corr^abs(matrix(1:b.corr.by, b.corr.by, b.corr.by, byrow = TRUE) - | |
matrix(1:b.corr.by, b.corr.by, b.corr.by)), b.corr.by, b.corr.by) | |
diag(covMat) <- 1 | |
x <- MASS::mvrnorm(n, rep(0, p), Sigma = covMat) | |
range(x) | |
eS <- eigen(covMat, symmetric = TRUE) # used by MASS::mvrnorm() | |
eS$vectors[which(eS$vectors[, 500] != 0), 500] # mac and Linux have opposite sign | |
packageVersion("MASS") | |
# Mac: [1] ‘7.3.49’ | |
# Linux: [1] ‘7.3.49’ | |
# Windows: [1] ‘7.3.47’ | |
R.version$version.string | |
# Mac: [1] "R version 3.4.3 (2017-11-30)" | |
# Linux: [1] "R version 3.4.4 (2018-03-15)" | |
# Windows: [1] "R version 3.4.3 (2017-11-30)" | |
## Method 2: mvtnorm::rmvnorm() | |
set.seed(1234) | |
x <- mvtnorm::rmvnorm(n=n, rep(0, p), sigma=covMat) | |
# NOTE: eigen is used by default | |
range(x) | |
# Mac: [1] -4.482566 4.459236 | |
# Linux: [1] -4.482566 4.459236 | |
## Method 3: mvnfast::rmvn() | |
set.seed(1234) | |
x <- mvnfast::rmvn(n, rep(0, p), covMat) | |
range(x) | |
# Mac: [1] -4.323585 4.355666 | |
# Linux: [1] -4.323585 4.355666 | |
library(microbenchmark) | |
library(MASS) | |
library(mvtnorm) | |
library(mvnfast) | |
microbenchmark(v1 <- rmvnorm(n=n, rep(0, p), sigma=covMat, "eigen"), | |
v2 <- rmvnorm(n=n, rep(0, p), sigma=covMat, "svd"), | |
v3 <- rmvnorm(n=n, rep(0, p), sigma=covMat, "chol"), | |
v4 <- rmvn(n, rep(0, p), covMat), | |
v5 <- mvrnorm(n, rep(0, p), Sigma = covMat)) | |
Unit: milliseconds | |
expr min lq | |
v1 <- rmvnorm(n = n, rep(0, p), sigma = covMat, "eigen") 296.55374 300.81089 | |
v2 <- rmvnorm(n = n, rep(0, p), sigma = covMat, "svd") 461.81867 466.98806 | |
v3 <- rmvnorm(n = n, rep(0, p), sigma = covMat, "chol") 118.33759 120.01829 | |
v4 <- rmvn(n, rep(0, p), covMat) 66.64675 69.89383 | |
v5 <- mvrnorm(n, rep(0, p), Sigma = covMat) 291.19826 294.88038 | |
mean median uq max neval cld | |
306.72485 301.99339 304.46662 335.6137 100 d | |
478.58536 470.44085 493.89041 571.7990 100 e | |
125.85427 121.26185 122.21361 151.1658 100 b | |
71.67996 70.52985 70.92923 100.2622 100 a | |
301.88144 296.76028 299.50839 346.7049 100 c |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
set.seed(1234); x <- mvrnorm(n, rep(0, p), Sigma = covMat) | |
debug(mvrnorm) | |
# eS --- macOS | |
# eS2 -- Linux | |
Browse[2]> range(abs(eS$values - eS2$values)) | |
# [1] 0.000000e+00 1.776357e-15 | |
Browse[2]> var(as.vector(eS$vectors)) | |
[1] 0.002000006 | |
Browse[2]> var(as.vector(eS2$vectors)) | |
[1] 0.001999987 | |
Browse[2]> all.equal(eS$values, eS2$values) | |
[1] TRUE | |
Browse[2]> which(eS$values != eS2$values) | |
[1] 6 7 8 9 10 11 12 13 14 20 22 23 24 25 26 27 28 29 | |
... | |
[451] 494 495 496 497 499 500 | |
Browse[2]> range(abs(eS$vectors - eS2$vectors)) | |
[1] 0.0000000 0.5636919 |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment