R Under development (unstable) (2024-09-01 r87083 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> 
> library("mvtnorm")
> 
> chk <- function(...) isTRUE(all.equal(...))
> 
> m <- 1:3
> 
> s <- diag(1:3)
> s[2,1] <- 1
> s[3,1] <- 2
> s[3,2] <- 3
> s <- s+t(s)
> 
> set.seed(1)
> 
> x <- rmvnorm(10000, m, s)
> stopifnot(chk(m, colMeans(x), tolerance=0.01))
> stopifnot(chk(s, var(x), tolerance=0.1))
> 
> x <- rmvnorm(10000, m, s, method="svd")
> stopifnot(chk(m, colMeans(x), tolerance=0.01))
> stopifnot(chk(s, var(x), tolerance=0.1))
> 
> x <- rmvnorm(10000, m, s, method="chol")
> stopifnot(chk(m, colMeans(x), tolerance=0.01))
> stopifnot(chk(s, var(x), tolerance=0.1))
> 
> ### suggested by Paul Johnson <pauljohn@ku.edu>
> set.seed(29)
> x <- rmvnorm(2, sigma = diag(2))
> set.seed(29)
> y <- rmvnorm(3, sigma = diag(2))[1:2,]
> stopifnot(chk(x, y))
> 
> ### Speed
> p <- 200
> set.seed(17)
> rcond(Sig <- cov(matrix(rnorm((p+p)*p), ncol = p))) # 0.00286, ok
[1] 0.002862034
> mu <- 1:p
> 
> set.seed(101)
> system.time(x <- rmvnorm(10000, mu, Sig))
   user  system elapsed 
   0.62    0.00    0.63 
> stopifnot(chk(mu, colMeans(x), tolerance= 0.001),
+           chk(Sig, cov(x), tolerance= 0.2))
> 
> set.seed(101)
> system.time(x <- rmvnorm(10000, mu, Sig, method="svd"))
   user  system elapsed 
   0.58    0.05    0.63 
> stopifnot(chk(mu, colMeans(x), tolerance= 0.001),
+           chk(Sig, cov(x), tolerance= 0.2))
> 
> set.seed(101)
> system.time(x <- rmvnorm(10000, mu, Sig, method="chol"))
   user  system elapsed 
   0.48    0.05    0.53 
> ## 'chol' is  5-10 % faster than the other two
> stopifnot(chk(mu, colMeans(x), tolerance= 0.001),
+           chk(Sig, cov(x), tolerance= 0.2))
> 
> proc.time()
   user  system elapsed 
   3.43    0.20    3.62