Generating and visualising multivariate random numbers in R

This post will present the wonderful pairs.panels function of the psych package [1] that I discovered recently to visualise multivariate random numbers.

Here is a little example with a Gaussian copula and normal and log-normal marginal distributions. I use pairs.panels to illustrate the steps along the way.

I start with standardised multivariate normal random numbers:

library(psych)
library(MASS)
Sig <- matrix(c(1, -0.7, -.5,
                -0.7, 1, 0.6,
                -0.5, 0.6, 1), 
                nrow=3)
X <- mvrnorm(1000, mu=rep(0,3), Sigma = Sig, 
             empirical = TRUE)  
pairs.panels(X)

Next, I map the random figures into the interval [0,1] using the distribution function pnorm, so I end up with multivariate uniform random numbers:

U <- pnorm(X)
pairs.panels(U)

Finally, I transform the uniform numbers into the desired marginals:

Z <- cbind(
  A=qlnorm(U[,1], meanlog=-2.5, sdlog=0.25),
  B=qnorm(U[,2], mean=1.70, sd=0.1),
  C=qnorm(U[,3], mean=0.63,sd=0.08)
  )
pairs.panels(Z)

Those steps can actually be shorten with functions of the copula package [2].

Unfortunately, I struggled to install the copula package on my local Mac, running Mavericks. No CRAN binaries are currently available and gfortran is playing up on my system to install it from source. Thus, I quickly fired up an Ubuntu virtual machine on Amazon’s EC2 Cloud and installed R. Within 15 minutes I was back in business - that is actually pretty amazing.

library(copula)
myCop=normalCopula(param=c(-0.7,-.5,0.6), dim = 3, dispstr = "un")
myMvd <- mvdc(copula=myCop, margins=c("lnorm", "norm", "norm"),
              paramMargins=list(list(meanlog=-2.5, sdlog=0.25),
                                list(mean=1.70, sd=0.1), 
                                list(mean=0.63,sd=0.08)) )
Z2 <- rmvdc(myMvd, 1000)
colnames(Z2) <- c("A", "B", "C")
pairs.panels(Z2)

References

[1] Revelle, W. (2014) psych: Procedures for Personality and Psychological Research, Northwestern University, Evanston, Illinois, USA, http://CRAN.R-project.org/package=psych Version = 1.4.5.

[2] Marius Hofert, Ivan Kojadinovic, Martin Maechler and Jun Yan (2014). copula: Multivariate Dependence with Copulas. http://CRAN.R-project.org/package=copula Version = 0.999-10

Session Info Local

R version 3.1.0 (2014-04-10)
Platform: x86_64-apple-darwin13.1.0 (64-bit)

locale:
[1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8

attached base packages:
[1] stats graphics  grDevices utils datasets methods base     

other attached packages:
[1] psych_1.4.5 MASS_7.3-31

loaded via a namespace (and not attached):
[1] tools_3.1.0

Session Info EC2

R version 3.0.2 (2013-09-25)
Platform: x86_64-pc-linux-gnu (64-bit)

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats graphics  grDevices utils datasets  methods   base     

other attached packages:
[1] psych_1.4.5 copula_0.999-10

loaded via a namespace (and not attached):
[1] ADGofTest_0.3     grid_3.0.2        gsl_1.9-10        lattice_0.20-24  
[5] Matrix_1.1-2      mvtnorm_0.9-99992 pspline_1.0-16    stabledist_0.6-6 
[9] stats4_3.0.2