# 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... – jbowman Jan 12 '18 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. – EliK Jan 12 '18 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. – jbowman Jan 12 '18 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$? – EliK Jan 12 '18 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! – jbowman Jan 12 '18 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? – EliK Jan 12 '18 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 '18 at 20:16
• Your insight is greatly appreciated! Thank you again. – EliK Jan 12 '18 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})$ – EliK Jan 12 '18 at 22:28
• Yes. The denominator is $1/12$, of course. – whuber Jan 12 '18 at 22:30

@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