I'm simulating a multivariate dataset using a modified form of the mvrnorm()
function (from the MASS package in R).
The problem is that I'm getting some negative values after the eigenvalue transformation (I assume) because I have many negative correlations and a number of small means.
Is there a special way to deal with this phenomenon? My idea is to just add some factor to all the datapoints (such as the 1st quartile) because I don't care so much what the exact means are, as long as the correlation structure remains intact.
A simple example:
require(MASS) # For mvrnorm()
require(Matrix) # For nearPD()
corr <- diag(5)
corr[5,1] <- .5
corr[1,5] <- .5
corr[4,2] <- -.5
corr[2,4] <- -.5
set.seed(1000)
mm <- mvrnorm(n=10, mu=rep(1,5), Sigma=nearPD(corr, corr=TRUE)$mat,
empirical=TRUE)
As you can see, mm has nonpositive examples. Since I want to model physical measurements, this makes no sense.
edit2: multivariate sampling from log-normal distribution still results in negative examples:
mvrlnorm <- function (n = 1, mu, Sigma, tol = 1e-06, empirical = TRUE) {
require(Matrix)
p <- length(mu)
if (!all(dim(Sigma) == c(p, p)))
stop("incompatible arguments")
eS <- eigen(Sigma, symmetric = TRUE, EISPACK = TRUE)
ev <- eS$values
if (!all(ev >= -tol * abs(ev[1L])))
stop("'Sigma' is not positive definite")
# HERE BE log-normal distribution
X <- matrix(rlnorm(p * n), nrow=n)
if (empirical) {
X <- scale(X, TRUE , FALSE)
X <- X %*% svd(X, nu = nrow(X), nv = ncol(X))$v
X <- scale(X, FALSE, TRUE)
}
retval <- eS$vectors %*% diag(sqrt(ev), length(ev)) %*% t(eS$vectors)
retval <- X %*% retval
retval <- sweep(retval, 2, mu, "+")
X <- retval
nm <- names(mu)
if (is.null(nm) && !is.null(dn <- dimnames(Sigma)))
nm <- dn[[1L]]
dimnames(X) <- list(nm, NULL)
if (n == 1)
drop(X)
else t(X)
}
rep(1,5)
instead of, say, $10$, which would make your problem go away? $\endgroup$ – whuber♦ Mar 28 '12 at 16:47