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 two rounds.

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

Define the random variables

$Y_1$: number of outcomes "head" in Round 1,

$X_2$: number of outcomes "head" in Round 2,

$Y_2=Y_1+X_2$.

I learned that

• $Y_1\sim\text{Poisson-Binomial}(\{p_1,\ldots,\,p_{100}\})$,
• $X_2\sim\text{Poisson-Binomial}(\{q_1,\ldots,\,q_{100}\})$, where $q_j=(1-p_j)\cdot p_j$, for all $j\in\{1,\ldots,100\}$, and
• $Y_2\sim\text{Poisson-Binomial}(\{r_1,\ldots,\,r_{100}\})$, where $r_j=p_j+q_j$, for all $j\in\{1,\ldots,100\}$.

Problem: Obtain efficiently the joint distribution of the random vector $(Y_1,\,Y_2)$ whose range is $\{(y_1,\,y_2)\in\{0,\,1,\ldots,\,100\}^2:\,y_2\geq y_1\}$.

Note that $\mathbb{P}(Y_1=y_1,\,Y_2=y_2)=\mathbb{P}(X_2=y_2-y_1 \mid Y_1=y_1) \cdot \mathbb{P}(Y_1=y_1)$.

If $y_1=0$, then $\mathbb{P}(Y_1=0,\,Y_2=y_2)=\left(\displaystyle\prod_{j=1}^{100}(1-p_j)\right) \cdot \mathbb{P}(X_2=y_2\mid Y_1=0)$,

where the second factor can be computed efficiently with the ${\tt R}$-command ${\tt dpoibin}$, for all $y_2\in\{0,\ldots,\,100\}$.

If $y_1=100$, then $\mathbb{P}(Y_1=100,\,Y_2=y_2)=\left(\displaystyle\prod_{j=1}^{100}p_j\right)\cdot 1$.

Troubles to compute $\mathbb{P}(Y_1=y_1,\,Y_2=y_2)$ arise when $y_1\in\{1,\,2,\ldots,\,99\}$ and $y_2\in\{y_1,\ldots,\,100\}$.

Does anyone know how to compute efficiently $\mathbb{P}(Y_1=y_1,\,Y_2=y_2)$, for all $(y_1,\,y_2)$? Thanks a lot for your help and suggestions.

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.

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

Answer 2

There is no efficient theoretical way to do this, but it is a straightforward dynamic programming problem for a computer. Here is sample code for it in Python.

#! /usr/bin/env python3
def transition_distribution (prob_vector):
    transitions = {(0,0): 1.0}
    for p_i in prob_vector:
        new_transitions = {}
        for entry, prob in transitions.items():
            entry_adjustments = {
                entry: prob * (1.0-p_i) * (1.0-p_i),
                tuple([entry[0], entry[1]+1]): prob * (1.0-p_i) * p_i,
                tuple([entry[0]+1, entry[1]+1]): prob * p_i,
                }
            for new_entry, new_prob in entry_adjustments.items():
                new_transitions[new_entry] = new_transitions.get(new_entry, 0.0) + new_prob
        transitions = new_transitions
    return transitions

print(transition_distribution([0.5, 0.3, 0.2]))

Note that with 100 coins you'll be throwing around dictionaries with 5,500 keys. But that isn't too hard for a computer to do.