Generate pairs of random numbers uniformly distributed and correlated

Question

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?

Answer 1

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.

Example
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 
set.seed(123)
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))
    return(U)
}

## 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$.

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
set.seed(921)
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")

Answer 2

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
library(simstudy)
# make vector of means=0 for all variables; note, that for every variable it can be different, this case is just for simplicity
mu<-rep(0,number_of_items)

# 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:
plot(df)
cor(df)

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

Answer 3

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

Answer 4

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