Sampling from multivariate Bernoulli Suppose you have a vector p drawn from a multivariate Beta distribution (not a Dirichlet), such as the one described here ( How to construct a multivariate Beta distribution? ) with a Gaussian copula.
How would you make a draw from a multivariate Bernoulli distribution parametrised in terms of p?

*

*Draw from your multivariate distribution with Beta margins.


*Use the marginal Beta distributions as $p$ parameters for the marginal Bernoulli variables.
From this construction, you wind up with dependent Bernoulli variables.
library(copula)
set.seed(2022)

# Multivariate sample size
#
N <- 1000

# Parameter of the Gaussian copula
#
rho <- 0.95 

# Define the Gaussian copula
#
cop <- copula::normalCopula(rho) 

# Define the multivariate distribution with the 
# Gaussian copula "cop" and Beta(1, 1) margins
#
mvbeta <- copula::mvdc(
  cop, c("beta", "beta"), 
  list(
    list(shape1 = 1, shape2 = 1), 
    list(shape1 = 1, shape2 = 1)
    )
  )

# Draw N points from the multivariate distribution
#
p <- copula::rMvdc(N, mvbeta)

# Extract the margins
#
p1 <- p[, 1]
p2 <- p[, 2]

# Define the two Bernoulli margins that, by this construction, are dependent
#
bernoulli_1 <- rbinom(N, 1, p1)
bernoulli_2 <- rbinom(N, 1, p2)

# Test the correlation between the Bernoulli varables.
# Significant correlation is evidence of dependence derived 
# from the Gaussian copula of "mvbeta" 
#
cor.test(bernoulli_1, bernoulli_2)

The Pearson correlation considered by cor.test isn't the only way to test that bernoulli_1 and bernoulli_2 are dependent, but it is a quick way.
