Efficiently computing poisson binomial sum

Question

Suppose there is a match between two teams where the first team to win a certain number of games wins the match. The match is handicapped, with Teams A needing to win $H_A$ games and Team B needing $H_B$ games. Without loss of generality, let $H_A \le H_B$. The odds of each game are different; let $p_i$ be the probability that Team A wins game $i$, where $1 \le i < H_A + H_B$. What is the probability that Team A wins the match (i.e. wins exactly $H_A$ games and loses fewer than $H_B$)?

I believe the chance of Team A winning is

$$\sum_{n=H_A}^{H_A+H_B-1} \Pr(X_{n-1}=H_A-1) \cdot p_n$$

where $X_n$ is a Poisson binomial distribution over $n$ trials and $\Pr(X_n=k)$ is the probability of $k$ successes. Explanation:

• If A wins, there must be at least $n = H_A$ games (A wins every game) and at most $H_A+H_B-1$ games (A wins $H_A$, B wins $H_B-1$).
• The winner always wins the last game, thus the $p_n$ term.
• In the previous $n-1$ games, A must have won $H_A-1$ times.

Questions:

1. Is the formula above correct?
2. How can I compute this efficiently? In my case, $H_A + H_B < 32$, and I need to compute this for a large number of sets of $\{p_i\}$. I saw something about using an FFT but haven't seen a simple explanation.

See also: How can I efficiently model the sum of Bernoulli random variables?

Answer 1

I believe my original formula is correct, but it can also be written as:

$$\Pr(\text{Team A wins}) = \sum_{n=H_A}^{H_A+H_B-1} \Pr(X_n \ge H_A)$$

This is actually very easy to implement. It is $O(n^2)$ and takes less than a second for $n \le 10000$.

import numpy as np

def poisson_binomial_pmf(p):
    """Returns the p.m.f. of the Poisson Binomial distribution.

    The Poisson Binomial distribution is the sum of `len(p)` independent yes/no
    trials, with trial `i` having success probability `p[i]`. (The special case
    of all probabilities being equal is the Binomial distribution.)

    Return value is a numpy array of length `len(p)+1`, where element `i` is the
    chance of `i` successes.

    >>> poisson_binomial_pmf([.25])
    array([ 0.75,  0.25])
    >>> poisson_binomial_pmf([.5, .25])
    array([ 0.375,  0.5  ,  0.125])
    >>> poisson_binomial_pmf([.5, .25, .125])
    array([ 0.328125,  0.484375,  0.171875,  0.015625])
    """
    n = len(p)
    pmf = np.zeros(n + 1, dtype=float)
    pmf[0] = 1.0
    for i in range(n):
        success = pmf[:i + 1] * p[i]
        pmf[:i + 2] *= (1 - p[i])  # failure
        pmf[1:i + 2] += success
    return pmf

def match_probability(game_probability, handicap_home, handicap_away):
    """Returns the probability of the home team winning."""
    max_games = handicap_home + handicap_away - 1
    game_chance = game_chance[:max_games]
    assert len(game_chance) == max_games
    return poisson_binomial_pmf(game_chance)[handicap_home:].sum()