I would like to generate pairs of random numbers with certain correlation. However, the usual approach of using a linear combination of two normal variables is not valid here, because a linear combination of uniform variables is not any more an uniformly distributed variable. I need the two variables to be uniform.

Any idea on how to generate pairs of uniform variables with a given correlation?

  • 7
    $\begingroup$ Closely related: stats.stackexchange.com/questions/30526. You also want to check out the copula tag--just click on the link here. A quick and dirty technique is to let $X$ be uniform$[0,1]$ and $Y=X$ when $X\le\alpha$ and $Y=1+\alpha-X$ otherwise. The correlation is $\rho=2(\alpha-1)^3+1$, whence $\alpha=1-((1-\rho)/2)^{1/3}$ does the trick. But copulas will give you more control... . $\endgroup$
    – whuber
    Commented Aug 5, 2013 at 20:32
  • $\begingroup$ Thanks for the comment, but yes, I think this method is really "dirty" $\endgroup$
    – Pythonist
    Commented Aug 6, 2013 at 12:46
  • 1
    $\begingroup$ My hope was that in seeing this approach you would recognize that you can (and ought to) provide additional criteria concerning the properties of your pairs of random numbers. If this is "dirty," then precisely what is wrong with the solution? Tell us so that we can provide more appropriate answers for your situation. $\endgroup$
    – whuber
    Commented Aug 6, 2013 at 14:13
  • $\begingroup$ This question was answered incidentally in the response to a closely related question: how to generate pairs of RVs with a linear regression relationship. Because the slope of the linear regression is related in a readily computed way to the correlation coefficient, and all possible slopes can be produced, it gives a way to produce exactly what you want. See stats.stackexchange.com/questions/257779/…. $\endgroup$
    – whuber
    Commented May 23, 2017 at 23:38
  • 1
    $\begingroup$ Please also see stats.stackexchange.com/questions/31771, which answers the generalization to three random uniforms. $\endgroup$
    – whuber
    Commented May 23, 2017 at 23:40

4 Answers 4


I'm not aware of a universal method to generate correlated random variables with any given marginal distributions. So, I'll propose an ad hoc method to generate pairs of uniformly distributed random variables with a given (Pearson) correlation. Without loss of generality, I assume that the desired marginal distribution is standard uniform (i.e., the support is $[0, 1]$).

The proposed approach relies on the following:
a) For standard uniform random variables $U_1$ and $U_2$ with respective distribution functions $F_1$ and $F_2$, we have $F_i(U_i) = U_i$, for $i = 1, 2$. Thus, by definition Spearman's rho is $$ \rho_{\rm S}(U_1, U_2) = {\rm corr}(F_1(U_1), F_2(U_2)) = {\rm corr}(U_1, U_2) . $$ So, Spearman's rho and Pearson's correlation coefficient are equal (sample versions might however differ).

b) If $X_1, X_2$ are random variables with continuous margins and Gaussian copula with (Pearson) correlation coefficient $\rho$, then Spearman's rho is $$ \rho_{\rm S}(X_1, X_2) = \frac{6}{\pi} \arcsin \left(\frac{\rho}{2}\right) . $$ This makes it easy to generate random variables that have a desired value of Spearman's rho.

The approach is to generate data from the Gaussian copula with an appropriate correlation coefficient $\rho$ such that the Spearman's rho corresponds to the desired correlation for the uniform random variables.

Simulation algorithm
Let $r$ denote the desired level of correlation, and $n$ the number of pairs to be generated. The algorithm is:

  1. Compute $\rho = 2\sin (r \pi/6)$.
  2. Generate a pair of random variables from the Gaussian copula (e.g., with this approach)
  3. Repeat step 2 $n$ times.

The following code is an example of implementation of this algorithm using R with a target correlation $r = 0.6$ and $n = 500$ pairs.

## Initialization and parameters 
r <- 0.6                            # Target (Spearman) correlation
n <- 500                            # Number of samples

