How to generate correlated Bernoulli variables? Let $X_1$ and $X_2$ be two Bernoulli random variables, with $P(X_1=1)=p_1$ and $P(X_2=1)=p_2$. The following discussion showed how to generate a pair of correlated Bernoulli variables with correlation coefficient $\rho$. I want to do the same but for 3 variables $X_1$, $X_2$, $X_3$ with respective probability $p_1$, $p_2$ and $p_3$ and identical correlation coefficient $\rho>0$ between $X_1$ and $X_2$, $X_1$ and $X_3$ and $X_2$ and $X_3$. Intuitively, it seems that this has multiple solutions depending on the higher order association between $X_1$, $X_2$, $X_3$, which I would then fix to 0. Is there a way to generate such distributions of Bernoulli variables ? By "generate", I don't mean produce with the help of an algorithm that would converge towards the desired distribution, but I mean more how we can construct it, i.e. calculate the probability of each combinations of the three Bernoulli variables.  If yes, can it be generalized to more than 3 variables (assuming null association between more than two variables)?
 A: To construct a measurable space let $\Omega=\left\{ 0,1\right\} ^{3}$
be equipped with its powerset as $\sigma$-algebra. 
It is handsome to think of an outcome $\omega\in\Omega$ as a function
$\left\{ 1,2,3\right\} \to\left\{ 0,1\right\} $.
For $i=1,2,3$ let $X_{i}:\Omega\to\mathbb{R}$ be prescribed by $\omega\mapsto\omega\left(i\right)$.
The outcome space has $2^{3}=8$ elements so that the $\sigma$-algebra contains $2^8$ events and we must find a probability measure $P$ on it that satisfies:


*

*$p_{1}=P\left(X_{1}^{-1}\left(\left\{ 1\right\} \right)\right)$ 

*$p_{2}=P\left(X_{2}^{-1}\left(\left\{ 1\right\} \right)\right)$ 

*$p_{3}=P\left(X_{3}^{-1}\left(\left\{ 1\right\} \right)\right)$ 

*$p_{1}p_{2}p_{3}=P\left(X_{1}^{-1}\left(\left\{ 1\right\} \right)\cap X_{2}^{-1}\left(\left\{ 1\right\} \right)\cap X_{3}^{-1}\left(\left\{ 1\right\} \right)\right)$

*$\rho\sqrt{p_{1}\left(1-p_{1}\right)p_{2}\left(1-p_{2}\right)}+p_{1}p_{2}=P\left(X_{1}^{-1}\left(\left\{ 1\right\} \right)\cap X_{2}^{-1}\left(\left\{ 1\right\} \right)\right)$

*$\rho\sqrt{p_{1}\left(1-p_{1}\right)p_{3}\left(1-p_{3}\right)}+p_{1}p_{3}=P\left(X_{1}^{-1}\left(\left\{ 1\right\} \right)\cap X_{3}^{-1}\left(\left\{ 1\right\} \right)\right)$

*$\rho\sqrt{p_{2}\left(1-p_{2}\right)p_{3}\left(1-p_{3}\right)}+p_{2}p_{3}=P\left(X_{2}^{-1}\left(\left\{ 1\right\} \right)\cap X_{3}^{-1}\left(\left\{ 1\right\} \right)\right)$
where the last $3$ equalities are based on the equality: $$\rho\left(X_{i},X_{j}\right)\sqrt{p_{i}\left(1-p_{i}\right)p_{j}\left(1-p_{j}\right)}=\rho\left(X_{i},X_{j}\right)\sigma_{X_{i}}\sigma_{X_{j}}=\mathsf{Cov}\left(X_{i},X_{j}\right)=\mathbb{E}X_{i}X_{j}-\mathbb{E}X_{i}\mathbb{E}X_{j}=$$$$\mathbb{E}X_{i}X_{j}-p_{i}p_{j}$$
There are $8$ disjoint sets of shape $X_{1}^{-1}\left(\left\{ x\right\} \right)\cap X_{2}^{-1}\left(\left\{ y\right\} \right)\cap X_{3}^{-1}\left(\left\{ z\right\} \right)$ involved (where $x,y,z\in\{0,1\}$) that cover the outcome space
and each of them has a probability. Giving these events a probability comes actually to the same as determining the probability measure. 
Denoting these probabilities as $a,u,v,w,r,s,t,z$ we meet the following
conditions:


*

*$1=a+u+v+w+r+s+t+z$

*$p_{1}=a+u+v+s=$

*$p_{2}=a+u+w+r$

*$p_{3}=a+v+w+t$

*$p_{1}p_{2}p_{3}=a$

*$\rho\sqrt{p_{1}\left(1-p_{1}\right)p_{2}\left(1-p_{2}\right)}+p_{1}p_{2}=a+u$

*$\rho\sqrt{p_{1}\left(1-p_{1}\right)p_{3}\left(1-p_{3}\right)}+p_{1}p_{3}=a+v$

*$\rho\sqrt{p_{2}\left(1-p_{2}\right)p_{3}\left(1-p_{3}\right)}+p_{2}p_{3}=a+w$
Here e.g. $a$ stands for the probability of $X_{1}^{-1}\left(\left\{ 1\right\} \right)\cap X_{2}^{-1}\left(\left\{ 1\right\} \right)\cap X_{3}^{-1}\left(\left\{ 1\right\} \right)$
and $u$ for the probability of $X_{1}^{-1}\left(\left\{ 1\right\} \right)\cap X_{2}^{-1}\left(\left\{ 1\right\} \right)\cap X_{3}^{-1}\left(\left\{ 0\right\} \right)$
and $s$ for the probability of $X_{1}^{-1}\left(\left\{ 1\right\} \right)\cap X_{2}^{-1}\left(\left\{ 0\right\} \right)\cap X_{3}^{-1}\left(\left\{ 0\right\} \right)$.
So we meet $8$ equalities on $8$ unknown variables. This however is not a guarantee that there is a proper solution because there are also constraints on base of inequalities. E.g. the probabilities cannot be negative.
