Simulating realizations of joint Bernoulli distribution

Question

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$?

Answer 1

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.