Bivariate Poisson-Binomial distribution

Question

Suppose you have $100$ coins whose probabilities of obtaining the outcome "head" are $p_1,\ldots,\,p_{100}$. These probabilities are not necessarily equal each other. Consider the following random experiment divided into rounds.

• Round 1: Throw simultaneously the $100$ coins and observe the number of heads.
• Round 2: Throw only those coins you obtained the outcome "tail" in Round 1, and observe the number of heads.
• Round 3: Throw only those coins you obtained the outcome "tail" in Round 2, and observe the number of heads.

... and so on. The experiment ends when no coins remain to throw. For all $k\in\{1,\,2,\,3,\ldots\}$, define the random variables

$X_k$: number of outcomes "head" in Round $k$, and

$Y_k=X_1+\cdots+X_k$.

I learned that

• $X_1\sim\text{Poisson-Binomial}(n=100,\,\{p_1,\ldots,\,p_{100}\})$ and that
• $Y_k\sim\text{Poisson-Binomial}(n=100,\,\{q_1,\ldots,\,q_{100}\})$, where $q_j=p_j+(1-p_j)\cdot p_j+\cdots+(1-p_j)^{k-1}\cdot p_j$, for all $j\in\{1,\ldots,\,100\}$.

I would like to obtain the joint probability distribution of the random vector $(Y_k,\,Y_{k+1})$.

For the case $k=1$, I noted that $\mathbb{P}(Y_1=y_1,\,Y_2=y_2)=\mathbb{P}(X_1=y_1,\,X_2=y_2-y_1)$. To compute this probability, I have to sum the probabilities of all ${100\choose y_1}\cdot{100-y_1\choose y_2-y_1}$ pairs $(X_1,\,X_2)$ that produce the event $\{X_1=y_1,\,X_2=y_2-y_1\}$. However, I believe this method is inefficient.

Question: Is there any efficient method or any approximation to obtain the probability distribution of the random vector $(Y_k,\,Y_{k+1})$, for all $k\in\{1,\,2,\,3,\ldots\}$?

Any suggestion will be very welcomed. Thanks a lot for your help.

Answer 1

I found a solution to my problem. This solution builds on the paper

Nelsen, R. B. (1987). Discrete bivariate distributions with given marginals and correlation. Communications in Statistics-Simulation and Computation, 16(1), 199-208.

Since $\mathbb{P}(Y_1=y_1,\,Y_2=y_2)=\mathbb{P}(Y_1=y_1,\,X_2=y_2-y_1)$, obtaining the joint probability distribution of the random vector $(Y_1,\,X_2)$ is enough.

Note that $\mathbb{C}\mbox{ov}(Y_1,\,X_2)=0.5\cdot(\mathbb{V}\mbox{ar}(Y_2)-\mathbb{V}\mbox{ar}(Y_1)-\mathbb{V}\mbox{ar}(X_2))$, with

• $\displaystyle\mathbb{V}\mbox{ar}(Y_1)=\sum_{j=1}^{100}(p_j\cdot(1-p_j))$,
• $\displaystyle\mathbb{V}\mbox{ar}(X_2)=\sum_{j=1}^{100}(q_j\cdot(1-q_j))$, and
• $\displaystyle\mathbb{V}\mbox{ar}(Y_2)=\sum_{j=1}^{100}(r_j\cdot(1-r_j))$.

Following notation of Nelsen (1987), one can

• obtain $f(x)$ and $f(y)$ efficiently with ${\tt dpoibin}$ in ${\tt R}$-software,
• obtain $F(x)$ and $F(y)$ efficiently with ${\tt ppoibin}$ in ${\tt R}$-software, and
• use $\rho=\frac{\mathbb{C}\text{ov}(Y_1,\,X_2)}{\sqrt{\mathbb{V}\mbox{ar}(Y_1)\cdot\mathbb{V}\mbox{ar}(X_2)}}$.

Thus, according to the sign for $\rho$, one can obtain the joint probability distribution of $(Y_1,\,X_2)$ with Nelsen's approach faster than with the approach I provided in my question.

Finally, one can use the same approach to obtain the joint probability distribution of $(Y_k,\,Y_{k+1})$, for all $k$.

If you have comments or corrections concerning this solution, please let me know. Thanks a lot.