CODEWITHSUNDEEP
r
Matata
Matata

Reputation: 271

create multi-normal data in R using mvrnorm

I was using mvrnorm to generate data from multi-normal distribution with the mean mu <- rep(0,4) and Sigma be some positive definite symmetric matrix. However, I found that the last element in the vector I am generating is always 0, any ideas about why is that?

> mvrnorm(n = 1, mu, Sigma, tol = 1e-6, empirical = FALSE, EISPACK = FALSE)
[1]  0.1813268 -0.8993918  0.7461007  0.0000000
> mvrnorm(n = 1, mu, Sigma, tol = 1e-6, empirical = FALSE, EISPACK = FALSE)
[1] 3.2539025 2.9855514 0.7313427 0.0000000
> mvrnorm(n = 1, mu, Sigma, tol = 1e-6, empirical = FALSE, EISPACK = FALSE)
[1] -0.8133201 -1.0011971 -0.3800518  0.0000000

Thanks in advance!

Edit: Thanks for the responses, yes, I checked the Sigma, there is something wrong with it.

Upvotes: 3

Views: 565

Answers (1)

Zheyuan Li
Zheyuan Li

Reputation: 73265

The reason is that your Sigma is not deemed full rank under tol = 1e-6. However, the way that mvrnorm does rank-detection is a bit odd. See inside MASS::mvrnorm:

eS <- eigen(Sigma, symmetric = TRUE)
ev <- eS$values
if (!all(ev >= -tol * abs(ev[1L]))) 
    stop("'Sigma' is not positive definite")
X <- matrix(rnorm(p * n), n)
#[...omitted...]
X <- drop(mu) + eS$vectors %*% diag(sqrt(pmax(ev, 0)), p) %*% t(X)

Instead of

ev >= tol * abs(ev[1L])

it does

ev >= -tol * abs(ev[1L])

Therefore, you must have negative eigenvalues to get rank-deficiency.

Upvotes: 3

Related Questions