How to generate correlated random numbers (given means, variances and degree of correlation)?

Question

I'm sorry if this seems a bit too basic, but I guess I'm just looking to confirm understanding here. I get the sense I'd have to do this in two steps, and I've started trying to grok correlation matrices, but it's just starting to seem really involved. I'm looking for a concise explanation (ideally with hints towards a pseudocode solution) of a good, ideally quick way to generate correlated random numbers.

Given two pseudorandom variables height and weight with known means and variances, and a given correlation, I think I'm basically trying to understand what this second step should look like:

   height = gaussianPdf(height.mean, height.variance)
   weight = gaussianPdf(correlated_mean(height.mean, correlation_coefficient), 
                        correlated_variance(height.variance, 
                        correlation_coefficient))

• How do I calculate the correlated mean and variance? But I want to confirm that's really the relevant problem here.
• Do I need to resort to matrix manipulation? Or do I have something else very wrong in my basic approach to this problem?

Answer 1

To answer your question on "a good, ideally quick way to generate correlated random numbers": Given a desired variance-covariance matrix $C$ that is by definition positive definite, the Cholesky decomposition of it is: $C$=$LL^T$; $L$ being lower triangular matrix.

If you now use this matrix $L$ to project an uncorrelated random variable vector $X$, the resulting projection $Y = LX$ will be that of correlated random variables.

You can find an concise explanation why this happens here.

Answer 2

+1 to @user11852, and @jem77bfp, these are good answers. Let me approach this from a different perspective, not because I think it's necessarily better in practice, but because I think it's instructive. Here are a few relevant facts that we know already:

1. $r$ is the slope of the regression line when both $X$ and $Y$ are standardized, i.e., $\mathcal N(0,1)$,

2. $r^2$ is the proportion of the variance in $Y$ attributable to the variance in $X$,

   
   
   (also, from the rules for variances):

3. the variance of a random variable multiplied by a constant is the constant squared times the original variance:
   $$\text{Var}[aX]=a^2\text{Var}[X]$$
4. variances add, i.e., the variance of the sum of two random variables (assuming they are independent) is the sum of the two variances:
   $$\text{Var}[X+\varepsilon]=\text{Var}[X]+\text{Var}[\varepsilon]$$

Now, we can combine these four facts to create two standard normal variables whose populations will have a given correlation, $r$ (more properly, $\rho$), although the samples you generate will have sample correlations that vary. The idea is to create a pseudorandom variable, $X$, that is standard normal, $\mathcal N(0,1)$, and then find a coefficient, $a$, and an error variance, $v_e$, such that $Y \sim\mathcal N(0,a^2+v_e)$, where $a^2+v_e=1$. (Note that $|a|$ must be $\le 1$ for this to work, and that, moreover, $a=r$.) Thus, you start with the $r$ that you want; that's your coefficient, $a$. Then you figure out the error variance that you will need, it's $1-r^2$. (If your software requires you to use the standard deviation, take the square root of that value.) Finally, for each pseudorandom variate, $x_i$, that you have generated, generate a pseudorandom error variate, $e_i$, with the appropriate error variance $v_e$, and compute the correlated pseudorandom variate, $y_i$, by multiplying and adding.

If you wanted to do this in R, the following code might work for you:

correlatedValue = function(x, r){
  r2 = r**2
  ve = 1-r2
  SD = sqrt(ve)
  e  = rnorm(length(x), mean=0, sd=SD)
  y  = r*x + e
  return(y)
}

set.seed(5)
x = rnorm(10000)
y = correlatedValue(x=x, r=.5)

cor(x,y)
[1] 0.4945964

