Skip to content

Instantly share code, notes, and snippets.

@arraytools
Created November 24, 2020 23:21
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
Star You must be signed in to star a gist
Save arraytools/0d7f0a02c233aefb9cefc6eb5f7b7754 to your computer and use it in GitHub Desktop.
R generate multivariate normal
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
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