fix: Stabilize beta CDF and extreme quantiles #447
There are no files selected for viewing
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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. |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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 |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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), | ||
|
Comment on lines
+42
to
+45
There was a problem hiding this comment. 📐 Maintainability & Code Quality | 🟡 Minor | ⚡ Quick win Rename The bit pattern 🤖 Prompt for AI Agents |
||
| ), | ||
| ("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); | ||
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -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) { | ||
| 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]` 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); | ||
| } | ||
|
Comment on lines
+708
to
+717
There was a problem hiding this comment. 📐 Maintainability & Code Quality | 🟡 Minor | ⚡ Quick win 🧩 Analysis chain🏁 Script executed: #!/bin/bash
# Description: Check whether other new beta tests use bit tolerance and locate the sf subnormal evaluation path.
set -euo pipefail
rg -n 'abs_diff\(.*\)\s*<=' --type=rust -g '**/beta*'
rg -n -C8 'checked_ln_beta_reg_complement' --type=rustRepository: statrs-dev/statrs Length of output: 155 🏁 Script executed: #!/bin/bash
set -euo pipefail
printf '%s\n' '--- beta.rs relevant symbols ---'
rg -n -C12 'test_sf_tiny_shape_preserves_representable_tail|checked_ln_beta_reg_complement|fn sf|impl Continuous|abs_diff|to_bits' src/distribution/beta.rs
printf '%s\n' '--- beta.rs file size ---'
wc -l src/distribution/beta.rs
printf '%s\n' '--- related tolerance assertions ---'
rg -n -C3 'abs_diff|to_bits\(\)|assert_relative_eq|assert_eq!' src/distribution --glob '*.rs' | head -n 300
printf '%s\n' '--- platform and math dependency configuration ---'
rg -n -C3 'libm|no_std|target|portable|rustversion|special|num-traits' Cargo.toml Cargo.lock .github src --glob '*.toml' --glob '*.yml' --glob '*.yaml' --glob '*.rs' | head -n 300Repository: statrs-dev/statrs Length of output: 45048 🏁 Script executed: #!/bin/bash
set -euo pipefail
printf '%s\n' '--- beta complement implementation ---'
rg -n -C20 'pub fn checked_ln_beta_reg_complement|fn checked_ln_beta_reg_complement|checked_ln_beta_reg_complement' src/function src
printf '%s\n' '--- beta function math calls ---'
rg -n -C4 '\.(ln|exp|log1p|expm1)\(\)|f64::(ln|exp)' src/function/beta src/function --glob '*.rs' | head -n 300
printf '%s\n' '--- feature-dependent floating-point operations ---'
rg -n -C5 '#\[cfg\(.*std|Float as _|use std|use core' src/function src/distribution/beta.rs --glob '*.rs' | head -n 300
printf '%s\n' '--- manifest and supported target jobs ---'
sed -n '1,90p' Cargo.toml
sed -n '1,145p' .github/workflows/test.yml
printf '%s\n' '--- beta-related recent changes ---'
git diff --stat -- src/distribution/beta.rs src/function/beta
git diff -- src/distribution/beta.rs src/function/beta | sed -n '1,260p'Repository: statrs-dev/statrs Length of output: 50375 🏁 Script executed: #!/bin/bash
set -euo pipefail
printf '%s\n' '--- complete complement dispatch ---'
sed -n '18,115p' src/function/beta/log_forward.rs
printf '%s\n' '--- beta test reference assertions ---'
sed -n '1,120p' src/function/beta/tests.rs
rg -n -C4 'abs_diff|to_bits|reference|500_digit' src/distribution/beta.rs src/function/beta/tests.rs
printf '%s\n' '--- supported test matrix ---'
sed -n '35,75p' .github/workflows/test.yml
printf '%s\n' '--- exact input branch conditions ---'
python3 - <<'PY'
import struct
def f(bits):
return struct.unpack(">d", bits.to_bytes(8, "big"))[0]
a = f(0x00000000000007e8)
b = f(0x4040000000000000)
x = f(0x01556e1fc2f8f359)
threshold = (a + 1.0) / (a + b + 2.0)
print(f"a={a!r}")
print(f"b={b!r}")
print(f"x={x!r}")
print(f"threshold={threshold!r}")
print(f"x<threshold: {x < threshold}")
print(f"x**a: {x**a!r}")
print(f"x**a > 0.5: {x**a > 0.5}")
PYRepository: statrs-dev/statrs Length of output: 50375 🏁 Script executed: #!/bin/bash
set -euo pipefail
printf '%s\n' '--- exact tiny-shape test and adjacent reference ---'
sed -n '1148,1202p' src/function/beta/tests.rs
sed -n '700,718p' src/distribution/beta.rs
printf '%s\n' '--- checked_ln_beta_reg dispatch ---'
rg -n 'pub\(super\) fn checked_ln_beta_reg|fn checked_ln_beta_reg|beta_small_shapes_series_log|beta_power_series_log_parts|use_beta_power_series' src/function/beta --glob '*.rs'
sed -n '190,275p' src/function/beta/log_forward.rs
sed -n '1,145p' src/function/beta/forward.rs
printf '%s\n' '--- test matrix only ---'
sed -n '35,70p' .github/workflows/test.yml
printf '%s\n' '--- math-operation traits and feature selection ---'
sed -n '1,90p' src/function/beta/mod.rs
sed -n '1,35p' src/distribution/beta.rs
rg -n -C3 'fn beta_small_shapes_series_log|fn beta_power_series_log_parts|fn beta_power_series_log' src/function/beta --glob '*.rs'Repository: statrs-dev/statrs Length of output: 16335 🏁 Script executed: #!/bin/bash
set -euo pipefail
printf '%s\n' '--- checked logarithmic beta dispatch ---'
sed -n '90,155p' src/function/beta/mod.rs
printf '%s\n' '--- power-series implementation ---'
sed -n '1,145p' src/function/beta/series.rs
printf '%s\n' '--- small-shape implementation ---'
sed -n '1,180p' src/function/beta/small_gamma.rs
printf '%s\n' '--- exact floating-point operations on the tested path ---'
rg -n -C4 'beta_power_series_log_parts_with_log_beta|beta_small_shapes_series_log_with_log_beta|accurate_ln|ln_beta_accurate_parts|\.ln\(\)|\.exp\(\)|\.powf\(' src/function/beta/{mod.rs,series.rs,small_gamma.rs,log_forward.rs,log_beta.rs}Repository: statrs-dev/statrs Length of output: 31783 🏁 Script executed: #!/bin/bash
set -euo pipefail
printf '%s\n' '--- accurate_ln implementation ---'
rg -n -C12 'pub\(super\) fn accurate_ln|fn accurate_ln|accurate_ln_one_plus_dd' src/function/beta/dd.rs
printf '%s\n' '--- complement finalization ---'
sed -n '120,190p' src/function/beta/log_forward.rs
printf '%s\n' '--- exact branch-condition probe ---'
python3 - <<'PY'
import struct, math
def f(bits):
return struct.unpack(">d", bits.to_bytes(8, "big"))[0]
a = f(0x00000000000007e8)
b = f(0x4040000000000000)
x = f(0x01556e1fc2f8f359)
y = 1.0 - x
scaled_b = b * x
threshold = (a + 1.0) / (a + b + 2.0)
symmetry_split = ((a + 1.0) / max(a + 1.0, b + 1.0)) / (
(a + 1.0) / max(a + 1.0, b + 1.0)
+ (b + 1.0) / max(a + 1.0, b + 1.0)
)
print("a_bits", hex(a.hex().__hash__() & 0)) # keep output deterministic below
print("a", a.hex(), "b", b.hex(), "x", x.hex())
print("x_lt_sf_split", x < threshold)
print("x_pow_a_gt_half", x ** a > 0.5)
print("x_pow_a_gt_09", x ** a > 0.9)
print("y_is_one", y == 1.0)
print("scaled_b", scaled_b.hex())
print("symmetry_split", symmetry_split.hex())
print("x_ge_symmetry_split", x >= symmetry_split)
print("special_complement_branch",
a <= math.sqrt(2.220446049250313e-16)
and b >= 32.0
and x ** a > 0.5)
PYRepository: statrs-dev/statrs Length of output: 5317 🌐 Web query:
💡 Result: In the Rust standard library, mathematical functions such as f64::ln, f64::exp, and f64::exp_m1 are documented as having "unspecified precision" [1][2][3]. The official documentation explicitly states that the precision of these functions is non-deterministic, meaning it can vary based on the platform, the Rust version, and even within the same execution across different invocations [1][3][4]. These functions are typically implemented by calling into the underlying platform's math library (e.g., libc's math functions or equivalent system-provided libraries) [5][6]. Because these system libraries vary across different operating systems and architectures, the exact results of floating-point operations can differ between platforms [5]. Key points regarding this behavior include: Non-deterministic Precision: The standard library documentation warns that because the implementation often relies on external, platform-specific libraries, it does not guarantee a specific level of precision or bit-for-bit identical results across different systems [1][3][4]. Implementation Source: The Rust standard library typically links against the system's C math library (libm) to provide these transcendental functions [3][5][6]. Consequently, observed discrepancies in output (e.g., between Linux and macOS) are often attributable to differences in the underlying libc math implementation [5]. Consistency Requirements: If a project requires strictly deterministic or reproducible floating-point results across different platforms, developers are often advised to use a portable, pure-Rust implementation of these functions, such as the libm crate, which provides consistent results regardless of the host operating system [5]. While the "unspecified precision" note in the documentation highlights that results are not strictly guaranteed to be identical across environments, once a program is compiled and run on identical hardware with the same environment, the operations themselves are generally deterministic [7]. Differences arise primarily when moving between different platforms or build environments where the underlying math library or optimization strategies may change [7][5]. Citations:
Use a ULP tolerance for this The tiny-shape path reaches 🤖 Prompt for AI Agents |
||
|
|
||
| #[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); | ||
|
|
||
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Oops, something went wrong.
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
📐 Maintainability & Code Quality | 🟡 Minor | ⚡ Quick win
🧩 Analysis chain
🏁 Script executed:
Repository: statrs-dev/statrs
Length of output: 3420
🌐 Web query:
Cargo manifest reference include README.md always included explicit include list💡 Result:
In Cargo, the README.md file is not automatically included in the package manifest's explicit
includelist, but its inclusion is handled through specific mechanisms [1][2]. To ensure a README file is included when publishing a crate, you should specify the path to it using thereadmefield in the[package]section of yourCargo.toml[1][3]: [package] #... readme = "README.md" When this field is set, Cargo automatically includes the specified file in the package, regardless of whether it is listed in theincludefield [1][4]. If you do not explicitly use theincludefield, Cargo defaults to including all files in the package root (excluding those ignored by version control) [1][3]. If you choose to use theincludefield to explicitly define which files are bundled with your crate, it acts as an allowlist [1][3]. Ifincludeis present, you must ensure that your README file (along with other necessary files like your source code) is explicitly included in that list [1][2]. For example: [package] include = ["/src", "README.md", "Cargo.toml"] You can always verify exactly which files will be included in your package by running the following command in your terminal [1][5]: cargo package --list This command provides an authoritative list of the files that will be bundled, allowing you to confirm that your README.md is present [1][6].Citations:
🏁 Script executed:
Repository: statrs-dev/statrs
Length of output: 13905
🌐 Web query:
site:doc.rust-lang.org/cargo/reference/manifest.html Cargo include readme field automatically included when include specified💡 Result:
No, the README file is not automatically included when the include field is specified in your Cargo.toml [1]. When you define an include field, you are explicitly specifying the set of files to be included in your package [1]. Unlike the default behavior where Cargo automatically discovers and includes files, specifying include overrides this and means that only the files or patterns you list will be part of the package [1]. While Cargo does have a list of files that are always included regardless of the include field (such as Cargo.toml, a minimized Cargo.lock, and any file specified by license-file), the README file is not on this list [1]. Therefore, if you use the include field, you must explicitly add your README file (e.g., "README.md") to the list of included patterns to ensure it is packaged [1].
Citations:
Add
README.mdto the package allowlist.The explicit
includelist does not automatically includeREADME.md. Add"README.md"so Cargo publishes the README.🤖 Prompt for AI Agents