Generating correlated positive random numbers (given means, variances and degree of correlation)? 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.
 A: 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.
A: 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)

