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.


  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?


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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.