How to simulate n Uniform[0,1] variables with a specified correlation?

Question

For a simulation study, I need to generate $n$ Uniform$[0,1]$ random variables with population correlation $\rho$. I'm not aware of any simple ways to do this.

I have considered sampling from a multivariate normal with mean 0 and exchangeable correlation matrix $ \Sigma= \left[ {\begin{array}{cccc} 1 & \rho &\dots & \rho \\ \rho & 1 & \dots & \rho \\ \vdots & \vdots & \vdots & 1 \end{array} } \right] \\ $ then using an inverse normal transformation, but the correlation will not be $\rho$.

Answer 1

Your method works.

Take a standard bivariate normal random variable $(X_1,X_2)$ with correlation $\tau$. Let $\Phi$ be the standard Normal distribution function. The copula $(U_1,U_2)$ given by $U_i = \Phi(X_i)$ has some correlation $\rho=f(\tau)$. Although this function $f$ cannot be determined analytically (as far as I can tell), numerical integration with a full suite of values $\tau$ spanning the interval $[-1,1]$ suggests its inverse $f^{-1}$, which determines what $\tau$ needs to be to yield a correlation $\rho$, can be approximated to about $0.000001$ throughout the interval by the function

$$\hat f^{-1}(\rho) = \frac{4118}{3163} \sin(\rho/3) + \frac{3183}{3149} \sin(2 \rho/3) - \frac{145}{2391} \sin(3 \rho/3).$$

Tests with the following R implementation support the claimed accuracy. (Actually, it would take a large simulation--greater than $10^{12}$ observations--to detect the error of approximation.)

Edit

Let $f_\tau$ be the standard bivariate Normal density,

$$f_\tau(x_1,x_2) = \frac{1}{2\pi \sqrt{1-\tau^2}} \exp\left(-\frac{1}{2(1-\tau^2)}(x_1^2+x_2^2-2\tau x_1x_2)\right).$$ Then the correlation of the $U_i$ is a function of the moments

$$m_\tau(i,j) = E[U_1^i U_2^j] = E[\Phi(X_1)^i\Phi(X_2)^j] = \iint f_\tau(x_1,x_2) \Phi(x_1)^i \Phi(x_2)^j dx_1 dx_2;$$

$$f(\tau)=\operatorname{Cor}(U_1,U_2) = \frac{m_\tau(1,1) - m_\tau(1,0)^2}{m_\tau(2,0)-m_\tau(1,0)^2}.$$

Most of these have easy analytical expressions based on univariate moments of the uniform distribution, because they do not depend on $\rho$: $m_\rho(1,0)=1/2,$ $m_\rho(2,0)=1/3.$ The one I integrated numerically is $m_\rho(1,1)$. I enforced the symmetry $f(-\tau)=-f(\tau)$ by numerically computing both $f(\pm \tau)$ and averaging their absolute values.

---

unif.to.norm <- function(rho) {
  4118/3163 * sin(rho/3) + 3183/3149 * sin(2 * rho/3) - 145/2391 * sin(rho)
}
# curve(unif.to.norm(rho)-rho, 0, 1, xname="rho") # Compare to a linear function
library(MASS) # mvrnorm
n <- 1e5      # observations to simulate
d <- 5        # dimensions
rho <- 0.6    # Must lie in [-1/(d-1), 1]
#
# Compute the correlation matrix for the multivariate Normal distribution.
#
rho.norm <- unif.to.norm(rho)
Sigma <- matrix(rho.norm, d, d) + diag(rep(1-rho.norm, d))
#
# Simulate data.
#
x <- pnorm(mvrnorm(n, rep(0, d), Sigma))
#
# Display the data.
# 
# pairs(x)               # The scatterplots
# par(mfrow=c(1, d))
# apply(x, 2, hist)      # The histograms (will be uniform)
# par(mfrow=c(1,1))
Sigma.hat <- cor(x)
prec <- ceiling(log10(n)/2)
signif(Sigma.hat, prec)  # Should almost match `rho` in all off-diagonal entries
signif(mean(Sigma.hat[lower.tri(Sigma.hat)]), prec+1)

Answer 2

Inspired by @whuber 's approach, I created my own approximation that I hope might be useful to future readers. The essential insight is that the function $f^{-1}(\rho)$ (which allows one to determine the appropriate value of $\tau$ to use for sampling in order to obtain a correlation $\rho$) must have $f^{-1}(-1) = -1$, $f^{-1}(0) = 0$, and $f^{-1}(1) = 1$. Additionally, it seems to be almost the identity (i.e., $|f^{-1}(\rho) - \rho| \lessapprox 0.02$. Therefore, it is appropriate to fit the function $f^{-1}(\rho) - \rho$ using a model which is zero at -1, 0, and 1. A cubic model works well: $$ f^{-1}(\rho) - \rho \approx a \rho (1+\rho)(1-\rho), $$ where $a$ is the only fit parameter. Fitting the model to numerical integration results using least squares gives: $$ a^\star = 0.0469815. $$ This is very close to the rational number $7/149 = 0.0469799$. On the interval $[-1,1]$, this approximation does not deviate from the actual curve by more than $6 \times 10^{-5}$. Below is a plot of the approximation and the curve computed using numerical integration.

Answer 3

@whuber: Great derivation, and highly useful for those creating copulas. I always wondered where some correlation gets lost when channeling the normal random deviates with desired correlation through the CDF:

## we arrive at rho = 1.91
set.seed(123)
X <- mvrnorm(100000, mu = c(0, 0, 0), Sigma = matrix(c(1, 0.2, 0, 0.2, 1, 0, 0, 0, 1), nrow = 3), empirical = T)
print(cor(X))
              [,1]         [,2]          [,3]
[1,]  1.000000e+00 2.000000e-01 -3.053512e-15
[2,]  2.000000e-01 1.000000e+00  2.217207e-15
[3,] -3.053512e-15 2.217207e-15  1.000000e+00
P <- pnorm(X)
print(cor(P))
             [,1]         [,2]         [,3]
[1,] 1.0000000000 0.1915424223 0.0004832619
[2,] 0.1915424223 1.0000000000 0.0009725174
[3,] 0.0004832619 0.0009725174 1.0000000000

## need to set rho a bit higher to obtain desired rho = 0.2
set.seed(123)
X <- mvrnorm(100000, mu = c(0, 0, 0), Sigma = matrix(c(1, 0.209, 0, 0.209, 1, 0, 0, 0, 1), nrow = 3), empirical = T)
print(cor(X))
             [,1]         [,2]          [,3]
[1,] 1.000000e+00  2.09000e-01  2.200462e-15
[2,] 2.090000e-01  1.00000e+00 -3.040320e-15
[3,] 2.200462e-15 -3.04032e-15  1.000000e+00
P <- pnorm(X)
print(cor(P))  
             [,1]         [,2]         [,3]
[1,] 1.0000000000 0.2001865963 0.0009690459
[2,] 0.2001865963 1.0000000000 0.0004781438
[3,] 0.0009690459 0.0004781438 1.0000000000