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

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

• Not all correlations are possible for all $n$, due to the fact that the variates have to sum to 1. Do you have any more information you could share? There are some special cases which aren't too hard to work with... Jan 12, 2018 at 18:24
• I'm sorry. I don't understand the issue that you're raising. Could you explain in a little more detail? I said $n$ and $\rho$ for generality but $n = 100, 1000,$ etc. is fine and a few values of $\rho$ are fine.
– Eli
Jan 12, 2018 at 18:35
• Well, if someone had told you $n = 1000$ and $\rho=0.9$, then you have a problem, but not one that is amenable to mathematical analysis (I believe) :) If you can pick $\rho$, that's a different matter. Jan 12, 2018 at 18:37
• Sure you can pick the $\rho$! The world's your oyster :). This is all preliminary analysis to test out an idea, anyways. Are you aware of methods for specific choices of $\rho$?
– Eli
Jan 12, 2018 at 19:15
• In the specific case of $\rho = -1/(n-1)$, you can use a Dirichlet distribution with all the parameters $\alpha_i = 1$. I'll look in a couple of texts when I get home. It's an interesting question! Jan 12, 2018 at 19:38

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)
• Thank you! How are you arriving at that approximation?
– Eli
Jan 12, 2018 at 20:14
• After performing the numerical integration, I had to fit a function to the $(\rho,\tau)$ pairs it generated. OLS regression with some simple forms--linear combinations of monomials, exponentials, and trig functions--indicated a linear combination of just three sine functions worked well. I stopped there. You actually won't go too far wrong just using $\tau$ for the value of $\rho$: it will be accurate to within $0.015$.
– whuber
Jan 12, 2018 at 20:16
• Your insight is greatly appreciated! Thank you again.
– Eli
Jan 12, 2018 at 20:17
• I'm going through this in detail now, and I'm unsure of how you're finding $\rho$. Is the numeric integration you're referring to the integration of the numerator of $Cov(U_1, U_2) / (\sigma_{U_1} \sigma_{U_2})$
– Eli
Jan 12, 2018 at 22:28
• Yes. The denominator is $1/12$, of course.
– whuber
Jan 12, 2018 at 22:30

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.

• +1 Very nice analysis!
– whuber
Nov 21, 2023 at 15:44

@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