1

I have a vector of means and a covariance matrix and I'm trying to create a data set that would like to generate strictly positive numbers that would fit the parameters. I have seen quite a few options but all of them resulted in creating some negative observations. I have implemented the following code using the MASS package, but I have also tried the mvtnorm as well as several other options which have not produced the needed results.

means <- c(0.1505, 0.1596, 0.1552)

cov_matrix <- matrix(c(0.23915, -0.04556, -0.04729,
                     -0.04556, 0.21460, -0.05269,
                     -0.04729, -0.05269, 0.20674), ncol = 3)


t <- MASS::mvrnorm(n = 50, mu = means, Sigma = cov_matrix, empirical = TRUE)


head(t)
        [,1]     [,2]     [,3]
[1,]  0.1432  0.51024  0.29832
[2,]  1.0483  0.29353 -1.02127
[3,]  0.3252  0.45140 -0.41820
[4,] -0.3459  0.08166  0.04431
[5,] -0.3144  0.29057  0.55096
[6,]  1.1736 -0.63583 -0.56427

I feel as if it should be quite simple and I'm just missing something.

Amidavid
  • 21
  • 3
  • 3
    If you want strictly positive numbers when why use the normal distribution? It has support for negative numbers whatever parameters you use. – Robert Long Jul 04 '20 at 11:50
  • 3
    It could be quite simple, but you need to tell us what joint distribution you want the samples to have. If you want them to be multivariate normal as per your code, then you would need the means to be large compared to their variance and maybe run the code with different seeds until you obtained strictly positive values. On the other hand you could choose a distribution that does not have support for negative numbers such as lognormal or gamma. – Robert Long Jul 04 '20 at 13:04

2 Answers2

1

Here is one suggested path.

Start by assuming that your data follows, for example, a Lognormal distribution, which provides a positive continuous probability function supported on the positive real numbers.

Required are the associated relationships in the Lognormal population between your target mean for, say variable Y, namely $Y_m$, and also its variance, $Var(Y)$, together with their relationship to the parameters of the generated lognormal deviate. The latter expressions for the mean is given by $\exp( \mu + \sigma^2/2)$ and the $Var(Y)$ is equal to $\exp(\sigma^2 -1)\exp(2\mu + \sigma^2)$. Also, per these two references, Wikipedia and this forum, we note:

$\text{Cov}(X,Y)=E(X)E(Y)[\exp(\rho\sigma_x\sigma_y)-1]$

So, upon dividing both sides of the expression by the product, $\sigma_X\sigma_Y$ (that is, the product of the respective standard deviations of X and Y in Lognormal population), we have a relationship connecting the parent bivariate Lognormal population target correlation with the $\rho$ associated with the correlation associated between the logs of X and Y, which I call simply ρ.

As the log of a Lognormal deviate is Normally distributed, one can now use relationships associated with the conditional distribution of $X$ given a value for $Y$, namely:

$E(X|Y=y) = \mu_X+ σ_X ρ(y−Y_m )/σ_Y $

$Var(X|Y=y) = σ_X^2 (1−ρ^2)$

Source, see, for example, Page 267 here.

So, the simulation procedure is to derive the corresponding parameter values of means, variances and correlation for the log of the lognonormal deviates. As these are Normally distributed, proceed to generate an appropriate random normal deviate for Y, and use the conditional Normal distribution expression to obtain the corresponding X. Repeat.

You now have a simulated population of bivariate Normal random deviates. Upon applying the exponential transform, you should now have a Lognormal population with the targetted parameter values.

AJKOER
  • 1,800
  • 1
  • 9
  • 9
1

Originally I preferred not to transform my distribution to log-normal, yet I have ended up doing so and I used the dmutate package for generating my numbers:

t <- rlmassnorm(n = 1000, mu = means, Sigma = cov_matrix, empirical = TRUE)
Amidavid
  • 21
  • 3