Implementation

The probability mass function for the final size distribution of Reed-Frost outbreaks is derived in Lefevre & Picard (1990) (equation 3.10):

\[ f(k; n, m, p) = n_{(k)} q^{(n-k)(m+k)} G_k(1 | q^{m+i}, i \in \mathbb{N}) \]

where:

• \(n_{(k)} = n!/(n-k)!\) is the falling factorial,
• \(k\) is the number of infections beyond the initial infections,
• \(n\) is the initial number susceptible (called \(S_0\) above),
• \(m\) is the initial number infected (called \(I_0\) above),
• \(p\) is the probability of effective contact (i.e., that any infected infects any susceptible in a time step),
• \(G_k\) are the Gontcharoff polynomials.

Gontcharoff polynomials

In general, the Gontcharoff polynomials are (equation 2.1):

\[ G_k(x | U) = \begin{cases} 1 & k = 0 \\ \frac{x^k}{k!} - \sum_{i=0}^{k-1} \frac{u_i^{k-i}}{(k-i)!} G_i(x | U) & k > 0 \end{cases} \]

For convenience, define \(g(k, q, m) \equiv G_k(1 | q^{m+i}, i \in \mathbb{N})\), a version of the Gontcharoff polynomials that are specific to this application, so that:

\[ g(k, q, m) = \begin{cases} 1 & k = 0 \\ \frac{1}{k!} - \sum_{i=0}^{k-1} \frac{q^{(m+i)(k-i)}}{(k-i)!} g(i, q, m) & k > 0 \end{cases} \]

Preventing overflow

In practice, the formulation for \(f\) and \(g\) above lead to numerical overflow. To avoid this, define \(h(k, q, m) = k! \times g(k, q, m)\) such that:

\[ h(k, q, m) = \begin{cases} 1 & k = 0 \\ 1 - \sum_{i=0}^{k-1} \binom{k}{i} q^{(m+i)(k-i)} h(i, q, m) & k > 0 \end{cases} \]

and

\[ f(k; n, m, p) = \binom{n}{k} q^{(n-k)(m+k)} h(k, q, m) \]

The use of binomial coefficients, rather than pure factorials, increases the range of values \(k\) and \(n\) before causing overflow.

Caveats

Large \(n\) can cause overflow, but in practice, numerical instability is a more pressing problem that leads to values of \(f\) greater than unity.