Numerical stability of beta model

tskit-dev · Feb 8, 2024 · 7ee1619 · 7ee1619
1 parent 804e036
commit 7ee1619
Show file tree

Hide file tree

Showing 2 changed files with 8 additions and 2 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -1,5 +1,11 @@
 # Changelog
 
+## [1.3.1] - 2024-XX-XX
+
+**Bug fixes**:
+
+- Change tolerance of polynomial approximation in Beta-coalescent acceptance probabilities ({issue}`2256`, {pr}`2257`, {user}`JereKoskela`)
+
 ## [1.3.0] - 2022-12-13
 
 **New features**

diff --git a/lib/msprime.c b/lib/msprime.c
@@ -7056,7 +7056,7 @@ msp_beta_common_ancestor_event(msp_t *self, population_id_t pop_id, label_id_t l
     beta_x = ran_inc_beta(self->rng, 2.0 - alpha, alpha, truncation_point);
 
     /* We calculate the probability of accepting the event */
-    if (beta_x > 1e-9) {
+    if (beta_x > 1e-5) {
         u = (n - 1) * log(1 - beta_x) + log(1 + (n - 1) * beta_x);
         u = exp(log(1 - exp(u)) - 2 * log(beta_x) - gsl_sf_lnchoose(n, 2));
     } else {
@@ -7065,7 +7065,7 @@ msp_beta_common_ancestor_event(msp_t *self, population_id_t pop_id, label_id_t l
         u = 0;
         for (j = 2; j <= n; j += 2) {
             increment = (j - 1) * exp(gsl_sf_lnchoose(n, j) + (j - 2) * log(beta_x));
-            if (increment / u < 1e-12) {
+            if (increment / (u + increment) < 1e-12) {
                 /* We truncate the expansion adaptively once the increment
                  * becomes negligible. */
                 break;