I am trying to generate random vectors with the same covariance, mean and sd as a set I have. I am using R's mvrnorm with Empirical = F because the underlying distribution is in general not normal. However when running this many times I get the impression that specific elements of the returned vectors are too similar to the original vector. Is there a way to verify this quantitatively? and if this is indeed the case, how can I generate vectors that are "completely different" from the original (while maintaining same sd and mean)?

For example:

> df
     A        B
1 1.138289 1.918876
2 1.621082 2.432379
3 1.295929 2.241403
4 1.325516 1.845258

mvrnorm(n = 4, mu = center, Sigma = as.matrix(cov(df)), empirical = F)
        A        B
[1,] 1.364986 2.235192
[2,] 1.116274 1.828175
[3,] 1.438028 2.256231
[4,] 1.586585 1.956803

Although the bias cannot be shown from this example (but only for a large number of samples), there are pairs which are the same as the original (up to 2 digits after the decimal point), which are returned significantly more than others.

  • 1
    $\begingroup$ While this is not a forum for debugging code, could you perhaps share the R-code to make it clear what you are doing, and what results you are seeing? $\endgroup$ – Repmat Feb 20 '18 at 7:13
  • 2
    $\begingroup$ (1) The empirical argument has nothing to do with normality of the original sample, the function will return normal draws in each case. (2) Why exactly you consider the draws to be strange..? How large are the samples you are taking about? If you'd toss two coins simultaneously, then at some point you'll toss two heads at random... $\endgroup$ – Tim Feb 20 '18 at 8:59

I'm assuming in my answer that center = colMeans(df).

You don't have to worry about mvrnorm "replicating" your data frame. The function doesn't even "know" what your data looks like (you only pass the mean and covariance matrix). The randomization in mvrnorm depends on rnorm which is very well tested. If the patterns you see don't seem random enough for you it is probably because the human eye is very bad at discerning randomness.

mvrnorm will always produce samples from a normal. The option empirical will scale the sample so that the sample mean and covariance equal exactly the numbers that you provided. It has nothing to do with accommodating non-normal distributions. To see this, run the following code:

df <- data.frame(A = c(1.138289, 1.621082, 1.295929, 1.325516), 
B = c(1.918876, 2.432379, 2.241403, 1.845258))
my_sample <- MASS::mvrnorm(n = 100, mu = colMeans(df), Sigma = as.matrix(cov(df)), 
  empirical = TRUE)

If you want to maintain mean and covariance structure then you want empirical = TRUE.

| cite | improve this answer | |
  • $\begingroup$ I would like to maintain mean and covariance per sample, which holds for both empirical = True or False. My suspicion that very close values are returned are corroborated by running the code 1000 times and seeing that a specific pair is returned more than others. But how can I check this suspicion? Is there a measure that will indicate "perfect randomness" for many samples? Something like comparing within cluster distance for the original vs. the sampled vector? $\endgroup$ – Noale Feb 20 '18 at 8:51
  • $\begingroup$ I don't think the randomization is the problem here. You are aware that you are generating highly correlated data? $\endgroup$ – Andreas Dzemski Feb 20 '18 at 9:06
  • $\begingroup$ Yes. I am starting to think that since the algorithm assumes underlying normality, the pair which is closest to the mean will be sampled more than a pair close to the tail. Could this be the issue? $\endgroup$ – Noale Feb 20 '18 at 9:10
  • $\begingroup$ That's right (where "closeness" will not be measured by the usual Euclidean metric because of the correlation). Maybe you are looking for a way to generate samples from your distribution nonparametrically? How to do this is probably the topic of a new question. The easiest way is resampling from the empirical distribution. $\endgroup$ – Andreas Dzemski Feb 20 '18 at 12:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.