Building a Multivariate Distribution
+As stated earlier, to generate multivariate data, we need a list of marginals (and their parameters), and a correlation structure (matrix). The marginal distributions can be built up as a list of lists, where each sublist contains the information for the target distribution.
+margins = list(
+ list("norm", mean = 3.14, sd = 0.1),
+ list("beta", shape1 = 1, shape2 = 4),
+ list("nbinom", size = 10, prob = 0.75)
+)The things to point out here are that in each sublist (marginal), the first item is an unnamed character string with the R name of the distribution without a letter prefix. E.g. instead of rnorm, we pass in just "norm". The second thing to note is that the remaining items are named arguments that go along with the distribution. A full list of built-in distributions is found in the appendix.
The next step is to define a correlation structure for the multivariate distribution. This correlation matrix can either come from observed data, or we can set it ourselves, or we can generate a random correlation matrix via bigsimr::rcor. Let’s create a simple correlation matrix where all off-diagonal elements are 0.5. Since we have 3 marginals, we need a \(3\times 3\) matrix.
rho <- matrix(0.5, nrow = 3, ncol = 3)
+diag(rho) <- 1.0
+rho
+#> [,1] [,2] [,3]
+#> [1,] 1.0 0.5 0.5
+#> [2,] 0.5 1.0 0.5
+#> [3,] 0.5 0.5 1.0Finally we can generate a random vector with our specified marginals and correlation structure. The last argument, type, is looking to know what kind of correlation matrix it is receiving. Right now it can handle Pearson, Spearman, or Kendal.
On my machine, there is no dedicated GPU, so I would see the following warning message once per session.
+#> Warning: warning.warn('No GPU/TPU found, falling back to CPU.')Taking a look at our random vector, we see that it hs 10 rows and 3 columns, one column for each marginal.
+x
+#> [,1] [,2] [,3]
+#> [1,] 3.150112 0.04986300 1
+#> [2,] 3.198523 0.12602385 4
+#> [3,] 3.092609 0.12625169 0
+#> [4,] 2.914857 0.02925966 1
+#> [5,] 3.213697 0.29187930 5
+#> [6,] 3.056378 0.22101881 2
+#> [7,] 3.131806 0.31518815 7
+#> [8,] 3.245636 0.14441806 3
+#> [9,] 3.358278 0.20572350 3
+#> [10,] 3.142494 0.39363244 3We can simulate many more samples and then check the histogram of each margin, as well as the estimated correlation between the columns.
+x = rvec(10000, rho = rho, params = margins, type = "pearson")
+
+par(mfrow=c(1,3))
+hist(x[,1], breaks = 30, xlab = "", main = "Normal")
+hist(x[,2], breaks = 30, xlab = "", main = "Beta")
+hist(x[,3], breaks = 30, xlab = "", main = "Negative Binomial")
cor(x)
+#> [,1] [,2] [,3]
+#> [1,] 1.0000000 0.4847622 0.4880023
+#> [2,] 0.4847622 1.0000000 0.4892170
+#> [3,] 0.4880023 0.4892170 1.0000000We can see that even wih 10,000 samples, the estimated correlation of the simulated data is not exactly the same as the target correlation. This can be explained by the fact that some correlations are simply not possible due to the discrete nature of certain distributions. Another possibility is that the copula algorith is biased and needs correction.
+