(Edit: I forgot to mention:) As I've described it, this procedure gives you two standard normal correlated variables. If you don't want standard normals, but want the variables to have some specific means (not 0) and SDs (not 1), you can transform them without affecting the correlation. Thus, you would subtract the observed mean to ensure that the mean is exactly $0$, multiply the variable by the SD you want and then add the mean you want. If you want the observed mean to fluctuate normally around the desired mean, you would add the initial difference back. Essentially, this is a z-score transformation in reverse. Because this is a linear transformation, the transformed variable will have the same correlation with the other variable as before.

Again, this, in it's simplest form, only lets you generate a pair of correlated variables (this could be scaled up, but gets ugly fast), and is certainly not the most convenient way to get the job done. In R, you would want to use ?mvrnorm in the MASS package, both because it's easier and because you can generate many variables with a given population correlation matrix. Nonetheless, I think it's worthwhile to have walked through this process to see how some basic principles play out in a simple way.

Answer 3

In general this not a simple thing to do, but I believe there are packages for multivariate normal variable generation (at least in R, see mvrnorm in the MASS package), where you just input a covariance matrix and a mean vector.

There is also one more "constructive" approach. Let's say we want to model a random vector $(X_1,X_2)$ and we have its distribution function $F(x_1,x_2)$. The first step is to get the marginal distribution function; i.e. integrate $F$ over all $x_2$: $$F_{X_1}(x_1)= \int_{-\infty}^{\infty} F(x_1,x_2)dx_2. $$ Then we find $F^{-1}_{X_1}$ - the inverse function of $F_{X_1}$ - and plug in a random variable $\xi_1$ which is uniformly distributed on the interval $[0,1]$. In this step we generate the first coordinate $\hat{x}_1=F^{-1}_{X_1}(\xi)$.

Now, since we have got one coordinate, we need to plug it in to our initial distribution function $F(x_1,x_2)$ and then get a conditional distribution function with condition $x_1=\hat{x}_1$: $$F(x_2 | X_1=\hat{x}_1)= \frac{F(\hat{x}_1,x_2)}{f_{X_1}(\hat{x}_1)}, $$ where $f_{X_1}$ is a probability density function of the marginal $X_1$ distribution; i.e. $F'_{X_1}(x_1)=f_{X_1}(x_1)$.

Then again you generate a uniformly distributed variable $\xi_2$ on $[0,1]$ (independent of $\xi_1$) and plug it in to the inverse of $F(x_2 | X_1=\hat{x}_1)$. Therefore you obtain $\hat{x}_2=(F(x_2 | X_1=\hat{x}_1))^{-1}(\xi)$; that is, $\hat x_2$ satisfies $F(\hat x_2 | X_1=\hat{x}_1) = \xi$. This method can be generalized to vectors with more dimensions, but its downside is that you have to calculate, analytically or numerically, many functions. The idea can be found in this article as well: http://www.econ-pol.unisi.it/dmq/pdf/DMQ_WP_34.pdf.

If you don't understand the meaning of plugging a uniform variable into an inverse probability distribution function, try to make a sketch of the univariate case and then remember what the geometric interpretation of the inverse function is.

Answer 4

If you are ready to give up efficiency, you can use a throw-away alogorithm. Its advantage is, that it allows for any kind of distributions (not only Gaussian).

Start by generating two uncorrelated sequences of random numbers $\{x_i\}_{i=1}^N$ and $\{y_i\}_{i=1}^N$ with any desired distributions. Let $C$ by the desired value of the correlation coefficient. Then do the following:

1) Compute correlation coefficient $c_{old}=corr(\{x_i\},\{y_i\})$

2) Generate two random munbers $n_1$ and $n_2: 1 \leq n_{1,2} \leq N$

3) Swap numbers $x_{n_1}$ and $x_{n_2}$

4) Compute new correlation $c_{new}=corr( \{x_i\},\{y_i\})$

5) If $|C-c_{new}| < |C-c_{old}|$ then keep the swap. Else undo the swap.

6) If $|C-c| < \epsilon$ stop, else goto 1)

Random swaps will not alter the marginal distribution of ${x_i}$.

Good luck!