bug: `Beta::cdf` loses accuracy above shapes of ~1e4 and returns values outside [0, 1] by 1e8

[FBruzzesi]

Same context and premise as in #432

The issue

Beta(s, s).cdf(0.5) is exactly 0.5 for every s, by the symmetry I_{1/2}(s, s) = 1/2, so this case needs no oracle. What I get:

s	Beta(s, s).cdf(0.5)	true	rel. err.
1e2	0.5000000000001004	0.5	2.0e-13
1e3	0.5000000000000433	0.5	8.7e-14
1e4	0.5000000000286743	0.5	5.7e-11
1e5	0.49999969504052444	0.5	6.1e-7
1e6	0.49121972699639316	0.5	1.8e-2
1e7	0.21285001452392216	0.5	5.7e-1
1e8	-1.1473179165955987	0.5	outside [0, 1]

The last row is negative, so cdf(0.5).ln() is NaN there. sf(0.5) returns the same number, so cdf + sf is not 1.

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() {
    for s in [1e2f64, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8] {
        let d = Beta::new(s, s).unwrap();
        println!("Beta({s:e}, {s:e}).cdf(0.5) = {:<22} ln = {:?}", d.cdf(0.5), d.cdf(0.5).ln());
    }
}

Beta(1e2, 1e2).cdf(0.5) = 0.5000000000001004     ln = -0.6931471805597446
Beta(1e3, 1e3).cdf(0.5) = 0.5000000000000433     ln = -0.6931471805598587
Beta(1e4, 1e4).cdf(0.5) = 0.5000000000286743     ln = -0.6931471805025967
Beta(1e5, 1e5).cdf(0.5) = 0.49999969504052444    ln = -0.6931477904790825
Beta(1e6, 1e6).cdf(0.5) = 0.49121972699639316    ln = -0.7108637421323526
Beta(1e7, 1e7).cdf(0.5) = 0.21285001452392216    ln = -1.547167518472497
Beta(1e8, 1e8).cdf(0.5) = -1.1473179165955987    ln = NaN

Expected

0.5 in every row.

Possibly related

I only tested Beta here. If Binomial reaches the incomplete beta for P(X > k), large n may land in the same range, but I have not checked that.

[RektPunk]

I think

statrs/src/function/beta.rs

Lines 191 to 232 in 02fdab6

	for m in 1..141 {
	let m = f64::from(m);
	let m2 = m * 2.0;
	let mut aa = m * (b - m) * x / ((qam + m2) * (a + m2));
	d = 1.0 + aa * d;
	if d.abs() < fpmin {
	d = fpmin;
	}
	c = 1.0 + aa / c;
	if c.abs() < fpmin {
	c = fpmin;
	}
	d = 1.0 / d;
	h = h * d * c;
	aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2));
	d = 1.0 + aa * d;
	if d.abs() < fpmin {
	d = fpmin;
	}
	c = 1.0 + aa / c;
	if c.abs() < fpmin {
	c = fpmin;
	}
	d = 1.0 / d;
	let del = d * c;
	h *= del;
	if (del - 1.0).abs() <= eps {
	return if symm_transform {
	Ok(1.0 - bt * h / a)
	} else {
	Ok(bt * h / a)
	};
	}
	}

is where the issue happens. It seems like the algorithm fails to converge for large a and b values. It silently returns an unconverged value.

[YeungOnion]

Resolved by #436 (commit included in PR #438)

[day01]

@YeungOnion it is not fully resolved yet, it has soem edge cases.

[YeungOnion]

Oh I closed the wrong thing.