Sample three Bernoulli variables given the 2 by 2 by 2 table I want to sample three Bernoulli variables given their 2 by 2 by 2 table:

where the last column is the probabilities with sum being equal to 1.
What I am thinking now is to use a multivariate normal distribution, then use cut-points to create Bernoulli variables.
Do you know any better way, directly sample from Bernoulli distributions?

Suppose that you want to sample a Bernoulli trial for $X$, then you have to know the probability of success $\mathbb{P}(X=1)$
You can calculate this probability by summing over the other two variables.
$$\mathbb{P}(X=1)=\sum_{Y_{0},Y_{1}=0,1}\mathbb{P}(X=1,Y_{0},Y_{1})$$
$$=\mathbb{P}(X=1,Y_{0}=0,Y_{1}=0)+\mathbb{P}(X=1],Y_{0}=1,Y_{1}=0) +$$
$$ \mathbb{P}(X=1,Y_{0}=0,Y_{1}=1) + \mathbb{P}(X=1,Y_{0}=1,Y_{1}=1)$$
$$= \pi_{100} + \pi_{110} + \pi_{101} +\pi_{111}$$
Hence, $X\sim Bern(\pi_{100} + \pi_{110} + \pi_{101} +\pi_{111})$ similarly for the rest.
If you want to take into account the association between the random variables, then I suppose you can use a Multinomial distribution.
Where for example the event $(0,0,0)$ will have probability $\pi_{000}$ similarly you can think for the rest. So,  (not 100% sure but seems logical to me) for sampling Bernoulli trials you can think in the following way:

*

*You conduct a single Multinomial sampling with probabilities $(\pi_{000},\pi_{001},\pi_{010},\pi_{100},\pi_{011},\pi_{101},\pi_{110},\pi_{111})$


*You identify the result that you produced, for example if it is $(1,0,1)$ you can regard it as $X=1$, $Y_{0}=0$ and $Y_{1}=1$.
In R you could do something like:
probs = rdirichlet(1,rep(1,8))
Mult = matrix(0,1000,8)
for(i in 1:1000){
  Mult[i,] = rmultinom(1,1,prob)
}
X = sum(apply(Mult,2,sum)[c(4,6,7,8)]) # how many 1 draws
Y_0 = sum(apply(Mult,2,sum)[c(3,5,7,8)])
Y_1 = sum(apply(Mult,2,sum)[c(2,5,6,8)])

