Generate a random variable with a defined correlation to an existing variable(s)

For a simulation study I have to generate random variables that show a prefined (population) correlation to an existing variable $Y$.

I looked into the R packages copula and CDVine which can produce random multivariate distributions with a given dependency structure. It is, however, not possible to fix one of the resulting variables to an existing variable.

Any ideas and links to existing functions are appreciated!

Conclusion: Two valid answers came up, with different solutions:

1. An R script by caracal, which calculates a random variable with an exact (sample) correlation to a predefined variable
2. An R function I found myself, which calculates a random variable with a defined population correlation to a predefined variable

[@ttnphns' addition: I took the liberty to expand the question title from single fixed variable case to arbitrary number of fixed variables; i.e. how to generate a variable having predefined corretation(s) with some fixed, existing variable(s)]

Here's another one: for vectors with mean 0, their correlation equals the cosine of their angle. So one way to find a vector $x$ with exactly the desired correlation $r$, corresponding to an angle $\theta$:

1. get fixed vector $x_1$ and a random vector $x_2$
2. center both vectors (mean 0), giving vectors $\dot{x}_{1}$, $\dot{x}_{2}$
3. make $\dot{x}_{2}$ orthogonal to $\dot{x}_{1}$ (projection onto orthogonal subspace), giving $\dot{x}_{2}^{\perp}$
4. scale $\dot{x}_{1}$ and $\dot{x}_{2}^{\perp}$ to length 1, giving $\bar{x}_{1}$ and $\bar{x}_{2}^{\perp}$
5. $\bar{x}_{2}^{\perp} + (1/\tan(\theta)) \cdot \bar{x}_{1}$ is the vector whose angle to $\bar{x}_{1}$ is $\theta$, and whose correlation with $\bar{x}_{1}$ thus is $r$. This is also the correlation to $x_1$ since linear transformations leave the correlation unchanged.

Here is the code:

n     <- 20                    # length of vector
rho   <- 0.6                   # desired correlation = cos(angle)
theta <- acos(rho)             # corresponding angle
x1    <- rnorm(n, 1, 1)        # fixed given data
x2    <- rnorm(n, 2, 0.5)      # new random data
X     <- cbind(x1, x2)         # matrix
Xctr  <- scale(X, center=TRUE, scale=FALSE)   # centered columns (mean 0)

Id   <- diag(n)                               # identity matrix
Q    <- qr.Q(qr(Xctr[ , 1, drop=FALSE]))      # QR-decomposition, just matrix Q
P    <- tcrossprod(Q)          # = Q Q'       # projection onto space defined by x1
x2o  <- (Id-P) %*% Xctr[ , 2]                 # x2ctr made orthogonal to x1ctr
Xc2  <- cbind(Xctr[ , 1], x2o)                # bind to matrix
Y    <- Xc2 %*% diag(1/sqrt(colSums(Xc2^2)))  # scale columns to length 1

x <- Y[ , 2] + (1 / tan(theta)) * Y[ , 1]     # final new vector
cor(x1, x)                                    # check correlation = rho For the orthogonal projection $P$, I used the $QR$-decomposition to improve numerical stability, since then simply $P = Q Q'$.

• I was trying to rewrite the code into SPSS syntax. I stumble over your QR decomposition which returns 20x1 column. In SPSS I have Gram-Schmidt orthonormalization (which is also a QR decomposition) but unable to replicate your resultant Q column. Can you chew over your QR action to me please. Or indicate some work-around to get the projection. Thanks. – ttnphns Aug 31 '11 at 21:28
• @caracal, P <- X %*% solve(t(X) %*% X) %*% t(X) doesn't produce r=0.6, so that's not the work-around. I'm still confused. (I'd be happy to mimic your expression Q <- qr.Q(qr(Xctr[ , 1, drop=FALSE])) in SPSS but don't know how.) – ttnphns Sep 3 '11 at 9:29
• @ttnphns Sorry for the confusion, my comment was for the general case. Applying it to situation in the example: Getting the projection matrix via QR-decomposition is just for numerical stability. You can get the projection matrix as $P=X(X'X)^{-1} X'$ if the subspace is spanned by the columns of matrix $X$. In R, you can here write Xctr[ , 1] %*% solve(t(Xctr[ , 1]) %*% Xctr[ , 1]) %*% t(Xctr[ , 1]) because the subspace is spanned by the first column of Xctr. The matrix for the projection onto the orthogonal complement is then I-P. – caracal Sep 4 '11 at 10:59
• Could anyone clarify how to perform something similar for a more than just two samples? Say, if I wanted 3 samples that are correlated pairwise by rho, how can I transform this solution to achieve that? – Andre Terra Jun 10 '16 at 17:24
• for the limit case rho=1 I found it useful to do something like this: if (isTRUE(all.equal(rho, 1))) rho <- 1-10*.Machine$double.eps, otherwise I was getting NaNs – PatrickT Nov 19 '17 at 14:33 I will describe the most general possible solution. Solving the problem in this generality allows us to achieve a remarkably compact software implementation: just two short lines of R code suffice. Pick a vector$X$, of the same length as$Y$, according to any distribution you like. Let$Y^\perp$be the residuals of the least squares regression of$X$against$Y$: this extracts the$Y$component from$X$. By adding back a suitable multiple of$Y$to$Y^\perp$, we may produce a vector having any desired correlation$\rho$with$Y$. Up to an arbitrary additive constant and positive multiplicative constant--which you are free to choose in any way--the solution is $$X_{Y;\rho} = \rho\, \operatorname{SD}(Y^\perp)Y + \sqrt{1-\rho^2}\,\operatorname{SD}(Y)Y^\perp.$$ ("$\operatorname{SD}$" stands for any calculation proportional to a standard deviation.) Here is working R code. If you don't supply$X$, the code will draw its values from the multivariate standard Normal distribution. complement <- function(y, rho, x) { if (missing(x)) x <- rnorm(length(y)) # Optional: supply a default if x is not given y.perp <- residuals(lm(x ~ y)) rho * sd(y.perp) * y + y.perp * sd(y) * sqrt(1 - rho^2) } To illustrate, I generated a random$Y$with$50$components and produced$X_{Y;\rho}$having various specified correlations with this$Y$. They were all created with the same starting vector$X=(1,2,\ldots, 50)$. Here are their scatterplots. The "rugplots" at the bottom of each panel show the common$Y$vector. There's a remarkable similarity among the plots, isn't there :-). If you would like to experiment, here is the code that produced these data and the figure. (I didn't bother to use the freedom to shift and scale the results, which are easy operations.) y <- rnorm(50, sd=10) x <- 1:50 # Optional rho <- seq(0, 1, length.out=6) * rep(c(-1,1), 3) X <- data.frame(z=as.vector(sapply(rho, function(rho) complement(y, rho, x))), rho=ordered(rep(signif(rho, 2), each=length(y))), y=rep(y, length(rho))) library(ggplot2) ggplot(X, aes(y,z, group=rho)) + geom_smooth(method="lm", color="Black") + geom_rug(sides="b") + geom_point(aes(fill=rho), alpha=1/2, shape=21) + facet_wrap(~ rho, scales="free") BTW, this method readily generalizes to more than one$Y$: if it's mathematically possible, it will find an$X_{Y_1,Y_2,\ldots,Y_k;\rho_1,\rho_2,\ldots,\rho_k}$having specified correlations with an entire set of$Y_i$. Just use ordinary least squares to take out the effects of all the$Y_i$from$X$and form a suitable linear combination of the$Y_i$and the residuals. (It helps to do this in terms of a dual basis for$Y$, which is obtained by computing a pseudo-inverse. The follownig code uses the SVD of$Y$to accomplish that.) Here's a sketch of the algorithm in R, where the$Y_i$are given as columns of a matrix y: y <- scale(y) # Makes computations simpler e <- residuals(lm(x ~ y)) # Take out the columns of matrix y y.dual <- with(svd(y), (n-1)*u %*% diag(ifelse(d > 0, 1/d, 0)) %*% t(v)) sigma2 <- c((1 - rho %*% cov(y.dual) %*% rho) / var(e)) return(y.dual %*% rho + sqrt(sigma2)*e) The following is a more complete implementation for those who would like to experiment. complement <- function(y, rho, x) { # # Process the arguments. # if(!is.matrix(y)) y <- matrix(y, ncol=1) if (missing(x)) x <- rnorm(n) d <- ncol(y) n <- nrow(y) y <- scale(y) # Makes computations simpler # # Remove the effects of y on x. # e <- residuals(lm(x ~ y)) # # Calculate the coefficient sigma of e so that the correlation of # y with the linear combination y.dual %*% rho + sigma*e is the desired # vector. # y.dual <- with(svd(y), (n-1)*u %*% diag(ifelse(d > 0, 1/d, 0)) %*% t(v)) sigma2 <- c((1 - rho %*% cov(y.dual) %*% rho) / var(e)) # # Return this linear combination. # if (sigma2 >= 0) { sigma <- sqrt(sigma2) z <- y.dual %*% rho + sigma*e } else { warning("Correlations are impossible.") z <- rep(0, n) } return(z) } # # Set up the problem. # d <- 3 # Number of given variables n <- 50 # Dimension of all vectors x <- 1:n # Optionally: specify x or draw from any distribution y <- matrix(rnorm(d*n), ncol=d) # Create d original variables in any way rho <- c(0.5, -0.5, 0) # Specify the correlations # # Verify the results. # z <- complement(y, rho, x) cbind('Actual correlations' = cor(cbind(z, y))[1,-1], 'Target correlations' = rho) # # Display them. # colnames(y) <- paste0("y.", 1:d) colnames(z) <- "z" pairs(cbind(z, y)) • This is indeed a nice solution. However, I failed to expand it myself to multiple$Y$variables (the fixed variables, in your answer) case. BTW, this method readily generalizes to more... Just use ordinary least squares... and form a suitable linear combination, you claim. Can you demonstrate it? Please, with annotated code readable by a not R user? – ttnphns Nov 11 '17 at 9:50 • @ttnphns I have done so. – whuber Nov 11 '17 at 17:41 • Thank you a lot! I see, and I've coded your approach today in SPSS for myself. Really great proposal of yours. I never thought of the notion of dual basis as applicable to solve the task. – ttnphns Nov 12 '17 at 11:04 • Is it possible to use a similar approach to come up with a uniformly distributed vector? That is, I have an existing vector x and want to generate a new vector y correlated with x but also want the y vector to be uniformly distributed. – Skumin Sep 18 '18 at 16:37 • @Skumin Consider using a copula for that so you can control the relationship between the two vectors. – whuber Sep 18 '18 at 17:08 Here's another computational approach (the solution is adapted from a forum post by Enrico Schumann). According to Wolfgang (see comments), this is computationally identical to the solution proposed by ttnphns. In contrast to caracal's solution it does not produce a sample with the exact correlation of$\rho$, but two vectors whose population correlation is equal to$\rho$. Following function can compute a bivariate sample distribution drawn from a population with a given$\rho$. It either computes two random variables, or it takes one existing variable (passed as parameter x) and creates a second variable with the desired correlation: # returns a data frame of two variables which correlate with a population correlation of rho # If desired, one of both variables can be fixed to an existing variable by specifying x getBiCop <- function(n, rho, mar.fun=rnorm, x = NULL, ...) { if (!is.null(x)) {X1 <- x} else {X1 <- mar.fun(n, ...)} if (!is.null(x) & length(x) != n) warning("Variable x does not have the same length as n!") C <- matrix(rho, nrow = 2, ncol = 2) diag(C) <- 1 C <- chol(C) X2 <- mar.fun(n) X <- cbind(X1,X2) # induce correlation (does not change X1) df <- X %*% C ## if desired: check results #all.equal(X1,X[,1]) #cor(X) return(df) } The function can also use non-normal marginal distributions by adjusting parameter mar.fun. Note, however, that fixing one variable only seems to work with a normally distributed variable x! (which might relate to Macro's comment). Also note that the "small correction factor" from the original post was removed as it seems to bias the resulting correlations, at least in the case of Gaussian distributions and Pearson correlations (also see comments). • It seems this is only an approximate solution, i.e., the empirical correlation is not exactly equal to$\rho$. Or am I missing something? – caracal Aug 31 '11 at 21:18 • It is easy to show that, except for that "small correction to rho" (whose purpose in this context eludes me), this is exactly the same as what ttnphns suggested earlier. The method is simply based on the Choleski decomposition of the correlation matrix to obtain the desired transformation matrix. See, for example: en.wikipedia.org/wiki/…. And yes, this will only give you two vectors whose population correlation is equal to rho. – Wolfgang Sep 1 '11 at 8:15 • The "small correction to rho" was in the original post and is described here. Actually, I don't really understand it; but an investigation of 50000 simulated correlations with rho = .3 shows that without the "small correction" an average of r's of .299 is produced, while with the correction an average of .312 (which is the value of the corrected rho) is produced. Therefore I removed that part from the function. – Felix S Sep 2 '11 at 6:53 • I know this is old, but I also want to note that this method won't work for non-positive definite correlation matrices. E.g - a correlation of -1. – zzk Oct 15 '13 at 1:07 • Thanks; I noticed that if x1 is not standardized mean=0, sd=1, and you'd rather not rescale it, you'll need to modify the line: X2 <- mar.fun(n) to X2 <- mar.fun(n,mean(x),sd(x)) to get the desired correlation between x1 and x2 – Dave M Mar 18 '16 at 3:55 Let$X$be your fixed variable and you want to generate$Y$variable that correlates with$X$by amount$r$. If$X$is standardized then (because$r$is beta coefficient in simple regression)$Y= rX+E$, where$E$is random variable from normal distribution having mean$0$and$\text{sd}=\sqrt{1-r^2}$. Observed correlation between$X$and$Y$data will be approximately$r$;$X$and$Y$can be seen as random samples from bivariate normal population (if$X$is from normal) with$\rho=r$. Now, if you want to attain the correlation in your bivariate sample exactly$r$, you need to provide that$E$has zero correlation with$X$. This tightening it to zero can be reached by modifying$E$iteratively. Well, with only two variables, one given ($X$) and one to generate ($Y$), the sufficient number of iterations is actually 1, but with multiple given variables ($X_1, X_2, X_3,...$) iterations will be needed. It should be noted that if$X$is normal then in the first procedure ("approximate$r$")$Y$will also be normal; however, in iterative fitting of$Y$to the "exact$r$"$Y$is likely to lose normality because the fitting exploits case values selectively. Update Nov 11, 2017. I've come across this old thread today and decided to expand my answer by showing the algorithm of the iterative fitting about which I was speaking initially. Here is an iterative solution how to train a randomly simulated or preexistent variable$Y$to correlate or covariate precisely as we desire (or very close to so - depending number of iterations) with a set of given variables$X$s (these cannot be modified). Disclamer: This iterative solution I've found inferior to the excellent one based on finding the dual basis and proposed by @whuber in this thread today. @whuber's solution is not iterative and, more importantly for me, it seems to be affecting the values of the input "pig" variable somewhat less than "my" algorithm (it'd be an asset then if the task is to "correct" the existing variable and not to generate random variate from scratch). Still, I'm publishing mine for curiosity and because it works (see also Footnote). So, we have given (fixed) variables$X_1, X_2,...,X_m$, and varible$Y$which is either just randomly generated "pig" of values or is an existent data variable which values we need to "correct" - to bring$Y$exactly to correlations (or it can be covariances)$r_1, r_2,...,r_m$with the$X$s. All data must be continuous; in other words, there should be a good deal of unique values. The idea: perform iterative fitting of residuals. Knowing the wanted (target) correlations/covariances, we may compute predicted values for the$Y$using the$X$s as multiple linear predictors. After obtaining the initial residuals (from the current$Y$and the ideal prediction), train them iteratively not to correlate with the predictors. In the end, regain$Y$with the residuals. (The procedure was my own experimental invention of the wheel many years ago when I knew none of the theory; I coded it then in SPSS.) 1. Convert the target$r$s to sums-of-crossproducts by multiplying them by$\text{df}=n-1$:$S_j=r_j \text{df}$. ($j$is a$X$variable index.) 2. Z-standardize all the variables (center each, then divide by the st. deviation computed on that above$\text{df}$).$Y$and$X$s are thus standard. Observed sums of squares are now =$\text{df}$. 3. Compute regressional coefficients predicting$Y$by$X$s according to the target$r$s:$\bf b=(X'X)^{-1} S$. 4. Compute predicted values for$Y$:$\hat{Y}=\bf Xb$. 5. Compute residuals$E=Y-\hat{Y}$. 6. Compute the needed (target) sum of squares for residuals:$SS_S=\text{df}-SS_{\hat {Y}}$. 7. (Begin to iterate.) Compute observed sums of crossproducts between current$E$and every$X_j$:$C_j= \sum_{i=1}^n E_i X_{ij}$8. Correct values of$E$in the aim to bring all$C$s closer to$0$($i$is a case index): $$E_i[\text{corrected}]=E_i-\frac{\sum_{j=1}^m C_j X_{ij}} {n\sum_{j=1}^m X_{ij}^2}$$ (the denominator doesn't change on iterations, compute it in advance) Or, alternatively, a more efficient formula additionally insures the mean of$E$becomes$0$. First, do center$E$at each iteration prior computation of the$C$s at step 7, then on this step 8 correct as: $$E_i[\text{corrected}]=E_i-\frac{\sum_{j=1}^m \frac{C_j X_{ij}^3}{\sum_{i=1}^n X_{ij}^2}} {\sum_{j=1}^m X_{ij}^2}$$ (again, denominators are known in advance)$^1$9. Bring$SS_E$to its target value:$E_i[\text{corrected}]=E_i \sqrt{SS_S/SS_E}$Go to step 7. (Do, say, 10-20 iterations; the greater is$m$the more iterations could be needed. If target$r$s were realistic,$SS_S$is positive, and if sample size$n$isn't too few, iterations always direct to convergence. End iterating.) 10. Ready: All the$C$s are almost zero now which means the residuals$E$has been trained to restore target$r$s. Compute the fitting$Y$:$Y[\text{corrected}]=\hat{Y}+E$. 11. The obtained$Y$is almost standardized. As a last stroke, you may want to standardize it precisely, again like you did it on step 2. 12. You may supply$Y$with any variance and mean you like. Actually, among the four statistics - min, max, mean, st. dev. - you may select any two values and linearly transform the variable so it posesses them without altering the$r$s (correlations) you've attained (it is all called linear rescaling). To warn again what was said above. With that pulling of$Y$exactly to the$r$, the output$Y$does not have to be normally distributed.$^1$The correction formula can be further sophisticated, for example, to insure greater homoscedasticity (in terms of sums-of-squares) of$Y$with every$X$as well, simultaneously with attaining the correlations, - I've implemented a code for that too. (I don't know if such "double" task is solvable via a more neat - noniterative - approach such as whuber's.) • Thanks for your answer. That is an empirical/ iterative solution I was thinking about as well. For my simulations, however, I need a more analytical solution without a costly fitting procedure. Fortunately, I just found a solution which I will post shortly ... – Felix S Aug 31 '11 at 11:56 • This works for generating bivariate normals but does not work for an arbitrary distribution (or any non-'additive' distribution) – Macro Aug 31 '11 at 11:57 • I don't see why you propose iteration when you can produce the entire cone of solutions directly. Is there some special purpose to this approach? – whuber Nov 10 '17 at 22:10 • Re your latest edit: Since I provide a simple formula for all solutions, one can achieve any desired objective such as "greater homoscedasticity" by minimizing an appropriate objective function over the set of all solutions. The approach is fully general. By extending the variable (or variables)$Y$to an orthogonal basis and exploiting the scale-invariance of correlation, the problem becomes one of optimizing a function defined on a sphere in a Euclidean space. – whuber Nov 19 '17 at 17:00 • @whuber, your comment is what I was waiting for; actually my answer (about heteroscedasticity, which I link to) was intended as a challenge for you: perhaps it is an invitation to post your solution - as thorough and brilliant as you usually do. – ttnphns Nov 19 '17 at 17:08 I felt like doing some programming, so I took @Adam's deleted answer and decided to write a nice implementation in R. I focus on using a functionally oriented style (i.e. lapply style looping). The general idea is to take two vectors, randomly permute one of the vectors until a certain correlation has been reached between them. This approach is very brute-force, but is simple to implement. First we create a function that randomly permutes the input vector: randomly_permute = function(vec) vec[sample.int(length(vec))] randomly_permute(1:100)  71 34 8 98 3 86 28 37 5 47 88 35 43 100 68 58 67 82  13 9 61 10 94 29 81 63 14 48 76 6 78 91 74 69 18 12  1 97 49 66 44 40 65 59 31 54 90 36 41 93 24 11 77 85  32 79 84 15 89 45 53 22 17 16 92 55 83 42 96 72 21 95  33 20 87 60 38 7 4 52 27 2 80 99 26 70 50 75 57 19  73 62 23 25 64 51 30 46 56 39 ...and create some example data vec1 = runif(100) vec2 = runif(100) ...write a function that permutes the input vector, and correlates it to a reference vector: permute_and_correlate = function(vec, reference_vec) { perm_vec = randomly_permute(vec) cor_value = cor(perm_vec, reference_vec) return(list(vec = perm_vec, cor = cor_value)) } permute_and_correlate(vec2, vec1)$vec
 0.79072381 0.23440845 0.35554970 0.95114398 0.77785348 0.74418811
 0.47871491 0.55981826 0.08801319 0.35698405 0.52140366 0.73996913
 0.67369873 0.85240338 0.57461506 0.14830718 0.40796732 0.67532970
 0.71901990 0.52031017 0.41357545 0.91780357 0.82437619 0.89799621
 0.07077250 0.12056045 0.46456652 0.21050067 0.30868672 0.55623242
 0.84776853 0.57217746 0.08626022 0.71740151 0.87959539 0.82931652
 0.93903143 0.74439384 0.25931398 0.99006038 0.08939812 0.69356590
 0.29254936 0.02674156 0.77182339 0.30047034 0.91790830 0.45862163
 0.27077191 0.74445997 0.34622648 0.58727094 0.92285322 0.83244284
 0.61397396 0.40616274 0.32203732 0.84003379 0.81109473 0.50573325
 0.86719899 0.45393971 0.19701975 0.63877904 0.11796154 0.26986325
 0.01581969 0.52571331 0.27087693 0.33821824 0.52590383 0.11261002
 0.89840404 0.82685046 0.83349287 0.46724807 0.15345334 0.60854785
 0.78854984 0.95770015 0.89193212 0.18885955 0.34303707 0.87332019
 0.08890968 0.22376395 0.02641979 0.43377516 0.58667068 0.22736077
 0.75948043 0.49734797 0.25235660 0.40125309 0.72147500 0.92423638
 0.27980561 0.71627101 0.07729027 0.05244047

$cor  0.1037542 ...and iterate a thousand times: n_iterations = lapply(1:1000, function(x) permute_and_correlate(vec2, vec1)) Note that R's scoping rules ensure that vec1 and vec2 are found in the global environment, outside the anonymous function used above. So, the permutations are all relative to the original test datasets we generated. Next, we find the maximum correlation: cor_values = sapply(n_iterations, '[[', 'cor') n_iterations[[which.max(cor_values)]]$vec
 0.89799621 0.67532970 0.46456652 0.75948043 0.30868672 0.83244284
 0.86719899 0.55623242 0.63877904 0.73996913 0.71901990 0.85240338
 0.81109473 0.52571331 0.82931652 0.60854785 0.19701975 0.26986325
 0.58667068 0.52140366 0.40796732 0.22736077 0.74445997 0.40125309
 0.89193212 0.52031017 0.92285322 0.91790830 0.91780357 0.49734797
 0.07729027 0.11796154 0.69356590 0.95770015 0.74418811 0.43377516
 0.55981826 0.93903143 0.30047034 0.84776853 0.32203732 0.25235660
 0.79072381 0.58727094 0.99006038 0.01581969 0.41357545 0.52590383
 0.27980561 0.50573325 0.92423638 0.11261002 0.89840404 0.15345334
 0.61397396 0.27077191 0.12056045 0.45862163 0.18885955 0.77785348
 0.23440845 0.05244047 0.25931398 0.57217746 0.35554970 0.34622648
 0.21050067 0.08890968 0.84003379 0.95114398 0.83349287 0.82437619
 0.46724807 0.02641979 0.71740151 0.74439384 0.14830718 0.82685046
 0.33821824 0.71627101 0.77182339 0.72147500 0.08801319 0.08626022
 0.87332019 0.34303707 0.45393971 0.47871491 0.29254936 0.08939812
 0.35698405 0.67369873 0.27087693 0.78854984 0.87959539 0.22376395
 0.02674156 0.07077250 0.57461506 0.40616274

$cor  0.3166681 ...or find the closest value to a correlation of 0.2: n_iterations[[which.min(abs(cor_values - 0.2))]]$vec
 0.02641979 0.49734797 0.32203732 0.95770015 0.82931652 0.52571331
 0.25931398 0.30047034 0.55981826 0.08801319 0.29254936 0.23440845
 0.12056045 0.89799621 0.57461506 0.99006038 0.27077191 0.08626022
 0.14830718 0.45393971 0.22376395 0.89840404 0.08890968 0.15345334
 0.87332019 0.92285322 0.50573325 0.40796732 0.91780357 0.57217746
 0.52590383 0.84003379 0.52031017 0.67532970 0.83244284 0.95114398
 0.81109473 0.35554970 0.92423638 0.83349287 0.34622648 0.18885955
 0.61397396 0.89193212 0.74445997 0.46724807 0.72147500 0.33821824
 0.71740151 0.75948043 0.52140366 0.69356590 0.41357545 0.21050067
 0.87959539 0.11796154 0.73996913 0.30868672 0.47871491 0.63877904
 0.22736077 0.40125309 0.02674156 0.26986325 0.43377516 0.07077250
 0.79072381 0.08939812 0.86719899 0.55623242 0.60854785 0.71627101
 0.40616274 0.35698405 0.67369873 0.82437619 0.27980561 0.77182339
 0.19701975 0.82685046 0.74418811 0.58667068 0.93903143 0.74439384
 0.46456652 0.85240338 0.34303707 0.45862163 0.91790830 0.84776853
 0.78854984 0.05244047 0.58727094 0.77785348 0.01581969 0.27087693
 0.07729027 0.71901990 0.25235660 0.11261002

$cor  0.2000199 To get a higher correlation, you need to increase the number of iterations. Let's solve a more general problem: given variable$Y_1$how to generate the random variables$Y_2,\dots,Y_n$with correlation matrix$R$? Solution: 1. get the cholesky decomposition of the correlation matrix$CC^T=R$2. create independent random vectors$X_2,\dots,X_n$of the same length as$Y_1$3. Use$Y_1$as the first column and append the generated randoms to it 4.$Y=CX$, where$Y_i$- the new random correlated numbers as required, note, that$Y_1$will not change Python code: import numpy as np import math from scipy.linalg import toeplitz, cholesky from statsmodels.stats.moment_helpers import cov2corr # create the large correlation matrix R p = 4 h = 2/p v = np.linspace(1,-1+h,p) R = cov2corr(toeplitz(v)) # create the first variable T = 1000; y = np.random.randn(T) # generate p-1 correlated randoms X = np.random.randn(T,p) X[:,0] = y C = cholesky(R) Y = np.matmul(X,C) # check that Y didn't change print(np.max(np.abs(Y[:,0]-y))) # check the correlation matrix print(R) print(np.corrcoef(np.transpose(Y))) Test Output: 0.0 [[ 1. 0.5 0. -0.5] [ 0.5 1. 0.5 0. ] [ 0. 0.5 1. 0.5] [-0.5 0. 0.5 1. ]] [[ 1. 0.50261766 0.02553882 -0.46259665] [ 0.50261766 1. 0.51162821 0.05748082] [ 0.02553882 0.51162821 1. 0.51403266] [-0.46259665 0.05748082 0.51403266 1. ]] • Could you clarify what "not that$Y_1$will not change" means? – whuber Nov 11 '17 at 20:24 • @whuber it was a typo – Aksakal Nov 12 '17 at 16:16 Generate normal variables with SAMPLING covariance matrix as given covsam <- function(nobs,covm, seed=1237) {; library (expm); # nons=number of observations, covm = given covariance matrix ; nvar <- ncol(covm); tot <- nvar*nobs; dat <- matrix(rnorm(tot), ncol=nvar); covmat <- cov(dat); a2 <- sqrtm(solve(covmat)); m2 <- sqrtm(covm); dat2 <- dat %*% a2 %*% m2 ; rc <- cov(dat2);}; cm <- matrix(c(1,0.5,0.1,0.5,1,0.5,0.1,0.5,1),ncol=3); cm; res <- covsam(10,cm) ; res; Generate normal variables with POPULATION covariance matrix as given covpop <- function(nobs,covm, seed=1237) {; library (expm); # nons=number of observations, covm = given covariance matrix; nvar <- ncol(covm); tot <- nvar*nobs; dat <- matrix(rnorm(tot), ncol=nvar); m2 <- sqrtm(covm); dat2 <- dat %*% m2; rc <- cov(dat2); }; cm <- matrix(c(1,0.5,0.1,0.5,1,0.5,0.1,0.5,1),ncol=3); cm; res <- covpop(10,cm); res • You need to learn to format the code in the answer! There is a specific option to mark text as code fragments, use it! – kjetil b halvorsen Feb 6 '17 at 10:43 Just create a random vector and sort until you get desired r. • In what situations would this be preferable to the above solutions? – Andy W Oct 5 '11 at 12:23 • A situation where a user wants a simple answer. I read a similar question on the r forum, and its the answer that was given. – Adam Oct 5 '11 at 14:22 • Unfortunately this solution is not only computationally inefficient and approximate, it will often fail altogether unless some analysis is first applied to determine an appropriate distribution for the "random vector.". I think there is merit to the underlying idea of just throwing some random numbers at the problem and randomly permuting them (not "sorting" them!) until an approximate$r\$ is attained (because this is quick and easy to program), but that idea is not clearly expressed in this short reply. – whuber Oct 5 '11 at 15:12
• If this answer was given on the r-help forum, I suspect it was either (a) ironic (ie, intended as a joke), or (b) offered by someone who isn't very statistically sophisticated. To put this more succinctly, this is a poor answer to the question. -1 – gung Sep 4 '13 at 18:52