bug: `Beta::sf` returns exactly 1.0 for tiny x when `shape_a < 1`

[FBruzzesi]

Premise

First and foremost, thanks for all the works in statrs! It's really nice to work with it!

Recently, I have been using it to build a polars expression plugin, which I hope to open source soon enough. One last bit before doing so was to stress test tail accuracy.

Important

Mandatory disclaimer: this is my first experience writing rust code

The issue

Beta::sf(x) forms the complement before the call:

beta::beta_reg(self.shape_b, self.shape_a, 1.0 - x)

For small values of x (< 1.1e-16) the argument 1.0 - x rounds to exactly 1.0, so sf returns 1.0 regardless of the shapes.

With shape_a < 1 the lower tail carries real mass and the answer is gone: Beta(0.05, 0.05).sf(1e-16) returns 1.0 where the value is 0.92045095674306394, a relative error of 8.6e-2.

1.0 - cdf(x) on the same object is accurate to 16 digits, because on that side of the split beta_reg never sees a rounded argument. Only the routing is wrong.

Environment: statrs 0.19.0 (default-features = false, features = ["std"]), rustc 1.97.1, aarch64-apple-darwin, --release.

Reproducer

use statrs::distribution::{Beta, ContinuousCDF};

fn main() {
    let d = Beta::new(0.05, 0.95).unwrap();
    for x in [1e-16, 1e-18, 1e-20, 1e-30, 1e-100] {
        let y = 1.0 - x;
        println!("x = {x:>8e}   (1-x) = {y:>20e}  sf = {:.17}   1 - cdf = {:.17}", d.sf(x), 1.0 - d.cdf(x));
    }
}

x =    1e-16   (1-x) = 9.999999999999999e-1  sf = 1.00000000000000000   1 - cdf = 0.84216163834910041
x =    1e-18   (1-x) =                  1e0  sf = 1.00000000000000000   1 - cdf = 0.87462453281806818
x =    1e-20   (1-x) =                  1e0  sf = 1.00000000000000000   1 - cdf = 0.90041072647564402
x =    1e-30   (1-x) =                  1e0  sf = 1.00000000000000000   1 - cdf = 0.96850710651415306
x =   1e-100   (1-x) =                  1e0  sf = 1.00000000000000000   1 - cdf = 0.99999004107264755```

Expected

Reference from mpmath at mp.dps = 60 (same as scipy.stats.beta(0.05, 0.05).sf(values)):

x	statrs sf	true	statrs 1 - cdf	rel. err. of sf
1e-16	1.0	0.92045095674306394	0.92045095674306399	8.6e-2
1e-18	1.0	0.93681194889571247	0.93681194889571251	6.8e-2
1e-20	1.0	0.94980794691066360	0.94980794691066361	5.3e-2
1e-30	1.0	0.98412787917976062	0.98412787917976063	1.6e-2
1e-100	1.0	0.99999498079469107	0.99999498079469107	5.0e-6

Script for the reference column

from mpmath import mp, mpf, betainc
mp.dps = 60
for x in ["1e-16", "1e-18", "1e-20", "1e-30", "1e-100"]:
    v = 1 - betainc(mpf("0.05"), mpf("0.05"), 0, mpf(x), regularized=True)
    print(x, mp.nstr(v, 17))

Suggested fix (Claude assisted, due to above disclaimer)

--- a/src/distribution/beta.rs
+++ b/src/distribution/beta.rs
@@ impl ContinuousCDF<f64, f64> for Beta @@ fn sf(&self, x: f64) -> f64
         } else if prec::ulps_eq!(self.shape_a, 1.0) && prec::ulps_eq!(self.shape_b, 1.0) {
             1. - x
+        } else if x < (self.shape_a + 1.0) / (self.shape_a + self.shape_b + 2.0) {
+            1.0 - beta::beta_reg(self.shape_a, self.shape_b, x)
         } else {
             beta::beta_reg(self.shape_b, self.shape_a, 1.0 - x)
         }

---

FYI I would be interested in trying opening a PR

[day01]

fixed in #438