# MASS for the multivariate normal RNG
require(MASS)
## Loading required package: MASS
# car package for scatterplot matrix
require(car)
## Loading required package: car
## Warning: package 'car' was built under R version 3.5.1
## Loading required package: carData
## Warning: package 'carData' was built under R version 3.5.1
# generate two correlated random normal variables (rho=.6)
d <- as.data.frame(mvrnorm(10000, mu=c(0,0), Sigma=matrix(c(1,.6,.6,1),2,2)))
names(d) <- c("X1","X2")

# view the bivariate distribution
scatterplotMatrix(~X1+X2, d, col=rgb(0,0,0,.1), regLine=F,
                  diagonal=list(method ="histogram", breaks="FD"))

# make them marginally uniform (retain correlation)
# by transforming with the cumulative distribution function of the
# appropriate normal distribution
d$X1 <- with(d, pnorm(X1, mean(X1), sd(X1)))
d$X2 <- with(d, pnorm(X2, mean(X2), sd(X2)))

# plot again to see that the marginals are now uniform
# but the correlation is still present
# one way to describe what has happened is that we separated the
# marginal disrtibution from the multivariate dependency
scatterplotMatrix(~X1+X2, d, col=rgb(0,0,0,.1), regLine=F,
                  diagonal=list(method ="histogram", breaks="FD"))

## now transform the maginals to the desired distribution
# poisson(2) (arbitrary choice)
d$X1 <- qpois(d$X1, lambda = 2)
# binomial(1, .5)
d$X2 <- qbinom(d$X2, 1, .5)

# the correlation is still there
# however, it's smaller, which seems related to conversion to binomial
# converting to other more continuous distributions results in correlations 
# similar to those generated in the normal
cor(d)
##           X1        X2
## X1 1.0000000 0.4574196
## X2 0.4574196 1.0000000
# look at the marginal distributions
par(mfrow=c(2,1))
barplot(table(d$X1), main="", xlab=expression(X[1]))
barplot(table(d$X2), main="", xlab=expression(X[2]))

# visualize the correlation
# (use some jitter and transparency, since these are discrete)
par(mfrow=c(1,1))
plot(jitter(d$X1),jitter(d$X2), pch=20, col=rgb(0,0,0,.1), 
     xlab=expression(X[1]), ylab=expression(X[2]))