1

I want to generate random samples from normal distribution such that :

$X \sim \mathcal N(u_1,s_1)$

$Y \sim \mathcal N(u_2,s_2)$

and $\mathrm{cor}(X,Y)=k$ {k is non zero}.

If k=0 then X and Y can be generated easily in R, like : 

X <- rnorm(n,mean, sd)
Y <- rnorm(n,mean, sd)

But I donot know how can I generate joint random samples from Normal distribution such that their correlation is not equal to zero.

It would be helpful if you provide procedure with R code.

Sven Hohenstein
  • 6,833
Neeraj
  • 2,320
  • See mvnorm in MASS library or https://cran.r-project.org/web/packages/mvtnorm/index.html , however as strictly software-related this kind of questions are off-topic on this site. – Tim Nov 10 '15 at 07:21
  • 1
    As a pure-R "what function should I call" this would be off topic. As a "what's a way to generate correlated multivariate normals (that I could implement in R)" question it would be on topic but has already been answered several times. The specific bivariate case has also been treated several times. – Glen_b Nov 10 '15 at 07:29
  • It is not pure R related question. I want to know the process in simple ways. I found similar question on [http://stats.stackexchange.com/questions/38856/how-to-generate-correlated-random-numbers-given-means-variances-and-degree-of] but it is difficult to understand as it involve matrix notation. I want answer in simple explanation. – Neeraj Nov 10 '15 at 07:33
  • 1
    For the multivariate case, the matrix version is the simple way to do it. Here's the bivariate case, which can be done easily without matrices. Does that cover what you need? – Glen_b Nov 10 '15 at 09:51

1 Answers1

1

Here I'm doing it the hard way, but if you go through it, it might be enlightening. The basic idea is I construct the covariance matrix $\Sigma$, then use the eigenvector decomposition to compute the matrix $\Sigma^{1/2}$. Then if $u = (u_1,u_2)$, we get the formula: $(X,Y) = u + \Sigma^{1/2}z$, where the $z$ are a pair of independent standard normal random variables.

For example:

s1 = 2
s2 = 4
u1 = 10
u2 = 12
r = 0.8
covar = r*s1*s2
Sigma = matrix(ncol=2,nrow=2,c(s1^2,covar,covar,s2^2))
temp = eigen(Sigma)
SqrtSigma = temp$vectors%*%diag(sqrt(temp$values))%*%t(temp$vectors)
XYvec = c(u1,u2) + SqrtSigma%*%rnorm(2)
XYvec
#           [,1]
# [1,]  9.265028
# [2,] 11.230126

I'll check it with a simulation:

x = rep(NA,1000)
y = rep(NA,1000)
for(i in 1:1000){
  XYvec = c(u1,u2) + SqrtSigma%*%rnorm(2)
  x[i] = XYvec[1]
  y[i] = XYvec[2]
}
var(x)
# [1] 4.060218
var(y)
# [1] 16.18099
cor(x,y)
# [1] 0.8002916

mean(x)
# [1] 9.935937
mean(y)
# [1] 11.94783
AlaskaRon
  • 2,277
  • I appreciate your answer. Can you use non matrix notation to simplify that? I have no idea about eigenvector decomposition. – Neeraj Nov 10 '15 at 10:26