Efficiently computing poisson binomial sum

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.

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 = 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()