## Functions
gen.gauss.cop <- function(r, n){
    rho <- 2 * sin(r * pi/6)        # Pearson correlation
    P <- toeplitz(c(1, rho))        # Correlation matrix
    d <- nrow(P)                    # Dimension
    ## Generate sample
    U <- pnorm(matrix(rnorm(n*d), ncol = d) %*% chol(P))

## Data generation and visualization
U <- gen.gauss.cop(r = r, n = n)
pairs(U, diag.panel = function(x){
          h <- hist(x, plot = FALSE)
          rect(head(h$breaks, -1), 0, tail(h$breaks, -1), h$counts/max(h$counts))})

In the figure below, the diagonal plots show histograms of variables $U_1$ and $U_2$, and off-diagonal plots show scatter plots of $U_1$ and $U_2$. enter image description here

By constuction, the random variables have uniform margins and a correlation coefficient (close to) $r$. But due to the effect of sampling, the correlation coefficient of the simulated data is not exactly equal to $r$.

cor(U)[1, 2]
# [1] 0.5337697

Note that the gen.gauss.cop function should work with more than two variables simply by specifying a larger correlation matrix.

Simulation study
The following simulation study repeated for target correlation $r= -0.5, 0.1, 0.6$ suggests that the distribution of the correlation coefficient converges to the desired correlation as the sample size $n$ increases.

## Simulation
r <- 0.6                                                # Target correlation
n <- c(10, 50, 100, 500, 1000, 5000); names(n) <- n     # Number of samples
S <- 1000                                               # Number of simulations

res <- sapply(n,
              function(n, r, S){
                   replicate(S, cor(gen.gauss.cop(r, n))[1, 2])
               r = r, S = S)
boxplot(res, xlab = "Sample size", ylab = "Correlation")
abline(h = r, col = "red")

enter image description here enter image description here enter image description here

  • 4
    $\begingroup$ The general method to generate correlated multivariate distributions with given marginal distributions is called a copula. $\endgroup$
    – whuber
    Commented Aug 6, 2013 at 2:25
  • $\begingroup$ @whuber, the use of copula allows to specify a dependence structure between random variables. The problem is that (Person) correlation is influenced by both the dependence structure and the margins. So, each choice of margins will require a corresponding choice of copula parameters, not to mention that some levels of correlation simply can't be attained for given margins (e.g., see here). If you are aware of a method that allows to 'control' the level of correlation for any choice of margins, I would love to know about it. $\endgroup$
    – QuantIbex
    Commented Aug 6, 2013 at 7:11
  • $\begingroup$ Thanks @QuantIbex. But I don't get why "a) implies that Spearman's rho and (Pearson's) correlation coefficient for random variables with standard uniform margins are approximately equal in large sample" $\endgroup$
    – Pythonist
    Commented Aug 6, 2013 at 12:45
  • 2
    $\begingroup$ Quantlbex, all you need is to create a continuous path of copulas from the lower to the upper Frechet-Hoeffding bounds. For identical marginals, the correlation coefficient will be a continuous function from that path into the interval $[-1,1]$. My "quick and dirty" example in a comment to the question is one such path, but obviously there are many others: copulas give you the fullest, most general way to create and describe such paths. What this shows is that the original question is (grossly) underdetermined: it ought to stipulate additional criteria for the solution. $\endgroup$
    – whuber
    Commented Aug 6, 2013 at 13:44
  • 1
    $\begingroup$ @Quantibex I took the liberty of adding a sentence that points out your gen.gauss.cop function will work for more than two variables with a (trivial) tweak. If you don't like the addition or wish to put it differently please revert or change as needed. $\endgroup$
    – Glen_b
    Commented May 24, 2017 at 0:26

From brillant answer to this topic by @QuantIbex some time has passed. Perhaps a fairly new module / library in R will allow someone to get the above a little easier. You can read the detailed documentation of the possibilities that offers simstudy. A definite advantage of this library is that the solution is generalized, that is, it allows for the determination of many variables with a given correlation.

# initiate variables
samples = 1000
number_of_items = 4
# get library
# make vector of means=0 for all variables; note, that for every variable it can be different, this case is just for simplicity

# for normal distribution params1 is mean, params2 is the variance 
# (that is, the square of the standard deviation); For most simulation 
# studies, we are interested in the value 1 (variance and deviation, 
# for a value of 1 it is the same). The vector mu + 1 will give the 
# vector as many values of 1 for the standard deviation as we want.
# so, code below generates distributions with a mean of 0 and 
# standard deviation of 1
df = genCorGen(samples, nvars = number_of_items, params1 = mu, 
params2 = (mu+1), dist = "normal", rho = .7, corstr = "cs", wide = TRUE)

# for uniform distribution params1 is minimum, params2 is maximum
df = genCorGen(samples, nvars = number_of_items, params1 = mu, 
params2 = (mu+1), dist = "uniform", rho = .7, corstr = "cs", wide = TRUE)
# this function generates a slightly unfriendly object that we will 
# simply convert into a "nicer" class
df = as.data.frame(df[,-1])

#one can check result with:

>           V1        V2        V3        V4
> V1 1.0000000 0.6896261 0.6916389 0.7098661
> V2 0.6896261 1.0000000 0.6858251 0.6965890
> V3 0.6916389 0.6858251 1.0000000 0.6755568
> V4 0.7098661 0.6965890 0.6755568 1.0000000


Of course, the library itself allows you to write better and shorter code, but I thought that a simple solution with a description might be useful to someone.


Intuitively, $u_1$ is $U(0,1)$ because $u_1$ equals $w_1$ [which is $U(0,1)$] if $I = 1$, and $u_1$ equals $w_2$ [which is $U(0,1)$] if $I = 0$, so $u_1$ is $U(0,1)$ in either case. The same for $u_2$. As for the correlation:

$E(u_1 u_2) = E[I w_1 + (1-I) w_2][I w_1 + (1-I) w_3]$

Expanding this, note first that $I(I-1)=0$, $I^2=I$, and $(1-I)^2=(1-I)$ because $I$ is always either $0$ or $1$. Note also that $I$ is independent of the $w$'s, which are also independent of each other. So:

$E(u_1 u_2) = E(I)E(w_1^2) + E(1-I)E(w_2)E(w_3)$ $=pE(w_1^2)+(1-p)/4$

From the fact that $V(w_1)=1/12$, we get $E(w_1^2)=1/3$, so $E(u_1 u_2) = p/12 + 1/4$, that is: $cov(u_1 u_2) = p/12$. Since $V(u_1)=V(u_2)=1/12$, we get finally that $cor(u_1, u_2) = p$.


Here is one easy method for positive correlation: Let $(u_1, u_2) = Iw_1 + (1-I) (w_2, w_3)$, where $w_1, w_2,$ and $w_3$ are independent $U(0,1)$ and $I$ is Bernoulli($p$). $u_1$ and $u_2$ will then have $U(0,1)$ distributions with correlation $p$. This extends immediately to $k$-tuples of uniforms with compound symmetric variance matrix.

If you want pairs with negative correlation, use $(u_1, u_2) = I(w_1, 1-w_1) + (1-I)(w_2, w_3)$, and the correlation will be $-p$.

  • $\begingroup$ Can you add a short proof of why this works? $\endgroup$ Commented Apr 28, 2018 at 13:30
  • 1
    $\begingroup$ if your want to be computationally efficient, $u_1=w_1$ also produces the same correlation (both positive and negative cases) $\endgroup$
    – Anvit
    Commented Apr 23, 2019 at 12:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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