fix: Stabilize beta CDF and extreme quantiles

Cargo.toml

@@ -3,14 +3,22 @@ name = "statrs"
 version = "0.19.1"
 authors = ["Michael Ma"]
 description = "Statistical computing library for Rust"
-license = "MIT"
+license = "MIT AND BSL-1.0"
 keywords = ["probability", "statistics", "stats", "distribution", "math"]
 categories = ["science"]
 homepage = "https://github.com/statrs-dev/statrs"
 repository = "https://github.com/statrs-dev/statrs"
 edition = "2024"
 
-include = ["CHANGELOG.md", "LICENSE.md", "src/", "tests/"]
+include = [
+    "CHANGELOG.md",
+    "LICENSE.md",
+    "LICENSE-BOOST.md",
+    "THIRD_PARTY_NOTICES.md",
+    "src/",
+    "tests/",
+    "benches/",
+]
 
 # When changing MSRV: Also update the README
 rust-version = "1.89.0"
@@ -29,6 +37,11 @@ name = "density"
 harness = false
 required-features = ["rand", "std", "kde"]
 
+[[bench]]
+name = "beta"
+harness = false
+required-features = ["std"]
+
 [features]
 default = ["std", "nalgebra", "rand"]
 std = ["approx/std", "num-traits/std", "nalgebra?/std", "rand?/std"]

LICENSE-BOOST.md

@@ -0,0 +1,23 @@
+Boost Software License - Version 1.0 - August 17th, 2003
+
+Permission is hereby granted, free of charge, to any person or organization
+obtaining a copy of the software and accompanying documentation covered by
+this license (the "Software") to use, reproduce, display, distribute,
+execute, and transmit the Software, and to prepare derivative works of the
+Software, and to permit third-parties to whom the Software is furnished to
+do so, all subject to the following:
+
+The copyright notices in the Software and this entire statement, including
+the above license grant, this restriction and the following disclaimer,
+must be included in all copies of the Software, in whole or in part, and
+all derivative works of the Software, unless such copies or derivative
+works are solely in the form of machine-executable object code generated by
+a source language processor.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
+SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
+FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
+ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+DEALINGS IN THE SOFTWARE.

README.md

@@ -1,7 +1,7 @@
 # statrs
 
 ![tests][actions-test-badge]
-[![MIT licensed][license-badge]](./LICENSE.md)
+[![MIT and BSL-1.0 licensed][license-badge]](./THIRD_PARTY_NOTICES.md)
 [![Crate][crates-badge]][crates-url]
 [![docs.rs][docsrs-badge]][docs-url]
 [![codecov-statrs][codecov-badge]][codecov-url]
@@ -10,7 +10,7 @@
 [actions-test-badge]: https://github.com/statrs-dev/statrs/actions/workflows/test.yml/badge.svg
 [crates-badge]: https://img.shields.io/crates/v/statrs.svg
 [crates-url]: https://crates.io/crates/statrs
-[license-badge]: https://img.shields.io/badge/license-MIT-blue.svg
+[license-badge]: https://img.shields.io/badge/license-MIT%20AND%20BSL--1.0-blue.svg
 [docsrs-badge]: https://img.shields.io/docsrs/statrs
 [docs-url]: https://docs.rs/statrs/*/statrs
 [codecov-badge]: https://codecov.io/gh/statrs-dev/statrs/graph/badge.svg?token=XtMSMYXvIf

THIRD_PARTY_NOTICES.md

