8
$\begingroup$

[On recent questions I was looking into generating random vectors in R, and I wanted to share that "research" as an independent Q&A on a specific point.]

Generating random data with correlation can be done using the Cholesky decomposition of the correlation matrix $C = LL^{T}$ here , as reflected on prior posts here and here.

The question that I want to address is how to use the Uniform distribution to generate correlated random numbers from different marginal distributions in R.

$\endgroup$
  • 2
    $\begingroup$ You appear to have rediscovered the Gaussian copula e.g. see related question here. There are many other copulas in popular use, but the Gaussian is quite convenient and can be quite suitable for some situations. $\endgroup$ – Glen_b May 4 '15 at 0:31
8
$\begingroup$

Since the question is

"how to use the Uniform distribution to generate correlated random numbers from different marginal distributions in $\mathbb{R}$"

and not only normal random variates, the above answer does not produce simulations with the intended correlation for an arbitrary pair of marginal distributions in $\mathbb{R}$.

The reason is that, for most cdfs $G_X$ and $G_Y$,$$\text{cor}(X,Y)\ne\text{cor}(G_X^{-1}(\Phi(X),G_Y^{-1}(\Phi(Y)),$$when$$(X,Y)\sim\text{N}_2(0,\Sigma),$$where $\Phi$ denotes the standard normal cdf.

To wit, here is a counter-example with an Exp(1) and a Gamma(.2,1) as my pair of marginal distributions in $\mathbb{R}$.

library(mvtnorm)
#correlated normals with correlation 0.7
x=rmvnorm(1e4,mean=c(0,0),sigma=matrix(c(1,.7,.7,1),ncol=2),meth="chol")
cor(x[,1],x[,2])
  [1] 0.704503
y=pnorm(x) #correlated uniforms
cor(y[,1],y[,2])
  [1] 0.6860069
#correlated Exp(1) and Ga(.2,1)
cor(-log(1-y[,1]),qgamma(y[,2],shape=.2))
  [1] 0.5840085

Another obvious counter-example is when $G_X$ is the Cauchy cdf, in which case the correlation is not defined.

To give a broader picture, here is an R code where both $G_X$ and $G_Y$ are arbitrary:

etacor=function(rho=0,nsim=1e4,fx=qnorm,fy=qnorm){
  #generate a bivariate correlated normal sample
  x1=rnorm(nsim);x2=rnorm(nsim)
  if (length(rho)==1){
    y=pnorm(cbind(x1,rho*x1+sqrt((1-rho^2))*x2))
    return(cor(fx(y[,1]),fy(y[,2])))
    }
  coeur=rho
  rho2=sqrt(1-rho^2)
  for (t in 1:length(rho)){
     y=pnorm(cbind(x1,rho[t]*x1+rho2[t]*x2))
     coeur[t]=cor(fx(y[,1]),fy(y[,2]))}
  return(coeur)
  }

enter image description here

Playing around with different cdfs led me to single out this special case of a $\chi^2_3$ distribution for $G_X$ and a log-Normal distribution for $G_Y$:

rhos=seq(-1,1,by=.01)
trancor=etacor(rho=rhos,fx=function(x){qchisq(x,df=3)},fy=qlnorm)
plot(rhos,trancor,ty="l",ylim=c(-1,1))
abline(a=0,b=1,lty=2)

which shows how far from the diagonal the correlation can be.

A final warning Given two arbitrary distributions $G_X$ and $G_Y$, the range of possible values of $\text{cor}(X,Y)$ is not necessarily $(-1,1)$. The problem may thus have no solution.

$\endgroup$
  • $\begingroup$ Fantastic! Ty! Is there any way that we can find an approximate segment where the departure is not to marked, like it seems to be the case with normals, to still be reasonable for practical applications? $\endgroup$ – Antoni Parellada May 4 '15 at 12:50
5
$\begingroup$

I wrote the correlate package. People said it is promising (worthy of a publish in Journal of Statistical Software), but I never wrote the paper for it because I chose not to pursue an academic career.

I believe the not maintained correlate package is still on CRAN.

When you install it, you can do the following:

require('correlate')
a <- rnorm(100)
b <- runif(100)
newdata <- correlate(cbind(a,b),0.5)

The result is that newdata will have a correlation of 0.5, without changing the univariate distributions of a and b (same values are there, they just get moved around until the multivariate 0.5 correlation has been reached.

I'll reply on questions here, sorry for the lack of documentation.

$\endgroup$
  • $\begingroup$ Bravo, this is the perfect answer! Do you have a way to detect values of the correlation that are impossible to reach? $\endgroup$ – Xi'an May 11 '15 at 17:00
  • $\begingroup$ @Xi'an There are some impossibilities, like few data points and a really specific correlation sought that just cannot be reached. e.g. only having 3 paired values. $\endgroup$ – PascalVKooten May 11 '15 at 18:45
  • $\begingroup$ Also note it is possible for more than 2 variables, e.g. for 3 variables you can define a 3x3 correlation matrix, 4 variables a 4x4. $\endgroup$ – PascalVKooten May 11 '15 at 18:47
  • $\begingroup$ Generally it will work as long as you don't want the impossible, but before you do serious work with it it is advised to do a couple of test runs. $\endgroup$ – PascalVKooten May 11 '15 at 18:49
  • $\begingroup$ People who were interested in it were using income data; loads of zeros and a gaussian-ish distribution for non-zero incomes. $\endgroup$ – PascalVKooten May 11 '15 at 18:49
1
$\begingroup$
  1. Generate two samples of correlated data from a standard normal random distribution following a predetermined correlation.

    As an example, let's pick a correlation r = 0.7, and code a correlation matrix such as:

    (C <- matrix(c(1,0.7,0.7,1), nrow = 2)) [,1] [,2] [1,] 1.0 0.7 [2,] 0.7 1.0

    We can use mvtnorm to generate now these two samples as a bivariate random vector:

    set.seed(0)

    SN <- rmvnorm(mean = c(0,0), sig = C, n = 1e5) resulting in two vector components distributed as ~ $N(0, 1)$ and with a cor(SN[,1],SN[,2])= 0.6996197 ~ 0.7. Both components can be extricated as follows:

    X1 <- SN[,1]; X2 <- SN[,2]

    Here's the plot with the overlapping regression line:

  2. Use the Probability Integral Transform here to obtain a bivariate random vector with marginal distributions ~ $U(0, 1)$ and the same correlation:

    U <- pnorm(SN) - so we are feeding into pnorm the SN vector to find $erf(SN)$ (or $\Phi(SN)$). In the process, we preserve the cor(U[,1], U[,2]) = 0.6816123 ~ 0.7 .

    Again we can decompose the vector U1 <- U[,1]; U2 <- U[,2] and produce a scatterplot with marginal distributions at the edges, clearly showing their uniform nature:

  3. Apply the inverse transform sampling method here to finally obtain the bivector of equally correlated points belonging to whichever distribution family we set out to reproduce.

    From here we can just generate two vectors distributed normally and with equal or different variances. For instance: Y1 <- qnorm(U1, mean = 8,sd = 10) and Y2 <- qnorm(U2, mean = -5, sd = 4), which will maintain the desired correlation, cor(Y1,Y2) = 0.6996197 ~ 0.7.

    Or opt for different distributions. If the distributions chosen are very dissimilar, the correlation may not be as precise. For instance, let's get U1 to follow a $t$ distribution with 3 d.f., and U2 an exponential with a $\lambda$=1: Z1 <- qt(U1, df = 3) and Z2 <- qexp(U2, rate = 1) The cor(Z1,Z2) [1] 0.5941299 < 0.7. Here are the respective histograms:

Here is an example of code for the entire process and normal marginals:

Cor_samples <- function(r, n, mean1, mean2, sd1, sd2){
C <- matrix(c(1,r,r,1), nrow = 2)
require(mvtnorm)
SN <- rmvnorm(mean = c(0,0), sig = C, n = n)
U <- pnorm(SN)
U1 <- U[,1]
U2 <- U[,2]

 Y1 <<- qnorm(U1, mean = mean1,sd = sd1) 
 Y2 <<- qnorm(U2, mean = mean2,sd = sd2) 

sample_measures <<- as.data.frame(c(mean(Y1), mean(Y2), sd(Y1), sd(Y2), cor(Y1,Y2)), names<-c("mean Y1", "mean Y2", "SD Y1", "SD Y2", "Cor(Y1,Y2)"))
sample_measures
}

For comparison, I've put together a function based on the Cholesky decomposition:

Cholesky_samples <- function(r, n, mean1, mean2, sd1, sd2){
C <- matrix(c(1,r,r,1), nrow = 2)
L <- chol(C)
X1 <- rnorm(n)
X2 <- rnorm(n)
X <- rbind(X1,X2)

Y <- t(L)%*%X
Y1 <- Y[1,]
Y2 <- Y[2,]

N_1 <<- Y[1,] * sd1 + mean1
N_2 <<- Y[2,] * sd2 + mean2

sample_measures <<- as.data.frame(c(mean(N_1), mean(N_2), sd(N_1), sd(N_2), cor(N_1, N_2)), 
                  names<-c("mean N_1", "mean N_2", "SD N_1", "SD N_2","cor(N_1,N_2)"))
sample_measures
}

Trying both methods to generate correlated (say, $r=0.7$) samples distributed ~ $N(97,23)$ and $N(32,8)$ we get, setting set.seed(99):

Using the Uniform:

cor_samples(0.7, 1000, 97, 32, 23, 8)
           c(mean(Y1), mean(Y2), sd(Y1), sd(Y2), cor(Y1, Y2))
mean Y1                                            96.5298821
mean Y2                                            32.1548306
SD Y1                                              22.8669448
SD Y2                                               8.1150780
cor(Y1,Y2)                                          0.7061308

and Using the Cholesky:

Cholesky_samples(0.7, 1000, 97, 32, 23, 8)
             c(mean(N_1), mean(N_2), sd(N_1), sd(N_2), cor(N_1, N_2))
mean N_1                                                   96.4457504
mean N_2                                                   31.9979675
SD N_1                                                     23.5255419
SD N_2                                                      8.1459100
cor(N_1,N_2)                                                0.7282176
$\endgroup$
  • $\begingroup$ Empirically, it seems like when you go from N(0,1) -> $$F^{-1}(X)$$ ~ Unif. -> $$f(F^{-1}(X))$$ ~ distributed according to the chosen distributions, the correlation does not change unless the last distribution is substantially different from the initial N(0,1). I included the values... In any case, do you see specific issues with the method itself for practical application? $\endgroup$ – Antoni Parellada May 3 '15 at 20:50
  • $\begingroup$ I changed the function at the end of the answer to include the correlation of the computed samples, so as to compare to the plugged-in number, and they seem to match. $\endgroup$ – Antoni Parellada May 3 '15 at 21:03
  • 2
    $\begingroup$ Whether there are issues with practical application depends on the practical application; for some things this is okay. Note that since the transformations are monotonic, nonparametric correlations like Spearman's rho and Kendall's tau won't be changed. $\endgroup$ – Glen_b May 4 '15 at 0:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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