Generating data from multivariate distributions

Besides learning from real data, SynthesizerPlus can generate data from fully specified multivariate distributions. This is useful for simulation studies, for teaching, and for benchmarking synthesizers against data with a known dependence structure. All generators accept a seed argument and need no additional packages.

Correlation matrices

make_corr() builds common correlation structures:

make_corr(4, "ar1", rho = 0.6)
#>       [,1] [,2] [,3]  [,4]
#> [1,] 1.000 0.60 0.36 0.216
#> [2,] 0.600 1.00 0.60 0.360
#> [3,] 0.360 0.60 1.00 0.600
#> [4,] 0.216 0.36 0.60 1.000
make_corr(3, "exchangeable", rho = 0.3, names = c("x", "y", "z"))
#>     x   y   z
#> x 1.0 0.3 0.3
#> y 0.3 1.0 0.3
#> z 0.3 0.3 1.0
R <- make_corr(5, "random")
round(R, 2)
#>       [,1] [,2]  [,3]  [,4] [,5]
#> [1,]  1.00 0.10 -0.54 -0.51 0.12
#> [2,]  0.10 1.00  0.16  0.26 0.10
#> [3,] -0.54 0.16  1.00  0.33 0.06
#> [4,] -0.51 0.26  0.33  1.00 0.10
#> [5,]  0.12 0.10  0.06  0.10 1.00

Classical multivariate distributions

Function	Distribution
`r_mvnorm()`	multivariate normal
`r_mvt()`	multivariate Student-t
`r_mvlnorm()`	multivariate log-normal
`r_mvskewnorm()`	multivariate skew-normal (Azzalini & Dalla Valle, 1996)
`r_dirichlet()`	Dirichlet (compositional data)
`r_mvmixture()`	finite mixture of multivariate normals

S <- make_corr(3, "ar1", 0.7)
x <- r_mvnorm(2000, mean = c(0, 5, 10), sigma = S, seed = 1)
round(colMeans(x), 2)
#>    X1    X2    X3 
#> -0.01  5.00  9.99
round(cor(x), 2)
#>      X1   X2   X3
#> X1 1.00 0.71 0.53
#> X2 0.71 1.00 0.70
#> X3 0.53 0.70 1.00

t3 <- r_mvt(2000, df = 3, sigma = S, seed = 1)
c(normal = max(abs(x[, 1])), t3 = max(abs(t3[, 1])))
#>    normal        t3 
#>  3.810277 18.231054

The skew-normal adds asymmetry through the shape vector alpha:

sn <- r_mvskewnorm(2000, omega = make_corr(2, "exchangeable", 0.5),
                   alpha = c(8, -3), seed = 1)
plot_pairs(as.data.frame(sn))

Mixtures produce multimodal data, e.g. latent clusters:

mix <- r_mvmixture(1500, weights = c(0.5, 0.3, 0.2),
                   means = list(c(0, 0), c(4, 4), c(-3, 5)),
                   sigmas = list(diag(2), make_corr(2, "exchangeable", 0.8),
                                 diag(c(0.5, 2))),
                   seed = 1)
plot_pairs(as.data.frame(mix))

Compositional data (shares that sum to one):

shares <- r_dirichlet(5, alpha = c(agriculture = 2, industry = 5, services = 10), seed = 1)
round(shares, 3)
#>      agriculture industry services
#> [1,]       0.053    0.397    0.550
#> [2,]       0.191    0.310    0.500
#> [3,]       0.242    0.271    0.486
#> [4,]       0.119    0.172    0.709
#> [5,]       0.209    0.176    0.614

Copulas

r_copula() samples from Gaussian, t, Clayton, Gumbel and Frank copulas (uniform margins). Archimedean copulas in more than two dimensions use the Marshall–Olkin frailty construction; their Kendall’s tau is a known function of theta:

Family	Kendall’s tau	Tail dependence
Clayton	\(\theta/(\theta+2)\)	lower
Gumbel	\(1 - 1/\theta\)	upper
Frank	(no closed form)	none

u_c <- r_copula(3000, "clayton", dim = 2, theta = 2, seed = 1)
u_g <- r_copula(3000, "gumbel", dim = 2, theta = 2, seed = 1)
c(clayton = cor(u_c, method = "kendall")[1, 2],
  gumbel = cor(u_g, method = "kendall")[1, 2])
#>   clayton    gumbel 
#> 0.4973071 0.4991010

Both have tau = 0.5, but extremes cluster in opposite tails:

both <- rbind(data.frame(u1 = u_c[, 1], u2 = u_c[, 2], copula = "Clayton (theta = 2)"),
              data.frame(u1 = u_g[, 1], u2 = u_g[, 2], copula = "Gumbel (theta = 2)"))
ggplot2::ggplot(both, ggplot2::aes(u1, u2)) +
  ggplot2::geom_point(alpha = 0.25, size = 0.6, colour = "#2a78d6") +
  ggplot2::facet_wrap(~copula) +
  ggplot2::coord_equal() +
  ggplot2::theme_minimal()

Arbitrary margins: `r_mvdist()`

r_mvdist() couples any copula with any marginal distributions, returning a data frame. Margins can be names of distributions with a q* function, lists of parameters, or quantile functions — including categorical margins from margin_categorical() and empirical margins learned from data with margin_empirical().

sim <- r_mvdist(
  1000,
  margins = list(
    income = list(dist = "lnorm", meanlog = 10, sdlog = 0.6),
    age = list(dist = "unif", min = 18, max = 80),
    children = list(dist = "pois", lambda = 1.5),
    sector = margin_categorical(c("agri", "industry", "services"), c(0.1, 0.3, 0.6)),
    eruption = margin_empirical(faithful$eruptions)
  ),
  copula = "t", corr = make_corr(5, "exchangeable", 0.4), df = 4, seed = 1
)
str(sim)
#> 'data.frame':    1000 obs. of  5 variables:
#>  $ income  : num  8888 23444 12856 60544 27701 ...
#>  $ age     : num  77.7 63.1 25.5 70.3 54.6 ...
#>  $ children: num  0 1 3 3 1 1 2 0 0 1 ...
#>  $ sector  : Factor w/ 3 levels "agri","industry",..: 3 3 3 3 1 3 3 2 3 3 ...
#>  $ eruption: num  1.75 4.29 4.35 2.23 1.79 ...
round(association_matrix(sim), 2)
#>          income  age children sector eruption
#> income     1.00 0.36     0.40   0.35     0.38
#> age        0.36 1.00     0.40   0.38     0.41
#> children   0.40 0.40     1.00   0.39     0.31
#> sector     0.35 0.38     0.39   1.00     0.37
#> eruption   0.38 0.41     0.31   0.37     1.00

Such data sets are convenient test beds for synthesizers, because the true dependence is known:

syn <- synthesize(sim, seed = 2)
summary(compare_synthetic(sim, syn, seed = 1))
#>                          metric      value ideal
#> 1        mean marginal distance 0.02480000     0
#> 2         max marginal distance 0.03600000     0
#> 3 mean |association difference| 0.02876745     0
#> 4                    pMSE ratio 1.60863755    ~1
#> 5             discriminator AUC 0.48352100   0.5
#> 6               exact copy rate 0.00000000     0
#> 7                     DCR ratio 1.08575375  >= 1

Generating data from multivariate distributions

Correlation matrices

Classical multivariate distributions

Copulas

Arbitrary margins: r_mvdist()

References

Arbitrary margins: `r_mvdist()`