@@ -0,0 +1,24 @@
+# Third-party notices
+
+## Boost.Math
+
+Portions of `src/function/beta/{bgrat,forward,fraction,log_forward,recurrence,series}.rs` are adapted from Boost.Math 1.90.0, `include/boost/math/special_functions/beta.hpp`.
+
+Portions of `src/function/beta/normal_tail.rs` are adapted from Boost.Math 1.90.0, `include/boost/math/special_functions/erf.hpp`.
+
+Copyright John Maddock 2006.
+Copyright Matt Borland 2024.
+
+The Boost-derived portions are licensed under the Boost Software License 1.0; see `LICENSE-BOOST.md`. Statrs modifications are licensed under MIT, so these files are subject to both licenses.
+
+Source: https://github.com/boostorg/math/blob/e0fcd19f7227d81391770ea46015acc3c80af810/include/boost/math/special_functions/beta.hpp
+
+Source: https://github.com/boostorg/math/blob/e0fcd19f7227d81391770ea46015acc3c80af810/include/boost/math/special_functions/erf.hpp
+
+## special 0.8.1
+
+The initial inverse-beta estimate in `src/function/beta/inverse/initial.rs` is adapted from `special` 0.8.1 under its MIT license option.
+
+Copyright 2014–2019 The special Developers.
+
+Source: https://github.com/stainless-steel/special/blob/c64902430bd50e8c8225c7c8b410334ffedf2f15/src/beta.rs

benches/beta.rs

@@ -0,0 +1,70 @@
+use criterion::{BenchmarkId, Criterion, criterion_group, criterion_main};
+use statrs::function::beta::{beta_reg, inv_beta_reg};
+use std::hint::black_box;
+
+fn bench_beta_reg(c: &mut Criterion) {
+    let mut group = c.benchmark_group("beta_reg");
+    for (name, a, b, x) in [
+        ("typical", 2.0, 5.0, 0.3),
+        (
+            "large_symmetric_adjacent",
+            1e8,
+            1e8,
+            f64::from_bits(0.5_f64.to_bits() + 1),
+        ),
+        (
+            "moderate_fraction",
+            25.32628846940565,
+            3.1028101710805442,
+            0.9276950604606229,
+        ),
+    ] {
+        group.bench_with_input(
+            BenchmarkId::new("cdf", name),
+            &(a, b, x),
+            |bencher, input| {
+                bencher
+                    .iter(|| beta_reg(black_box(input.0), black_box(input.1), black_box(input.2)));
+            },
+        );
+    }
+    group.finish();
+}
+
+fn bench_inv_beta_reg(c: &mut Criterion) {
+    let mut group = c.benchmark_group("inv_beta_reg");
+    for (name, a, b, probability) in [
+        ("typical", 2.0, 5.0, 0.3),
+        ("nontermination_regression", 200.0, 2.0, 1e-60),
+        ("panic_regression", 200.0, 2.0, 1e-165),
+        ("tiny_quantile", 0.1, 500.0, 1e-30),
+        (
+            "subnormal_shape_two",
+            0.5,
+            2.0,
+            f64::from_bits(0x1e72_f942_2c23_c47c),
+        ),
+        ("subnormal_unit_shape", 1.0, 10.0, f64::from_bits(5)),
+        ("shape_two_large_median", 2.0, 1e308, 0.5),
+        (
+            "shape_two_zero_cell",
+            2.0,
+            1e200,
+            f64::from_bits(0x0c8b_4ec7_f919_73ff),
+        ),
+    ] {
+        group.bench_with_input(
+            BenchmarkId::new("quantile", name),
+            &(a, b, probability),
+            |bencher, input| {
+                bencher.iter(|| {
+                    inv_beta_reg(black_box(input.0), black_box(input.1), black_box(input.2))
+                });
+            },
+        );
+    }
+    group.finish();
+}
+
+criterion_group!(benches, bench_beta_reg, bench_inv_beta_reg);
+criterion_main!(benches);

src/distribution/beta.rs

