Generating and visualising multivariate random numbers in R

[This article was first published on mages' blog, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

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

To leave a comment for the author, please follow the link and comment on their blog: mages' blog.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)