# Simulating realizations of joint Bernoulli distribution

Let $$X$$ and $$Y$$ be Bernoulli random variables with success probability $$p$$ and $$q$$ respectively, i.e., \begin{align*} X = \begin{cases} 1 & \text{with probability p} \\ 0 & \text{with probability 1-p} \end{cases} \qquad,\qquad Y = \begin{cases} 1 & \text{with probability q} \\ 0 & \text{with probability 1-q} \end{cases} \end{align*} Let $$\mathbb{P}( X = 1 \text{ and } Y= 1 ) = \tfrac 13$$.

Suppose you have a function Bernoulli$$(r)$$ that simulates the outcomes of a Bernoulli random variable for any success probability $$r$$, i.e., it will generate a $$1$$ with probability $$r$$ and $$0$$ with probability $$1-r$$. Explain how you can use this function to generate a realization of $$X$$ and $$Y$$ above for $$p=\tfrac 12$$, $$q=\tfrac 23$$.

Since $$X$$ and $$Y$$ are independent for these specific probabilities (because $$Cov(X,Y) = 0$$ and uncorrelatedness implies independence for Bernoulli R.V.s), can we simply generate $$N$$ outcomes using Bernoulli($$p$$) and $$N$$ outcomes for Bernoulli($$q$$), and then just pair the i-th outcome of $$X$$ with the i-th outcome of $$Y$$?

Yes, when they're independent you can just simulate from each of the Bernoullis "individually" without regard for the other; which observation you pair them with can be chosen arbitrarily (as long as it's without regard to the values they take), so you can just pair up $$i$$th values happily, or indeed scramble them up any other way you like.
When they're not independent, you could simulate from the four joint probabilities $$p_{00}, p_{01}, p_{10}$$ and $$p_{11}$$ according to their specific values.