@@ -129,6 +129,10 @@ impl ContinuousCDF<f64, f64> for Beta {
     /// Calculates the cumulative distribution function for the beta
     /// distribution at `x`.
     ///
+    /// # Panics
+    ///
+    /// If the numerical method does not converge.
+    ///
     /// # Formula
     ///
     /// ```text
@@ -151,6 +155,10 @@ impl ContinuousCDF<f64, f64> for Beta {
 
     /// Calculates the survival function for the beta distribution at `x`.
     ///
+    /// # Panics
+    ///
+    /// If the numerical method does not converge.
+    ///
     /// # Formula
     ///
     /// ```text
@@ -167,12 +175,9 @@ impl ContinuousCDF<f64, f64> for Beta {
         } else if self.shape_a == 1.0 && self.shape_b == 1.0 {
             1. - x
         } else if x < (self.shape_a + 1.0) / (self.shape_a + self.shape_b + 2.0) {
-            // Below the continued fraction split point of `beta_reg`,
-            // `beta_reg(b, a, 1 - x)` reduces to `1 - beta_reg(a, b, x)`;
-            // computing the complement here instead avoids `1.0 - x`
-            // rounding to 1.0 for tiny x (< ~1.1e-16), which would lose
-            // the lower tail entirely. See #432
-            1.0 - beta::beta_reg(self.shape_a, self.shape_b, x)
+            beta::checked_ln_beta_reg_complement(self.shape_a, self.shape_b, x)
+                .unwrap()
+                .exp()
         } else {
             beta::beta_reg(self.shape_b, self.shape_a, 1.0 - x)
         }
@@ -183,7 +188,7 @@ impl ContinuousCDF<f64, f64> for Beta {
     ///
     /// # Panics
     ///
-    /// If x is not in `[0, 1]`.
+    /// If x is not in `[0, 1]` or the numerical method does not converge.
     ///
     /// # Formula
     ///
@@ -208,6 +213,10 @@ impl ContinuousCDF<f64, f64> for Beta {
     ///
     /// If x is not in `[0, 1]`.
     ///
+    /// # Panics
+    ///
+    /// If the numerical method does not converge.
+    ///
     /// # Formula
     ///
     /// ```text
@@ -651,6 +660,18 @@ mod tests {
         }
     }
 
+    #[test]
+    fn test_cdf_large_symmetric_shapes() {
+        for shape in [1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8] {
+            let dist = Beta::new(shape, shape).unwrap();
+            let cdf = dist.cdf(0.5);
+            let sf = dist.sf(0.5);
+            assert_eq!(cdf, 0.5);
+            assert_eq!(sf, 0.5);
+            assert_eq!(cdf + sf, 1.0);
+        }
+    }
+
     #[test]
     fn test_sf() {
         let sf = |arg: f64| move |x: Beta| x.sf(arg);
@@ -684,6 +705,17 @@ mod tests {
         }
     }
 
+    #[test]
+    fn test_sf_tiny_shape_preserves_representable_tail() {
+        let distribution = Beta::new(
+            f64::from_bits(0x00000000000007e8),
+            f64::from_bits(0x4040000000000000),
+        )
+        .unwrap();
+        let x = f64::from_bits(0x01556e1fc2f8f359);
+        assert_eq!(distribution.sf(x).to_bits(), 0x0000000000155101);
+    }
+
     #[test]
     fn test_inverse_cdf() {
         // let inverse_cdf = |arg: f64| move |x: Beta| x.inverse_cdf(arg);
@@ -705,6 +737,18 @@ mod tests {
         }
     }
 
+    #[test]
+    fn test_inverse_cdf_extreme_lower_tail() {
+        let dist = Beta::new(200.0, 2.0).unwrap();
+        let actual = dist.inverse_cdf(1e-170);
+        let expected = 0.13765877485659653;
+        let relative_error = ((actual - expected) / expected).abs();
+        assert!(
+            relative_error <= 5e-13,
+            "actual {actual}, expected {expected}"
+        );
+    }
+
     #[test]
     fn test_cdf_input_lt_0() {
         let cdf = |arg: f64| move |x: Beta| x.cdf(arg);

src/distribution/binomial/mod.rs

@@ -514,7 +514,7 @@ mod tests {
         test_absolute(0.3, 3, 0.657, 1e-14, sf(0));
         test_absolute(0.3, 3, 0.216, 1e-15, sf(1));
         test_exact(0.3, 3, 0.0, sf(3));
-        test_absolute(0.3, 10, 0.9717524751000001, 1e-16, sf(0));
+        test_absolute(0.3, 10, 0.9717524751, 1e-16, sf(0));
         test_absolute(0.3, 10, 0.850691654100002, 1e-14, sf(1));
         test_exact(0.3, 10, 0.0, sf(10));
         test_exact(1.0, 1, 1.0, sf(0));