library(MASS)
# set correlation coefficient and make matrix
r <- 0.9
(corr_matrix <- matrix(c(1, r, r, 1), nrow = 2))
## [,1] [,2]
## [1,] 1.0 0.9
## [2,] 0.9 1.0
# generate multivariate normal distribution
tmp <- mvrnorm(
n = 5000,
mu = c(0,0),
Sigma = corr_matrix
)
# get probability that random variable with N(0,1) is below the value given
dat <- pnorm(tmp)
# check that the variables are uniformly distributed between 0 and 1
par(mfrow = c(1,2))
hist(dat[,1])
hist(dat[,2])

# check that the variables are correlated
par(mfrow = c(1,1))
plot(dat[,1], dat[,2], pch = 20, cex = 0.05)
