Add a couple tests where 2F1 is singular

Author: nspope

tests/test_hypergeo.py

@@ -240,3 +240,28 @@ def test_2f1_grad(self, muts, hyp2f1_func, pars):
         offset = self._2f1_validate(*pars)
         check = self._2f1_grad_validate(*pars, offset=offset)
         assert np.allclose(grad, check)
+
+
+@pytest.mark.parametrize(
+    "pars",
+    [
+        # taken from examples in issues tsdate/286, tsdate/289
+        [1.104, 0.0001125, 118.1396, 0.009052, 1.0, 0.001404],
+        [2.7481, 0.001221, 344.94083, 0.02329, 3.0, 0.00026624],
+    ],
+)
+class TestSingular2F1:
+    """
+    Test detection of cases where 2F1 is close to singular and DLMF 15.8.3
+    suffers from catastrophic cancellation: in these cases, use DLMF 15.8.1
+    even though it takes much longer to converge.
+    """
+
+    def test_dlmf1583_throws_exception(self, pars):
+        with pytest.raises(Exception, match="is singular"):
+            hypergeo._hyp2f1_dlmf1583(*pars)
+
+    def test_exception_uses_dlmf1581(self, pars):
+        v1, *_ = hypergeo._hyp2f1(*pars)
+        v2, *_ = hypergeo._hyp2f1_dlmf1581(*pars)
+        assert np.isclose(v1, v2)

tsdate/approx.py

@@ -121,9 +121,6 @@ def sufficient_statistics(a_i, b_i, a_j, b_j, y_ij, mu_ij):
         a_i, b_i, a_j, b_j, y_ij, mu_ij
     )
 
-    if sign_f <= 0:
-        raise hypergeo.Invalid2F1("Singular hypergeometric function")
-
     logconst = (
         log_f + hypergeo._betaln(y_ij + 1, b) + hypergeo._gammaln(a) - a * np.log(t)
     )
@@ -139,6 +136,10 @@ def sufficient_statistics(a_i, b_i, a_j, b_j, y_ij, mu_ij):
         - hypergeo._digamma(c)
     )
 
+    # check that Jensen's inequality holds
+    assert np.log(t_i) > ln_t_i
+    assert np.log(t_j) > ln_t_j
+
     return logconst, t_i, ln_t_i, t_j, ln_t_j
 
 

tsdate/hypergeo.py

@@ -388,7 +388,7 @@ def _hyp2f1_dlmf1583(a_i, b_i, a_j, b_j, y, mu):
         )
         or sign <= 0
     ):
-        raise Invalid2F1("Hypergeometric series did not converge")
+        raise Invalid2F1("Hypergeometric series is singular")
 
     return val, sign, da_i, db_i, da_j, db_j