Use h function

docs/index.md

@@ -45,3 +45,18 @@ g(k, q, m) = \begin{cases}
   \frac{1}{k!} - \sum_{i=0}^{k-1} \frac{q^{(m+i)(k-i)}}{(k-i)!} g(i, q, m) & k > 0
 \end{cases}
 ```
+
+Further define $h(k, q, m) = k! \times g(k, q, m)$ such that:
+
+```math
+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
+
+```math
+f(k; n, m, p) = \binom{n}{k} q^{(n-k)(m+k)} h(k, q, m)
+```

src/reedfrost/__init__.py

@@ -4,14 +4,13 @@
 import scipy.stats
 from numpy.typing import NDArray
 from scipy.optimize import brentq
-from scipy.special import gamma, rgamma
+from scipy.special import comb
 
 
 @functools.cache
-def _gontcharoff1(k: int, q: float, m: int) -> float:
+def _kgontcharoff1(k: int, q: float, m: int) -> float:
     """
-    Gontcharoff polynomials, specific to the Lefevre & Picard
-    formulations of Reed-Frost final outbreak size pmf calculations
+    Gontcharoff polynomials times k!, for a single value of k
 
     See Lefevre & Picard 1990 (doi:10.2307/1427595) equation 2.1
 
@@ -26,20 +25,20 @@ def _gontcharoff1(k: int, q: float, m: int) -> float:
     if k == 0:
         return 1.0
     else:
-        return rgamma(k + 1) - sum(
+        return 1.0 - sum(
             [
-                q ** ((m + i) * (k - i)) * rgamma(k - i + 1) * _gontcharoff1(i, q, m)
+                comb(k, i) * q ** ((m + i) * (k - i)) * _kgontcharoff1(i, q, m)
                 for i in range(0, k)
             ]
         )
 
 
-def _gontcharoff(
+def _kgontcharoff(
     k: int | NDArray[np.int64], q: float, m: int
 ) -> float | NDArray[np.float64]:
     """
     Gontcharoff polynomials, specific to the Lefevre & Picard
-    formulations of Reed-Frost final outbreak size pmf calculations
+    formulation, times k!
 
     See Lefevre & Picard 1990 (doi:10.2307/1427595) equation 2.1
 
@@ -52,9 +51,9 @@ def _gontcharoff(
         float, or float array: value of the polynomial
     """
     if isinstance(k, int):
-        return _gontcharoff1(k=k, q=q, m=m)
+        return _kgontcharoff1(k=k, q=q, m=m)
     else:
-        return np.array([_gontcharoff1(k=kk, q=q, m=m) for kk in k])
+        return np.array([_kgontcharoff1(k=kk, q=q, m=m) for kk in k])
 
 
 def pmf(
@@ -76,12 +75,7 @@ def pmf(
     """
     q = 1.0 - p
 
-    return (
-        gamma(n + 1)
-        * rgamma(n - k + 1)
-        * q ** ((n - k) * (m + k))
-        * _gontcharoff(k, q, m)
-    )
+    return comb(n, k) * q ** ((n - k) * (m + k)) * _kgontcharoff(k, q, m)
 
 
 def _theta_fun(w: float, lambda_: float) -> float: