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?
up vote 41 down vote accepted

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.

  • Thanks! This was enormously helpful. I think I at least have a better sense of what I need to look at next. – Joseph Weissman Oct 8 '12 at 1:51
  • 7
    Does this method apply only for Gaussian distributions (as specified in the question), or can it be used for generating correlated variables that follow other distributions? If not, are you aware of a method that could be used in that case? – user000001 Mar 3 '13 at 15:03
  • 1
    @Michael: Yes. Having said that given $C$ is a valid covariance matrix the Cholesky decomposition is the fastest way. You could also get the (symmetric) square root $X$ matrix of $C$ by using SVD (so $C = XX = XX^T$, where $X = U S^{0.5} V^T$ from $C = USV^T$) but that would be more expensive too. – usεr11852 Apr 1 '16 at 18:29
  • 1
    @Michael: Of course. Their covariance will be (approximately) the same, not the numbers themselves. – usεr11852 Apr 3 '16 at 19:39
  • 1
    @Sid: Any continuous distribution not supported on the whole real line will immediately fail. For example if we use a uniform $U[0,1]$ we cannot guarantee that the "correlated numbers" will be in $[0,1]$; similarly for a Poisson we will end up with non-discrete numbers. In addition, any distribution where the sum of the distributions is not still the same distribution (eg. summing $t$-distribution does not result in $t$-distributions) will also fail. In all cases mentioned, the numbers produced will be correlated according the $C$ but they will not correspond to the distribution we started. – usεr11852 Jun 24 '17 at 12:55

+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.

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.

  • Smart idea! Has simple intuitive appeal. But yes seems expensive computationally. – MichaelChirico Sep 28 '16 at 16:27
  • (+1) very good point. It would be better at the start saying $f_{X,Y}(x,y)=f_X(x)\cdot f_{Y|X}(y)$, then it flows more natural to first generate one mariginal distribution and then the conditional distribution. Very excellent! – KevinKim Jan 18 '17 at 22:04

Your Answer

 
discard

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

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