From 743fa96d00388ec61445cb470b7c4226cf5b5880 Mon Sep 17 00:00:00 2001 From: Ciaran Ryan-Anderson Date: Mon, 6 Jul 2026 11:07:05 -0600 Subject: [PATCH 01/11] Add libm/puruspe deps and scipy oracle fixtures for Jeffreys interval work --- Cargo.lock | 22 +++++ Cargo.toml | 3 + crates/pecos-num/Cargo.toml | 7 ++ .../fixtures/generate_jeffreys_fixtures.py | 87 +++++++++++++++++++ .../tests/fixtures/jeffreys_scipy.csv | 76 ++++++++++++++++ 5 files changed, 195 insertions(+) create mode 100644 crates/pecos-num/tests/fixtures/generate_jeffreys_fixtures.py create mode 100644 crates/pecos-num/tests/fixtures/jeffreys_scipy.csv diff --git a/Cargo.lock b/Cargo.lock index 2398038d6..8102b0bca 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -3207,6 +3207,16 @@ version = "0.20.2" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "507460a910eb7b32ee961886ff48539633b788a36b65692b95f225b844c82553" +[[package]] +name = "lambert_w" +version = "1.2.34" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f5f0846ee4f0299ca4c5b9ca06ff55cf88b3430a763bf591474cc734479c9b24" +dependencies = [ + "num-complex 0.4.6", + "num-traits", +] + [[package]] name = "lazy_static" version = "1.5.0" @@ -4538,12 +4548,14 @@ dependencies = [ name = "pecos-num" version = "0.2.0-dev.0" dependencies = [ + "libm", "log", "nalgebra", "ndarray 0.17.2", "num-complex 0.4.6", "num-traits", "pecos-random", + "puruspe", "rand 0.10.1", "rustworkx-core", "serde", @@ -5481,6 +5493,16 @@ dependencies = [ "syn 2.0.117", ] +[[package]] +name = "puruspe" +version = "0.4.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2c9d62b4ad8b456f2ac8a171567bd475cae8cecd53ff5a46fce4f261eef17648" +dependencies = [ + "lambert_w", + "num-complex 0.4.6", +] + [[package]] name = "py_literal" version = "0.4.0" diff --git a/Cargo.toml b/Cargo.toml index e659998e4..a4a05d278 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -89,8 +89,11 @@ pest_derive = "2" regex = "1" # --- Numerical computing --- +libm = "0.2" nalgebra = "0.35" num = "0.4" +# Dev-only differential oracle for special-function tests (see design/jeffreys-pecos-num.md) +puruspe = "0.4" num-complex = "0.4" num-traits = "0.2" num-bigint = { version = "0.4", features = ["serde"] } diff --git a/crates/pecos-num/Cargo.toml b/crates/pecos-num/Cargo.toml index b34dab9ae..9157d0bac 100644 --- a/crates/pecos-num/Cargo.toml +++ b/crates/pecos-num/Cargo.toml @@ -38,5 +38,12 @@ rustworkx-core.workspace = true serde = { workspace = true, features = ["derive"] } serde_json.workspace = true +# Deterministic libm (musl-derived) for special functions; see design/jeffreys-pecos-num.md +libm.workspace = true + +[dev-dependencies] +# Differential oracle for special-function tests only -- never a runtime dependency +puruspe.workspace = true + [lints] workspace = true diff --git a/crates/pecos-num/tests/fixtures/generate_jeffreys_fixtures.py b/crates/pecos-num/tests/fixtures/generate_jeffreys_fixtures.py new file mode 100644 index 000000000..fabe3d6ab --- /dev/null +++ b/crates/pecos-num/tests/fixtures/generate_jeffreys_fixtures.py @@ -0,0 +1,87 @@ +# Copyright 2026 The PECOS Developers +# +# Licensed under the Apache License, Version 2.0 (the "License"); you may +# not use this file except in compliance with the License. You may obtain a +# copy of the License at +# +# https://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, WITHOUT +# WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the +# License for the specific language governing permissions and limitations +# under the License. +"""Offline generator for the Jeffreys-interval oracle fixture table. + +scipy.stats.beta.ppf (TOMS-708 lineage) is the oracle of record per +design/jeffreys-pecos-num.md in the pecos-docs vault. Run offline with: + + uv run --with scipy python generate_jeffreys_fixtures.py > jeffreys_scipy.csv + +scipy is never a runtime dependency; this script exists so the pinned CSV is +reproducible. The adversarial case grid follows the design note's v3 fixture +table (endpoint boundaries, n=1e8 band, typical-LER regime, alpha extremes +within the supported regime alpha in [1e-6, 0.5]). +""" + +from __future__ import annotations + +import sys + +from scipy.stats import beta + + +def jeffreys_row(k: int, n: int, alpha: float) -> str: + a = k + 0.5 + b = n - k + 0.5 + lo = 0.0 if k == 0 else float(beta.ppf(alpha / 2, a, b)) + hi = 1.0 if k == n else float(beta.ppf(1 - alpha / 2, a, b)) + median = float(beta.ppf(0.5, a, b)) + return f"{k},{n},{alpha!r},{lo!r},{hi!r},{median!r}" + + +def cases() -> list[tuple[int, int, float]]: + out: list[tuple[int, int, float]] = [] + # Minimum and smallest non-trivial experiments + for n in (1, 2): + for k in range(n + 1): + out.append((k, n, 0.05)) + # Endpoint boundaries k=0 and k=n across the full supported regime + for n in (100, 1_000, 10_000, 100_000, 1_000_000, 10_000_000, 100_000_000): + for alpha in (0.01, 0.05, 0.1): + out.append((0, n, alpha)) + out.append((n, n, alpha)) + # Near-boundary shapes a=1.5 / b=1.5 at large n + for n in (1_000_000, 100_000_000): + out.append((1, n, 0.05)) + out.append((n - 1, n, 0.05)) + # n=1e8 boundary band (v3 addition), all in-regime alphas + n8 = 100_000_000 + for k in (0, 1, 10, n8 - 10, n8 - 1, n8): + for alpha in (0.01, 0.05, 0.1, 1e-6): + out.append((k, n8, alpha)) + # Typical logical-error-rate working regime + for n in (1_000, 10_000, 100_000, 1_000_000): + out.append((max(1, n // 1000), n, 0.05)) + # Wide-CI alpha extreme and best-conditioned symmetric peak + for n in (1_000, 1_000_000): + out.append((n // 2, n, 1e-6)) + out.append((n // 2, n, 0.05)) + # Deduplicate, preserving order + seen: set[tuple[int, int, float]] = set() + unique = [] + for case in out: + if case not in seen: + seen.add(case) + unique.append(case) + return unique + + +def main() -> None: + print("k,n,alpha,lo,hi,median") + for k, n, alpha in cases(): + print(jeffreys_row(k, n, alpha)) + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/crates/pecos-num/tests/fixtures/jeffreys_scipy.csv b/crates/pecos-num/tests/fixtures/jeffreys_scipy.csv new file mode 100644 index 000000000..44f86df43 --- /dev/null +++ b/crates/pecos-num/tests/fixtures/jeffreys_scipy.csv @@ -0,0 +1,76 @@ +k,n,alpha,lo,hi,median +0,1,0.05,0.0,0.8532536836904248,0.16319398540839275 +1,1,0.05,0.14674631630957513,1.0,0.8368060145916072 +0,2,0.05,0.0,0.6668217544012089,0.09552581803782106 +1,2,0.05,0.06083027592009732,0.9391697240799026,0.5 +2,2,0.05,0.3331782455987909,1.0,0.9044741819621789 +0,100,0.01,0.0,0.03853633289959721,0.002266429228530002 +100,100,0.01,0.9614636671004028,1.0,0.99773357077147 +0,100,0.05,0.0,0.02474527001526989,0.002266429228530002 +100,100,0.05,0.9752547299847301,1.0,0.99773357077147 +0,100,0.1,0.0,0.01897688770307349,0.002266429228530002 +100,100,0.1,0.9810231122969265,1.0,0.99773357077147 +0,1000,0.01,0.0,0.003930987519119225,0.00022738549453875173 +1000,1000,0.01,0.9960690124808808,1.0,0.9997726145054613 +0,1000,0.05,0.0,0.0025081643414222978,0.00022738549453875173 +1000,1000,0.05,0.9974918356585777,1.0,0.9997726145054613 +0,1000,0.1,0.0,0.001918406712256031,0.00022738549453875173 +1000,1000,0.1,0.9980815932877439,1.0,0.9997726145054613 +0,10000,0.01,0.0,0.00039388448646258923,2.2745993797440672e-05 +10000,10000,0.01,0.9996061155135374,1.0,0.9999772540062025 +0,10000,0.05,0.0,0.0002511564843832975,2.2745993797440672e-05 +10000,10000,0.05,0.9997488435156167,1.0,0.9999772540062025 +0,10000,0.1,0.0,0.00019204969529015292,2.2745993797440672e-05 +10000,10000,0.1,0.9998079503047098,1.0,0.9999772540062025 +0,100000,0.01,0.0,3.939631833459792e-05,2.2746738418241375e-06 +100000,100000,0.01,0.9999606036816654,1.0,0.9999977253261582 +0,100000,0.05,0.0,2.5119052649257954e-05,2.2746738418241375e-06 +100000,100000,0.05,0.9999748809473508,1.0,0.9999977253261582 +0,100000,0.1,0.0,1.920706162724842e-05,2.2746738418241375e-06 +100000,100000,0.1,0.9999807929383727,1.0,0.9999977253261582 +0,1000000,0.01,0.0,3.939710542701222e-06,2.2746812882186074e-07 +1000000,1000000,0.01,0.9999960602894573,1.0,0.9999997725318712 +0,1000000,0.05,0.0,2.5119393107467838e-06,2.2746812882186074e-07 +1000000,1000000,0.05,0.9999974880606892,1.0,0.9999997725318712 +0,1000000,0.1,0.0,1.9207270855660624e-06,2.2746812882186074e-07 +1000000,1000000,0.1,0.9999980792729144,1.0,0.9999997725318712 +0,10000000,0.01,0.0,3.9397184137489613e-07,2.2746820328599193e-08 +10000000,10000000,0.01,0.9999996060281586,1.0,0.9999999772531797 +0,10000000,0.05,0.0,2.5119427153660024e-07,2.2746820328599193e-08 +10000000,10000000,0.05,0.9999997488057285,1.0,0.9999999772531797 +0,10000000,0.1,0.0,1.9207291778687745e-07,2.2746820328599193e-08 +10000000,10000000,0.1,0.9999998079270822,1.0,0.9999999772531797 +0,100000000,0.01,0.0,3.939719200854972e-08,2.27468210732407e-09 +100000000,100000000,0.01,0.999999960602808,1.0,0.9999999977253179 +0,100000000,0.05,0.0,2.5119430558282962e-08,2.27468210732407e-09 +100000000,100000000,0.05,0.9999999748805695,1.0,0.9999999977253179 +0,100000000,0.1,0.0,1.920729387099232e-08,2.27468210732407e-09 +100000000,100000000,0.1,0.9999999807927061,1.0,0.9999999977253179 +1,1000000,0.05,1.0789766246541867e-07,4.6741920467298275e-06,1.1829865382054428e-06 +999999,1000000,0.05,0.9999953258079533,0.9999998921023375,0.9999988170134618 +1,100000000,0.05,1.0789764152348355e-09,4.674201704692767e-08,1.1829869381478459e-08 +99999999,100000000,0.05,0.999999953257983,0.9999999989210235,0.9999999881701306 +0,100000000,1e-06,0.0,1.2631909533627897e-07,2.27468210732407e-09 +1,100000000,0.01,3.58608873764682e-10,6.419078043324197e-08,1.1829869381478459e-08 +1,100000000,0.1,1.759231591596988e-09,3.907363885056535e-08,1.1829869381478459e-08 +1,100000000,1e-06,7.616416785435591e-13,1.6047017918167062e-07,1.1829869381478459e-08 +10,100000000,0.01,4.016826820241171e-08,2.0700531226423903e-07,1.0168613747779586e-07 +10,100000000,0.05,5.1414490033077726e-08,1.7739437222048644e-07,1.0168613747779586e-07 +10,100000000,0.1,5.795652711755923e-08,1.633528611217682e-07,1.0168613747779586e-07 +10,100000000,1e-06,1.3296244027164594e-08,3.452114945345266e-07,1.0168613747779586e-07 +99999990,100000000,0.01,0.9999997929946878,0.9999999598317318,0.9999998983138625 +99999990,100000000,0.05,0.9999998226056278,0.99999994858551,0.9999998983138625 +99999990,100000000,0.1,0.9999998366471389,0.9999999420434729,0.9999998983138625 +99999990,100000000,1e-06,0.9999996547885055,0.999999986703756,0.9999998983138625 +99999999,100000000,0.01,0.9999999358092195,0.9999999996413911,0.9999999881701306 +99999999,100000000,0.1,0.9999999609263611,0.9999999982407685,0.9999999881701306 +99999999,100000000,1e-06,0.9999998395298209,0.9999999999992384,0.9999999881701306 +100000000,100000000,1e-06,0.9999998736809047,1.0,0.9999999977253179 +1,1000,0.05,0.00010791880468097968,0.004664458808791372,0.0011825830506087232 +10,10000,0.05,0.0005142569824508165,0.001773213021453329,0.0010168273845347388 +100,100000,0.05,0.0008182225087040382,0.0012106486647282626,0.0010016652994178298 +1000,1000000,0.05,0.0009394857355159988,0.00106340264158502,0.0010001663530183856 +500,1000,1e-06,0.42313580559571984,0.5768641944045314,0.5 +500,1000,0.05,0.46904771269911394,0.5309522873008861,0.5 +500000,1000000,1e-06,0.49755419600449113,0.5024458039955171,0.5 +500000,1000000,0.05,0.49902001919386424,0.5009799808061357,0.5 From f2ff279f4ed56184e6945c9d71087a7ca2e24dbf Mon Sep 17 00:00:00 2001 From: Ciaran Ryan-Anderson Date: Mon, 6 Jul 2026 13:16:47 -0600 Subject: [PATCH 02/11] Extend Jeffreys oracle fixtures with GL-asymmetric and CF-band rows from adversarial review --- .../fixtures/generate_jeffreys_fixtures.py | 13 ++++++++ .../tests/fixtures/jeffreys_scipy.csv | 32 +++++++++++++++++++ 2 files changed, 45 insertions(+) diff --git a/crates/pecos-num/tests/fixtures/generate_jeffreys_fixtures.py b/crates/pecos-num/tests/fixtures/generate_jeffreys_fixtures.py index fabe3d6ab..e6cad4280 100644 --- a/crates/pecos-num/tests/fixtures/generate_jeffreys_fixtures.py +++ b/crates/pecos-num/tests/fixtures/generate_jeffreys_fixtures.py @@ -67,6 +67,19 @@ def cases() -> list[tuple[int, int, float]]: for n in (1_000, 1_000_000): out.append((n // 2, n, 1e-6)) out.append((n // 2, n, 0.05)) + # Asymmetric Gauss-Legendre-region shapes (both Beta shapes > 3000), including + # the k=3016, n=1e6 tail-orientation regression and a large-n asymmetric case + for k in (3000, 3016, 3100, 5000, 10_000): + for alpha in (0.01, 0.05, 1e-6): + out.append((k, 1_000_000, alpha)) + out.append((15_881, 20_000, 0.05)) + out.append((10_000, 100_000_000, 0.05)) + out.append((67_867_393, 100_000_000, 0.05)) + # Continued-fraction accuracy band: moderate a with huge b, straddling the + # asymptotic-branch gate at a = 64.5 + for k in (63, 64, 65, 100, 300, 1_000, 3_000): + for alpha in (0.01, 0.05): + out.append((k, 100_000_000, alpha)) # Deduplicate, preserving order seen: set[tuple[int, int, float]] = set() unique = [] diff --git a/crates/pecos-num/tests/fixtures/jeffreys_scipy.csv b/crates/pecos-num/tests/fixtures/jeffreys_scipy.csv index 44f86df43..a534da554 100644 --- a/crates/pecos-num/tests/fixtures/jeffreys_scipy.csv +++ b/crates/pecos-num/tests/fixtures/jeffreys_scipy.csv @@ -74,3 +74,35 @@ k,n,alpha,lo,hi,median 500,1000,0.05,0.46904771269911394,0.5309522873008861,0.5 500000,1000000,1e-06,0.49755419600449113,0.5024458039955171,0.5 500000,1000000,0.05,0.49902001919386424,0.5009799808061357,0.5 +3000,1000000,0.01,0.0028614808676996543,0.003143247046939454,0.003000165673185513 +3000,1000000,0.05,0.002894242329240523,0.0031086345614359566,0.003000165673185513 +3000,1000000,1e-06,0.002740507165346054,0.0032756789001165063,0.003000165673185513 +3016,1000000,0.01,0.0028771067966790976,0.0031596209664586815,0.003016165667817254 +3016,1000000,0.05,0.0029099577001540575,0.003124919098145593,0.003016165667817254 +3016,1000000,1e-06,0.002755796789951582,0.0032923887625396072,0.003016165667817254 +3100,1000000,0.01,0.0029591591419661813,0.003245567825669476,0.0031001656396398066 +3100,1000000,0.05,0.0029924757389128663,0.0032104005743676463,0.0031001656396398066 +3100,1000000,1e-06,0.0028360981029479115,0.003380084755062128,0.0031001656396398066 +5000,1000000,0.01,0.004820663512914193,0.005184045428870807,0.005000165003885435 +5000,1000000,0.05,0.004863184941630663,0.0051396803899975965,0.005000165003885435 +5000,1000000,1e-06,0.004662988700966737,0.005353132838454261,0.005000165003885435 +10000,1000000,0.01,0.009746033453896267,0.010258627962139022,0.010000163335243898 +10000,1000000,0.05,0.009806401413221992,0.010196434989008803,0.010000163335243898 +10000,1000000,1e-06,0.009521236717447784,0.010494722430337888,0.010000163335243898 +15881,20000,0.05,0.7884031796533509,0.7996117205822614,0.7940450990562469 +10000,100000000,0.05,9.80545708146841e-05,0.00010197436631891206,0.0001000016663530796 +67867393,100000000,0.05,0.6785823973179447,0.6787654523393154,0.6786739294044202 +63,100000000,0.01,4.48522072486598e-07,8.589803669099042e-07,6.316697931915235e-07 +63,100000000,0.05,4.884924603417635e-07,8.004289916647383e-07,6.316697931915235e-07 +64,100000000,0.01,4.5691245771068727e-07,8.705909692320621e-07,6.416697443785606e-07 +64,100000000,0.05,4.972659791123161e-07,8.116558085874771e-07,6.416697443785606e-07 +65,100000000,0.01,4.6531526399032663e-07,8.821891201850533e-07,6.516696970632925e-07 +65,100000000,0.05,5.060489905109922e-07,8.22873122558363e-07,6.516696970632925e-07 +100,100000000,0.01,7.655619354259218e-07,1.2819639049330916e-06,1.0016686359541935e-06 +100,100000000,0.05,8.181502132297306e-07,1.210779186771337e-06,1.0016686359541935e-06 +300,100000000,0.01,2.577273345719969e-06,3.4702789793604664e-06,3.001667314799613e-06 +300,100000000,0.05,2.674810323215624e-06,3.3541280237101947e-06,3.001667314799613e-06 +1000,100000000,0.01,9.209040467852042e-06,1.0838520545443545e-05,1.0001666830830989e-05 +1000,100000000,0.05,9.394578061817165e-06,1.0634363055929011e-05,1.0001666830830989e-05 +3000,100000000,0.01,2.8612849477460176e-05,3.143471230750312e-05,3.0001666632500787e-05 +3000,100000000,0.05,2.8940912004756113e-05,3.1088028863089894e-05,3.0001666632500787e-05 From 7039aa874e3b15f201e264b6723eeaca65ad92d9 Mon Sep 17 00:00:00 2001 From: Ciaran Ryan-Anderson Date: Mon, 6 Jul 2026 14:21:24 -0600 Subject: [PATCH 03/11] Add mid-band CF and median-gate oracle fixtures from round-2 adversarial review --- .../fixtures/generate_jeffreys_fixtures.py | 11 +++++++ .../tests/fixtures/jeffreys_scipy.csv | 31 +++++++++++++++++++ 2 files changed, 42 insertions(+) diff --git a/crates/pecos-num/tests/fixtures/generate_jeffreys_fixtures.py b/crates/pecos-num/tests/fixtures/generate_jeffreys_fixtures.py index e6cad4280..8cf283fde 100644 --- a/crates/pecos-num/tests/fixtures/generate_jeffreys_fixtures.py +++ b/crates/pecos-num/tests/fixtures/generate_jeffreys_fixtures.py @@ -80,6 +80,17 @@ def cases() -> list[tuple[int, int, float]]: for k in (63, 64, 65, 100, 300, 1_000, 3_000): for alpha in (0.01, 0.05): out.append((k, 100_000_000, alpha)) + # Mid-band b in [1e7, 5e7): small-a upper bounds routed through the CF path + # missed the first asymptotic gate (round-2 review); mirrors pin the lo side + for n in (10_000_000, 20_000_000, 40_000_000, 49_999_990): + for k in (10, 64, 300): + for alpha in (0.01, 0.05): + out.append((k, n, alpha)) + out.append((n - 64, n, 0.05)) + # Median-shortcut gate edges (min posterior shape near 1e5 and 1e6) + out.append((100_000, 1_000_000, 0.05)) + out.append((100_000, 100_000_000, 0.05)) + out.append((1_000_000, 10_000_000, 0.05)) # Deduplicate, preserving order seen: set[tuple[int, int, float]] = set() unique = [] diff --git a/crates/pecos-num/tests/fixtures/jeffreys_scipy.csv b/crates/pecos-num/tests/fixtures/jeffreys_scipy.csv index a534da554..4f1e88fc6 100644 --- a/crates/pecos-num/tests/fixtures/jeffreys_scipy.csv +++ b/crates/pecos-num/tests/fixtures/jeffreys_scipy.csv @@ -106,3 +106,34 @@ k,n,alpha,lo,hi,median 1000,100000000,0.05,9.394578061817165e-06,1.0634363055929011e-05,1.0001666830830989e-05 3000,100000000,0.01,2.8612849477460176e-05,3.143471230750312e-05,3.0001666632500787e-05 3000,100000000,0.05,2.8940912004756113e-05,3.1088028863089894e-05,3.0001666632500787e-05 +10,10000000,0.01,4.016827811364708e-07,2.0700520792861373e-06,1.0168613441830198e-06 +10,10000000,0.05,5.141450011725289e-07,1.773943064471412e-06,1.0168613441830198e-06 +64,10000000,0.01,4.569128238785155e-06,8.705900462656166e-06,6.416697251195756e-06 +64,10000000,0.05,4.9726628732008384e-06,8.116551633597828e-06,6.416697251195756e-06 +300,10000000,0.01,2.5772781904541334e-05,3.4702715572205535e-05,3.0016672247407065e-05 +300,10000000,0.05,2.674814177267023e-05,3.354122603137294e-05,3.0016672247407065e-05 +9999936,10000000,0.05,0.9999918834483664,0.9999950273371268,0.9999935833027488 +10,20000000,0.01,2.0084136303702057e-07,1.0350263294642324e-06,5.084306805901041e-07 +10,20000000,0.05,2.570724725746622e-07,8.869717149394389e-07,5.084306805901041e-07 +64,20000000,0.01,2.28456310225881e-06,4.352952795124541e-06,3.2083486790950578e-06 +64,20000000,0.05,2.4863305804670033e-06,4.058277609098738e-06,3.2083486790950578e-06 +300,20000000,0.01,1.2886377494627323e-05,1.73513784032059e-05,1.5008336373867161e-05 +300,20000000,0.05,1.3374060180599496e-05,1.6770628072870463e-05,1.5008336373867161e-05 +19999936,20000000,0.05,0.9999959417223909,0.9999975136694196,0.9999967916513209 +10,40000000,0.01,1.0042067463570894e-07,5.175132371874155e-07,2.5421534241970054e-07 +10,40000000,0.05,1.2853622928443265e-07,4.4348590314565433e-07,2.5421534241970054e-07 +64,40000000,0.01,1.14228129684639e-06,2.1764770385110207e-06,1.6041743529218244e-06 +64,40000000,0.05,1.2431650762004882e-06,2.0291392526240094e-06,1.6041743529218244e-06 +300,40000000,0.01,6.443185382925531e-06,8.675694355856318e-06,7.50416824947449e-06 +300,40000000,0.05,6.6870274138833214e-06,8.385317800713864e-06,7.50416824947449e-06 +39999936,40000000,0.05,0.9999979708607474,0.9999987568349238,0.9999983958256471 +10,49999990,0.01,8.033655467463092e-08,4.140106841449082e-07,2.033723149501669e-07 +10,49999990,0.05,1.0282900287288403e-07,3.5478880078242726e-07,2.033723149501669e-07 +64,49999990,0.01,9.138251795569253e-07,1.7411820815969277e-06,1.283339741145293e-06 +64,49999990,0.05,9.945322256216035e-07,1.623311798453432e-06,1.283339741145293e-06 +300,49999990,0.01,5.15454879894997e-06,6.9405576974750375e-06,6.00333581025329e-06 +300,49999990,0.05,5.349622572806005e-06,6.708256184505317e-06,6.00333581025329e-06 +49999926,49999990,0.05,0.9999983766882016,0.9999990054677744,0.9999987166602589 +100000,1000000,0.05,0.09941316855795058,0.10058914688223708,0.10000013333347238 +100000,100000000,0.05,0.0009938195820597677,0.0010062093030989007,0.0010000016633353022 +1000000,10000000,0.05,0.09981417726360871,0.1001860542808535,0.10000001333333473 From 480b345f6385be5fd82230dd8fea26b9109a7f07 Mon Sep 17 00:00:00 2001 From: Ciaran Ryan-Anderson Date: Mon, 6 Jul 2026 14:54:46 -0600 Subject: [PATCH 04/11] Add small-a below-gate pocket oracle fixtures from round-3 adversarial review --- .../fixtures/generate_jeffreys_fixtures.py | 7 ++++++ .../tests/fixtures/jeffreys_scipy.csv | 25 +++++++++++++++++++ 2 files changed, 32 insertions(+) diff --git a/crates/pecos-num/tests/fixtures/generate_jeffreys_fixtures.py b/crates/pecos-num/tests/fixtures/generate_jeffreys_fixtures.py index 8cf283fde..e369b1ad8 100644 --- a/crates/pecos-num/tests/fixtures/generate_jeffreys_fixtures.py +++ b/crates/pecos-num/tests/fixtures/generate_jeffreys_fixtures.py @@ -91,6 +91,13 @@ def cases() -> list[tuple[int, int, float]]: out.append((100_000, 1_000_000, 0.05)) out.append((100_000, 100_000_000, 0.05)) out.append((1_000_000, 10_000_000, 0.05)) + # Small-a upper-tail pocket below the round-4 BGRAT gate (round-3 review): + # CF served hi endpoints here at up to 1.1e-10; rows span the candidate + # boundary region so both sides of wherever the final gate lands are pinned + for n in (1_000_000, 3_000_000, 6_000_000, 8_900_002, 8_999_990): + for k in (2, 5, 20, 64): + out.append((k, n, 0.05)) + out.append((n - 2, n, 0.05)) # Deduplicate, preserving order seen: set[tuple[int, int, float]] = set() unique = [] diff --git a/crates/pecos-num/tests/fixtures/jeffreys_scipy.csv b/crates/pecos-num/tests/fixtures/jeffreys_scipy.csv index 4f1e88fc6..ece3930a6 100644 --- a/crates/pecos-num/tests/fixtures/jeffreys_scipy.csv +++ b/crates/pecos-num/tests/fixtures/jeffreys_scipy.csv @@ -137,3 +137,28 @@ k,n,alpha,lo,hi,median 100000,1000000,0.05,0.09941316855795058,0.10058914688223708,0.10000013333347238 100000,100000000,0.05,0.0009938195820597677,0.0010062093030989007,0.0010000016633353022 1000000,10000000,0.05,0.09981417726360871,0.1001860542808535,0.10000001333333473 +2,1000000,0.05,4.15606032083811e-07,6.416235225085568e-06,2.175729360444771e-06 +5,1000000,0.05,1.9078765988486266e-06,1.0959989229536586e-05,5.170497303788971e-06 +20,1000000,0.05,1.2607302768672983e-05,3.028012265219301e-05,2.016764072793671e-05 +64,1000000,0.05,4.972693694625552e-05,8.116487110228729e-05,6.416695325297659e-05 +999998,1000000,0.05,0.9999935837647749,0.999999584393968,0.9999978242706395 +2,3000000,0.05,1.3853529395225863e-07,2.138748579900197e-06,7.252432835044664e-07 +5,3000000,0.05,6.359583167867743e-07,3.653337610058977e-06,1.7234994864645048e-06 +20,3000000,0.05,4.2024246007293395e-06,1.0093410487460378e-05,6.7225484054094275e-06 +64,3000000,0.05,1.657556954609152e-05,2.705511635140547e-05,2.138898917296257e-05 +2999998,3000000,0.05,0.9999978612514201,0.9999998614647061,0.9999992747567165 +2,6000000,0.05,6.926764071667458e-08,1.0693747280595378e-06,3.6262166217176255e-07 +5,6000000,0.05,3.179790897065049e-07,1.8266697883904397e-06,8.617497913824508e-07 +20,6000000,0.05,2.1012110934354853e-06,5.046709777484324e-06,3.3612743897168966e-06 +64,6000000,0.05,8.287775260381557e-06,1.352757809021066e-05,1.0694495180894354e-05 +5999998,6000000,0.05,0.9999989306252719,0.9999999307323593,0.9999996373783379 +2,8900002,0.05,4.66972740076304e-08,7.209267171947198e-07,2.4446399146336406e-07 +5,8900002,0.05,2.1436786238388728e-07,1.231462717971291e-06,5.809547955642139e-07 +20,8900002,0.05,1.4165462211133705e-06,3.402276485853845e-06,2.266027210322075e-06 +64,8900002,0.05,5.587260875286885e-06,9.119717891657e-06,7.209770542464782e-06 +8900000,8900002,0.05,0.9999992790732828,0.999999953302726,0.9999997555360085 +2,8999990,0.05,4.617847706288215e-08,7.129173748609068e-07,2.4174804792778066e-07 +5,8999990,0.05,2.119862800809379e-07,1.217781430541664e-06,5.745005099555548e-07 +20,8999990,0.05,1.4008086838723163e-06,3.364477930796822e-06,2.2408521245384656e-06 +64,8999990,0.05,5.525187532107589e-06,9.018399839345175e-06,7.129671507878304e-06 +8999988,8999990,0.05,0.9999992870826251,0.999999953821523,0.9999997582519521 From 261b43a683c07bd22d4925289756da9cc58a2c2e Mon Sep 17 00:00:00 2001 From: Ciaran Ryan-Anderson Date: Mon, 6 Jul 2026 21:19:33 -0600 Subject: [PATCH 05/11] Implement Jeffreys interval and special functions in pecos-num per design note --- crates/pecos-num/src/lib.rs | 5 +- crates/pecos-num/src/prelude.rs | 8 +- crates/pecos-num/src/special.rs | 1580 +++++++++++++++++++++------- crates/pecos-num/src/stats.rs | 364 ++++--- crates/pecos-num/tests/jeffreys.rs | 572 ++++++++++ 5 files changed, 2001 insertions(+), 528 deletions(-) create mode 100644 crates/pecos-num/tests/jeffreys.rs diff --git a/crates/pecos-num/src/lib.rs b/crates/pecos-num/src/lib.rs index 538ef4e85..df965af51 100644 --- a/crates/pecos-num/src/lib.rs +++ b/crates/pecos-num/src/lib.rs @@ -57,5 +57,6 @@ pub use graph::Graph; pub use linalg::{matrix_exp, matrix_log}; pub use optimize::{BrentqOptions, NewtonOptions, OptimizeError, brentq, newton}; pub use polynomial::{Poly1d, PolynomialError, polyfit}; -pub use special::{betainc_inv, betainc_reg, ln_gamma}; -pub use stats::{jeffreys_interval, mean}; +pub use stats::{ + JeffreysError, JeffreysEstimator, JeffreysInterval, jeffreys_interval, jeffreys_point, mean, +}; diff --git a/crates/pecos-num/src/prelude.rs b/crates/pecos-num/src/prelude.rs index 2dafd93d5..55eccc4f7 100644 --- a/crates/pecos-num/src/prelude.rs +++ b/crates/pecos-num/src/prelude.rs @@ -31,13 +31,11 @@ pub use crate::random; // Re-export statistical functions pub use crate::stats::{ - jackknife_resamples, jackknife_stats, jackknife_stats_axis, jackknife_weighted, - jeffreys_interval, mean, mean_axis, std, std_axis, weighted_mean, + JeffreysError, JeffreysEstimator, JeffreysInterval, jackknife_resamples, jackknife_stats, + jackknife_stats_axis, jackknife_weighted, jeffreys_interval, jeffreys_point, mean, mean_axis, + std, std_axis, weighted_mean, }; -// Re-export special functions -pub use crate::special::{betainc_inv, betainc_reg, ln_gamma}; - // Re-export mathematical traits (use these for polymorphism!) pub use crate::math::{ Abs, Acos, Acosh, Asin, Asinh, Atan, Atan2, Atanh, Ceil, Cos, Cosh, Exp, Floor, Ln, LogBase, diff --git a/crates/pecos-num/src/special.rs b/crates/pecos-num/src/special.rs index 214657bbf..ac53df4b2 100644 --- a/crates/pecos-num/src/special.rs +++ b/crates/pecos-num/src/special.rs @@ -12,457 +12,1313 @@ // See the License for the specific language governing permissions and // limitations under the License. -//! Special functions: log-gamma and the regularized incomplete beta -//! function with its inverse. +//! Special functions used by PECOS numerical routines. //! -//! These are the primitives behind Beta-distribution quantiles, which in -//! turn back the Jeffreys binomial interval in [`crate::stats`]. -//! -//! Algorithms follow Numerical Recipes, 3rd edition, section 6.1 -//! (Lanczos log-gamma) and section 6.4 (Lentz continued fraction for the -//! incomplete beta; Halley iteration with the Abramowitz & Stegun 26.5.22 -//! initial guess for the inverse). Differential tests against `SciPy` -//! reference values live in the test module below. - -/// Convergence threshold for the inverse incomplete beta Halley iteration. -const BETAINC_INV_EPS: f64 = 1.0e-8; -/// Maximum iterations for the incomplete-beta continued fraction. -const BETACF_MAX_ITER: u32 = 10_000; -/// Convergence threshold for the continued fraction. -const BETACF_EPS: f64 = 3.0e-14; -/// Smallest representable ratio guard for Lentz's method. -const BETACF_FPMIN: f64 = f64::MIN_POSITIVE / BETACF_EPS; - -/// Natural logarithm of the gamma function for `x > 0`. +//! The incomplete-beta implementation follows the Numerical Recipes-style +//! `ln_gamma` / continued-fraction / regularized-beta / inverse-beta chain, +//! with the reflection branch selected before forming complements. + +use std::fmt; + +const BETACF_MAX_ITERATIONS: usize = 10_000; +const BETACF_EPSILON: f64 = f64::EPSILON; +const BETA_POWER_SERIES_MAX_ITERATIONS: usize = 10_000; +const BETA_POWER_SERIES_EPSILON: f64 = f64::EPSILON; +const GAMMA_MAX_ITERATIONS: usize = 10_000; +const GAMMA_EPSILON: f64 = f64::EPSILON; +/// Minimum `b` for the large-b beta expansion. +/// +/// The effective accuracy bound is the scale test `a*a/b <= 0.1`; this floor +/// only prevents very small-count probes from entering the asymptotic path. +/// Round-5 validation swept the scipy pocket fixtures from `n = 1e6` through +/// `n = 8_999_990` for `a <= 64.5`, plus below-floor ordered probes at +/// `b ~= 1e4` and `b ~= 4.1e4`. +const BETA_LARGE_SHAPE_MIN: f64 = 10_000.0; +const BETA_LARGE_SHAPE_MODERATE_MAX: f64 = 3_001.0; +/// Maximum beta argument for large-b expansion probes. +/// +/// Jeffreys quantiles in the validated large-b branch are well below this; the +/// wider guard keeps inverse-solver boundary probes out of the reflected-CF +/// underflow zone. +const BETA_LARGE_SHAPE_X_MAX: f64 = 1.0e-2; +/// Expansion scale bound for the BGRAT large-b series. +/// +/// The value was retained only after fixture and boundary sweeps covering +/// `a <= 3000.5`, `b <= 1e8`, and the Round-5 small-a pocket rows. +const BETA_LARGE_SHAPE_EXPANSION_PARAMETER_MAX: f64 = 0.1; +const BETA_LARGE_B_EXPANSION_TERMS: usize = 48; +const BETA_LARGE_B_EXPANSION_LEN: usize = BETA_LARGE_B_EXPANSION_TERMS + 1; +const BETA_GAUSS_LEGENDRE_SWITCH: f64 = 3_000.0; +const BETA_GAUSS_LEGENDRE_PANELS: usize = 2; +const BETA_GAUSS_LEGENDRE_PROBABILITY_FLOOR: f64 = 5.0e-15; +const BETA_ASYMPTOTIC_NORMALIZER_MIN: f64 = 100_000.0; +const BETA_MEDIAN_ASYMPTOTIC_MIN: f64 = 1_000_000.0; +const GAMMA_HALF_INTEGER_MAX_X: f64 = 700.0; +const GAMMA_HALF_INTEGER_MAX_STEPS: usize = 128; +const FPMIN: f64 = f64::MIN_POSITIVE / f64::EPSILON; +const HALF_LN_TWO_PI: f64 = 0.918_938_533_204_672_8; +const LN_EPSILON: f64 = -36.043_653_389_117_15; +const LN_MIN_POSITIVE: f64 = -708.396_418_532_264_1; +const SQRT_PI: f64 = 1.772_453_850_905_516; +const GAUSS_LEGENDRE_Y: [f64; 64] = [ + 0.000_347_479_132_113_914_8, + 0.001_829_941_614_022_39, + 0.004_493_314_261_627_856_5, + 0.008_331_873_057_687_012, + 0.013_336_586_105_044_512, + 0.019_495_600_173_973_116, + 0.026_794_312_570_798_617, + 0.035_215_413_934_030_215, + 0.044_738_931_460_748_59, + 0.055_342_277_002_442_93, + 0.067_000_300_922_953_56, + 0.079_685_351_873_709_79, + 0.093_367_342_438_601_23, + 0.108_013_820_528_329_31, + 0.123_590_046_369_734_03, + 0.140_059_074_914_194_56, + 0.157_381_843_472_883_36, + 0.175_517_264_372_671_34, + 0.194_422_322_413_803_36, + 0.214_052_176_898_683, + 0.234_360_267_990_052_72, + 0.255_298_427_146_473_55, + 0.276_816_991_373_268, + 0.298_864_921_018_004_2, + 0.321_389_920_831_165_9, + 0.344_338_564_004_894_5, + 0.367_656_418_895_616_3, + 0.391_288_178_129_996_44, + 0.415_177_789_788_003_57, + 0.439_268_590_351_939_7, + 0.463_503_439_106_100_5, + 0.487_824_853_668_287_76, + 0.512_175_146_331_712_2, + 0.536_496_560_893_899_5, + 0.560_731_409_648_060_2, + 0.584_822_210_211_996_4, + 0.608_711_821_870_003_5, + 0.632_343_581_104_383_7, + 0.655_661_435_995_105_5, + 0.678_610_079_168_834_1, + 0.701_135_078_981_995_8, + 0.723_183_008_626_732, + 0.744_701_572_853_526_5, + 0.765_639_732_009_947_3, + 0.785_947_823_101_317_1, + 0.805_577_677_586_196_7, + 0.824_482_735_627_328_7, + 0.842_618_156_527_116_7, + 0.859_940_925_085_805_4, + 0.876_409_953_630_266, + 0.891_986_179_471_670_7, + 0.906_632_657_561_398_8, + 0.920_314_648_126_290_3, + 0.932_999_699_077_046_4, + 0.944_657_722_997_557, + 0.955_261_068_539_251_4, + 0.964_784_586_065_969_8, + 0.973_205_687_429_201_4, + 0.980_504_399_826_026_9, + 0.986_663_413_894_955_5, + 0.991_668_126_942_313, + 0.995_506_685_738_372_1, + 0.998_170_058_385_977_6, + 0.999_652_520_867_886_1, +]; +const GAUSS_LEGENDRE_W: [f64; 64] = [ + 0.000_891_640_360_848_133_9, + 0.002_073_516_630_281_260_4, + 0.003_252_228_984_489_184, + 0.004_423_379_913_181_967, + 0.005_584_069_730_065_538, + 0.006_731_523_948_359_368, + 0.007_863_015_238_012_236, + 0.008_975_857_887_848_613, + 0.010_067_411_576_765_09, + 0.011_135_086_904_191_606, + 0.012_176_351_284_355_466, + 0.013_188_734_857_527_322, + 0.014_169_836_307_129_73, + 0.015_117_328_536_201_213, + 0.016_028_964_177_425_803, + 0.016_902_580_918_570_807, + 0.017_736_106_628_441_18, + 0.018_527_564_270_120_034, + 0.019_275_076_589_307_8, + 0.019_976_870_566_360_18, + 0.020_631_281_621_311_774, + 0.021_236_757_561_826_813, + 0.021_791_862_264_661_736, + 0.022_295_279_081_878_29, + 0.022_745_813_963_709_06, + 0.023_142_398_290_657_198, + 0.023_484_091_408_104_996, + 0.023_770_082_857_415_158, + 0.023_999_694_298_229_166, + 0.024_172_381_117_401_45, + 0.024_287_733_720_751_697, + 0.024_345_478_504_569_865, + 0.024_345_478_504_569_865, + 0.024_287_733_720_751_697, + 0.024_172_381_117_401_45, + 0.023_999_694_298_229_166, + 0.023_770_082_857_415_158, + 0.023_484_091_408_104_996, + 0.023_142_398_290_657_198, + 0.022_745_813_963_709_06, + 0.022_295_279_081_878_29, + 0.021_791_862_264_661_736, + 0.021_236_757_561_826_813, + 0.020_631_281_621_311_774, + 0.019_976_870_566_360_18, + 0.019_275_076_589_307_8, + 0.018_527_564_270_120_034, + 0.017_736_106_628_441_18, + 0.016_902_580_918_570_807, + 0.016_028_964_177_425_803, + 0.015_117_328_536_201_213, + 0.014_169_836_307_129_73, + 0.013_188_734_857_527_322, + 0.012_176_351_284_355_466, + 0.011_135_086_904_191_606, + 0.010_067_411_576_765_09, + 0.008_975_857_887_848_613, + 0.007_863_015_238_012_236, + 0.006_731_523_948_359_368, + 0.005_584_069_730_065_538, + 0.004_423_379_913_181_967, + 0.003_252_228_984_489_184, + 0.002_073_516_630_281_260_4, + 0.000_891_640_360_848_133_9, +]; + +/// Error type for special-function routines. +#[derive(Debug, Clone, PartialEq)] +pub enum SpecialError { + /// Input was outside the mathematical domain of the requested function. + InvalidInput { + /// Explanation of the invalid input. + message: String, + }, + /// A fixed iteration budget was exhausted before convergence. + MaxIterations { + /// Name of the routine that exhausted its budget. + function: &'static str, + /// Number of iterations attempted. + iterations: usize, + }, + /// A non-finite or otherwise unusable intermediate value was encountered. + NumericalIssue { + /// Explanation of the numerical issue. + message: String, + }, +} + +impl fmt::Display for SpecialError { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + match self { + Self::InvalidInput { message } => write!(f, "Invalid input: {message}"), + Self::MaxIterations { + function, + iterations, + } => write!( + f, + "Maximum iterations ({iterations}) exceeded in {function}" + ), + Self::NumericalIssue { message } => write!(f, "Numerical issue: {message}"), + } + } +} + +impl std::error::Error for SpecialError {} + +/// Options for the bracket-preserving inverse regularized-beta solver. +#[derive(Debug, Clone, Copy)] +pub struct InverseBetaOptions { + /// Maximum number of Newton/Halley iterations. + pub max_newton: usize, + /// Maximum number of bisection-only iterations after Newton is exhausted. + pub max_bisect: usize, + /// Absolute probability-scale residual tolerance. + pub probability_tolerance: f64, + /// Relative-x step tolerance. + pub relative_x_tolerance: f64, +} + +impl Default for InverseBetaOptions { + fn default() -> Self { + Self { + max_newton: 100, + max_bisect: 200, + probability_tolerance: 1.0e-15, + relative_x_tolerance: 1.0e-12, + } + } +} + +/// Inverse-beta quantile and its complement. +/// +/// For roots at or below `0.5`, [`invbetai`] solves directly for [`Self::x`] +/// and derives [`Self::complement`]. For roots above `0.5`, it solves the +/// swapped complement problem directly for [`Self::complement`] and derives +/// [`Self::x`]. This reconciles the design-note requirement to avoid losing the +/// well-conditioned side to `1 - x` cancellation while preserving the +/// user-facing `x`/`complement` pair. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct BetaQuantile { + /// The quantile `x` such that `I_x(a,b) = p`. + pub x: f64, + /// The complement `1 - x`, returned with the quantile for tail-sensitive callers. + pub complement: f64, +} + +/// Lower/upper tail probabilities for a regularized distribution function. +/// +/// The pair is used for both beta `I_x(a,b)`/`1 - I_x(a,b)` and gamma +/// `P(a,x)`/`Q(a,x)` helper paths. +#[derive(Debug, Clone, Copy)] +struct TailPair { + lower: f64, + upper: f64, +} + +/// Natural logarithm of the gamma function. +/// +/// This uses `libm::lgamma` explicitly, matching the Jeffreys interval +/// determinism contract for special-function paths. The v3 design note's +/// implementation plan named a Lanczos approximation; delegating to +/// `libm::lgamma` is an intentional reconciliation with the same note's +/// mandatory libm-by-name determinism contract. /// -/// Drop-in replacement for `scipy.special.gammaln` on the positive real -/// axis. Uses the 14-term Lanczos approximation (Numerical Recipes 3rd -/// ed., section 6.1), accurate to roughly machine precision. +/// # Errors /// -/// Returns NaN for `x <= 0`. +/// Returns an error if `x` is not positive and finite, or if `libm::lgamma` +/// returns a non-finite value. /// /// # Examples /// /// ``` /// use pecos_num::special::ln_gamma; /// -/// // Gamma(1) = 1, Gamma(2) = 1 -/// assert!(ln_gamma(1.0).abs() < 1e-14); -/// assert!(ln_gamma(2.0).abs() < 1e-14); -/// // Gamma(0.5) = sqrt(pi) -/// let half = std::f64::consts::PI.sqrt().ln(); -/// assert!((ln_gamma(0.5) - half).abs() < 1e-14); +/// let value = ln_gamma(5.0).unwrap(); +/// assert!((value - 3.178_053_830_347_945_8).abs() < 1e-14); +/// ``` +pub fn ln_gamma(x: f64) -> Result { + ensure_positive("x", x)?; + let value = libm::lgamma(x); + if value.is_finite() { + Ok(value) + } else { + Err(SpecialError::NumericalIssue { + message: "libm::lgamma returned a non-finite value".to_string(), + }) + } +} + +/// Regularized incomplete beta function `I_x(a,b)`. +/// +/// The implementation evaluates the continued fraction in the better +/// conditioned tail and uses the reflection identity +/// `I_x(a,b) = 1 - I_{1-x}(b,a)` when the complement is the stable branch. +/// +/// # Errors +/// +/// Returns an error if `a <= 0`, `b <= 0`, `x` is outside `[0, 1]`, a +/// non-finite value is encountered, or the continued-fraction iteration budget +/// is exhausted. +/// +/// # Examples +/// +/// ``` +/// use pecos_num::special::betai; +/// +/// let value = betai(2.0, 2.0, 0.5).unwrap(); +/// assert!((value - 0.5).abs() < 1e-14); +/// ``` +pub fn betai(a: f64, b: f64, x: f64) -> Result { + Ok(betai_pair(a, b, x)?.lower) +} + +/// Inverse regularized incomplete beta function. +/// +/// Returns both `x` and `1 - x` for the solution to `I_x(a,b) = p`. When the +/// root lies above `0.5`, this solves the swapped problem for `1 - x` directly +/// and derives `x`; otherwise it solves `x` directly and derives the +/// complement. For `p = 0.5` with both shapes at least `1_000_000`, it uses the +/// Peizer-Pratt/Kerman large-shape beta-median asymptotic because the GL +/// forward CDF noise dominates further Newton refinement. See [`BetaQuantile`] +/// for the exact-side contract. +/// +/// # Errors +/// +/// Returns an error if the inputs are outside the function domain, a +/// non-finite value is encountered, or the fixed iteration budget is exhausted. +/// +/// # Examples +/// +/// ``` +/// use pecos_num::special::invbetai; +/// +/// let quantile = invbetai(0.5, 2.0, 2.0).unwrap(); +/// assert!((quantile.x - 0.5).abs() < 1e-14); /// ``` -#[must_use] -pub fn ln_gamma(x: f64) -> f64 { - const COF: [f64; 14] = [ - 57.156_235_665_862_92, - -59.597_960_355_475_49, - 14.136_097_974_741_746, - -0.491_913_816_097_620_2, - 3.399_464_998_481_189e-5, - 4.652_362_892_704_858e-5, - -9.837_447_530_487_956e-5, - 1.580_887_032_249_125e-4, - -2.102_644_417_241_049e-4, - 2.174_396_181_152_126e-4, - -1.643_181_065_367_639e-4, - 8.441_822_398_385_274e-5, - -2.619_083_840_158_141e-5, - 3.689_918_265_953_162e-6, - ]; - const LANCZOS_G: f64 = 5.242_187_5; - const SQRT_2PI: f64 = 2.506_628_274_631_000_5; +pub fn invbetai(p: f64, a: f64, b: f64) -> Result { + invbetai_with_options(p, a, b, InverseBetaOptions::default()) +} + +/// Inverse regularized incomplete beta function with explicit iteration options. +/// +/// This is primarily useful for tests and diagnostics that need to assert the +/// fixed-budget error behavior. +/// +/// # Errors +/// +/// Returns an error if the inputs or options are invalid, a non-finite value is +/// encountered, or the configured iteration budget is exhausted. +pub fn invbetai_with_options( + p: f64, + a: f64, + b: f64, + options: InverseBetaOptions, +) -> Result { + ensure_unit_interval("p", p)?; + ensure_positive("a", a)?; + ensure_positive("b", b)?; + ensure_positive("probability_tolerance", options.probability_tolerance)?; + ensure_positive("relative_x_tolerance", options.relative_x_tolerance)?; + + if p <= 0.0 { + return Ok(BetaQuantile { + x: 0.0, + complement: 1.0, + }); + } + if p >= 1.0 { + return Ok(BetaQuantile { + x: 1.0, + complement: 0.0, + }); + } + if let Some(quantile) = large_shape_median_quantile(p, a, b) { + return Ok(quantile); + } + + let half = betai_pair(a, b, 0.5)?; + if p > half.lower { + let complement = invbetai_lower_tail(1.0 - p, b, a, options)?; + return Ok(BetaQuantile { + x: complement.complement, + complement: complement.x, + }); + } + + invbetai_lower_tail(p, a, b, options) +} + +fn invbetai_lower_tail( + p: f64, + a: f64, + b: f64, + options: InverseBetaOptions, +) -> Result { + let total_iterations = options.max_newton + options.max_bisect; + if total_iterations == 0 { + return Err(SpecialError::MaxIterations { + function: "invbetai", + iterations: 0, + }); + } + + let mut lo = 0.0; + let mut hi = 1.0; + let mut x = initial_invbetai_guess(p, a, b); + if !is_strictly_inside(x, lo, hi) { + x = midpoint(lo, hi); + } + + let mut previous_x = x; + let mut have_previous_step = false; + let mut lo_established = false; + let mut hi_established = false; + let effective_probability_tolerance = effective_probability_tolerance(options, a, b); + + for iteration in 0..total_iterations { + let residual = invbetai_residual(p, a, b, x)?; + let abs_residual = residual.abs(); + + if same_float(abs_residual, 0.0) { + return Ok(direct_quantile(x)); + } + + if residual < 0.0 { + lo = x; + lo_established = true; + } else { + hi = x; + hi_established = true; + } + + let bracket_established = lo_established && hi_established; + let step_ok = have_previous_step + && (x - previous_x).abs() <= options.relative_x_tolerance * previous_x.abs(); + let bracket_limited = bracket_established && bracket_is_representation_limited(lo, hi, x); + if abs_residual <= effective_probability_tolerance && (step_ok || bracket_limited) { + return Ok(direct_quantile(x)); + } + if bracket_limited { + return Ok(direct_quantile(midpoint(lo, hi))); + } + + let newton = (iteration < options.max_newton) + .then(|| newton_candidate(a, b, x, residual)) + .flatten() + .filter(|&candidate| is_strictly_inside(candidate, lo, hi)); + let mut candidate = + newton.unwrap_or_else(|| fallback_candidate(residual, x, lo, hi, bracket_established)); + + if same_float(candidate, x) { + if abs_residual <= effective_probability_tolerance { + return Ok(direct_quantile(x)); + } + candidate = fallback_candidate(residual, x, lo, hi, bracket_established); + } + + if same_float(candidate, x) && bracket_established { + let bisected = midpoint(lo, hi); + if same_float(bisected, x) { + if bracket_limited { + return Ok(direct_quantile(bisected)); + } + return Err(SpecialError::MaxIterations { + function: "invbetai", + iterations: iteration + 1, + }); + } + candidate = bisected; + } + + if same_float(candidate, x) { + return Err(SpecialError::MaxIterations { + function: "invbetai", + iterations: iteration + 1, + }); + } + + previous_x = x; + x = candidate; + have_previous_step = true; + } + + Err(SpecialError::MaxIterations { + function: "invbetai", + iterations: total_iterations, + }) +} + +fn effective_probability_tolerance(options: InverseBetaOptions, a: f64, b: f64) -> f64 { + if a > BETA_GAUSS_LEGENDRE_SWITCH && b > BETA_GAUSS_LEGENDRE_SWITCH { + // This is a local solver tolerance, not a global bound on GL forward + // error against scipy. Large symmetric CDF probes can differ by up to + // roughly `(a + b) * eps` from log-difference rounding; measured probes + // were about `0.26 * (a + b) * eps`. Newton only needs a deterministic + // zero of this fixed forward function, and the bracket-collapse check + // plus huge local density keep the returned x stable. + options + .probability_tolerance + .max(BETA_GAUSS_LEGENDRE_PROBABILITY_FLOOR) + } else { + options.probability_tolerance + } +} + +fn fallback_candidate(residual: f64, x: f64, lo: f64, hi: f64, bracket_established: bool) -> f64 { + if bracket_established { + return midpoint(lo, hi); + } + + if residual < 0.0 { + let expanded = (2.0 * x).min(0.5).max(midpoint(x, 0.5)); + if is_strictly_inside(expanded, x, hi) { + expanded + } else { + midpoint(x, hi) + } + } else { + let contracted = 0.5 * x; + if is_strictly_inside(contracted, lo, x) { + contracted + } else { + midpoint(lo, x) + } + } +} + +fn large_shape_median_quantile(p: f64, a: f64, b: f64) -> Option { + if !same_float(p, 0.5) || a < BETA_MEDIAN_ASYMPTOTIC_MIN || b < BETA_MEDIAN_ASYMPTOTIC_MIN { + return None; + } + + // Peizer-Pratt/Kerman large-shape beta-median asymptotic. In this regime + // the GL forward CDF noise dominates further Newton refinement; the + // approximation is the Boost/Temme-style asymptotic inverse used for the + // median-only path. + let denominator = a + b - 2.0 / 3.0; + if a > b { + let complement = (b - 1.0 / 3.0) / denominator; + Some(BetaQuantile { + x: 1.0 - complement, + complement, + }) + } else { + let x = (a - 1.0 / 3.0) / denominator; + Some(BetaQuantile { + x, + complement: 1.0 - x, + }) + } +} + +#[inline] +fn same_float(left: f64, right: f64) -> bool { + left.to_bits() == right.to_bits() +} + +fn betai_pair(a: f64, b: f64, x: f64) -> Result { + ensure_positive("a", a)?; + ensure_positive("b", b)?; + ensure_unit_interval("x", x)?; + + if x <= 0.0 { + return Ok(TailPair { + lower: 0.0, + upper: 1.0, + }); + } + if x >= 1.0 { + return Ok(TailPair { + lower: 1.0, + upper: 0.0, + }); + } + + if let Some(cdf) = beta_large_shape_pair(a, b, x)? { + return Ok(cdf); + } + if a > BETA_GAUSS_LEGENDRE_SWITCH && b > BETA_GAUSS_LEGENDRE_SWITCH { + return betai_gauss_legendre_pair(a, b, x); + } + + let prefactor = beta_prefactor(a, b, x)?; + let threshold = (a + 1.0) / (a + b + 2.0); + if x <= threshold || (x <= 0.1 && b * x <= 5.0) { + let lower = lower_beta_tail(a, b, x, prefactor)?; + ensure_finite("regularized beta lower tail", lower)?; + Ok(TailPair { + lower, + upper: 1.0 - lower, + }) + } else { + let upper = lower_beta_tail(b, a, 1.0 - x, prefactor)?; + ensure_finite("regularized beta upper tail", upper)?; + Ok(TailPair { + lower: 1.0 - upper, + upper, + }) + } +} + +fn lower_beta_tail(a: f64, b: f64, x: f64, prefactor: f64) -> Result { + if b * x <= 1.0 && x <= 0.95 { + beta_power_series(a, b, x) + } else { + Ok(prefactor * betacf(a, b, x)? / a) + } +} + +fn beta_prefactor(a: f64, b: f64, x: f64) -> Result { + let value = libm::exp(log_beta_normalizer(a, b)? + a * libm::log(x) + b * libm::log1p(-x)); + ensure_finite("regularized beta prefactor", value)?; + Ok(value) +} + +fn beta_power_series(a: f64, b: f64, x: f64) -> Result { + let reciprocal_a = 1.0 / a; + let mut term = (1.0 - b) * x; + let mut value = term / (a + 1.0); + let mut sum = value; + let tolerance = BETA_POWER_SERIES_EPSILON * reciprocal_a; + let mut n = 2.0; + + for _iteration in 2..=BETA_POWER_SERIES_MAX_ITERATIONS { + term *= (n - b) * x / n; + value = term / (a + n); + sum += value; + if value.abs() <= tolerance { + let log_prefactor = log_beta_normalizer(a, b)? + a * libm::log(x); + let result = (sum + reciprocal_a) * libm::exp(log_prefactor); + ensure_finite("regularized beta power series", result)?; + return Ok(result); + } + n += 1.0; + } + + Err(SpecialError::MaxIterations { + function: "beta_power_series", + iterations: BETA_POWER_SERIES_MAX_ITERATIONS, + }) +} + +fn betai_gauss_legendre_pair(a: f64, b: f64, x: f64) -> Result { + let a1 = a - 1.0; + let b1 = b - 1.0; + let mu = a / (a + b); + let ln_mu = libm::log(mu); + let ln_mu_complement = libm::log1p(-mu); + let sigma = libm::sqrt(a * b / ((a + b) * (a + b) * (a + b + 1.0))); + let integrate_upper = x > mu; + + let xu = if integrate_upper { + if x >= 1.0 { + return Ok(TailPair { + lower: 1.0, + upper: 0.0, + }); + } + (mu + 10.0 * sigma).max(x + 5.0 * sigma).min(1.0) + } else { + if x <= 0.0 { + return Ok(TailPair { + lower: 0.0, + upper: 1.0, + }); + } + (mu - 10.0 * sigma).min(x - 5.0 * sigma).max(0.0) + }; + + let span = if integrate_upper { xu - x } else { x - xu }; + let panel_width = span / usize_to_f64(BETA_GAUSS_LEGENDRE_PANELS); + let start = if integrate_upper { x } else { xu }; + let mut sum = 0.0; + for panel in 0..BETA_GAUSS_LEGENDRE_PANELS { + let panel_start = start + usize_to_f64(panel) * panel_width; + for (&y, &weight) in GAUSS_LEGENDRE_Y.iter().zip(GAUSS_LEGENDRE_W.iter()) { + let t = panel_start + panel_width * y; + sum += weight + * libm::exp( + a1 * (libm::log(t) - ln_mu) + b1 * (libm::log1p(-t) - ln_mu_complement), + ); + } + } + + let tail = sum + * panel_width + * libm::exp(a1 * ln_mu + b1 * ln_mu_complement + log_beta_normalizer(a, b)?); + ensure_finite("large-shape beta quadrature", tail)?; + if integrate_upper { + Ok(TailPair { + lower: 1.0 - tail, + upper: tail, + }) + } else { + Ok(TailPair { + lower: tail, + upper: 1.0 - tail, + }) + } +} + +fn beta_large_shape_pair(a: f64, b: f64, x: f64) -> Result, SpecialError> { + if let Some(pair) = beta_large_b_saturation_pair(a, b, x)? { + return Ok(Some(pair)); + } + + if beta_large_b_gate(a, b, x) { + return Ok(Some(beta_large_b_asymptotic_pair(a, b, x)?)); + } + + let x_complement = 1.0 - x; + if let Some(swapped) = beta_large_b_saturation_pair(b, a, x_complement)? { + return Ok(Some(TailPair { + lower: swapped.upper, + upper: swapped.lower, + })); + } + + if beta_large_b_gate(b, a, x_complement) { + let swapped = beta_large_b_asymptotic_pair(b, a, x_complement)?; + return Ok(Some(TailPair { + lower: swapped.upper, + upper: swapped.lower, + })); + } + + Ok(None) +} + +fn beta_large_b_gate(a: f64, b: f64, x: f64) -> bool { + b >= BETA_LARGE_SHAPE_MIN + && a <= BETA_LARGE_SHAPE_MODERATE_MAX + && x <= BETA_LARGE_SHAPE_X_MAX + && a * a / b <= BETA_LARGE_SHAPE_EXPANSION_PARAMETER_MAX +} + +fn beta_large_b_saturation_pair(a: f64, b: f64, x: f64) -> Result, SpecialError> { + if b < BETA_LARGE_SHAPE_MIN + || a > BETA_LARGE_SHAPE_MODERATE_MAX + || a * a / b > BETA_LARGE_SHAPE_EXPANSION_PARAMETER_MAX + { + return Ok(None); + } + + let w = -b * libm::log1p(-x); + if w <= a { + return Ok(None); + } + + if gamma_upper_is_below_precision(a, w)? { + Ok(Some(TailPair { + lower: 1.0, + upper: 0.0, + })) + } else { + Ok(None) + } +} + +// DiDonato and Morris, ACM TOMS 18(3), 1992, Algorithm 708, uses the BGRAT +// large-b expansion to reduce the beta tail to incomplete-gamma tails. With +// `w = -b log1p(-x)`, the beta integrand is a gamma density multiplied by the +// fixed series `((1 - exp(-z)) / z)^(a - 1)`, `z = w / b`. The fixed expansion +// budget is enabled only for the validated Jeffreys regime `b >= 1e4`, +// `x <= 1e-2`, `a <= 3001`, and `a*a/b <= 0.1`. The scale bound, not the fixed +// floor, is the active lower limit for the Round-5 small-a pocket rows; far-tail +// probes with gamma-Q below machine precision saturate before reaching the +// reflected continued fraction. +fn beta_large_b_asymptotic_pair(a: f64, b: f64, x: f64) -> Result { + let w = -b * libm::log1p(-x); + if w <= 0.0 { + return Ok(TailPair { + lower: 0.0, + upper: 1.0, + }); + } + + let coefficients = large_b_expansion_coefficients(a); + let mut rising = 1.0; + let mut b_power = 1.0; + let mut lower_sum = 0.0; + let mut upper_sum = 0.0; + + for (index, coefficient) in coefficients.iter().enumerate() { + if index > 0 { + rising *= a + usize_to_f64(index - 1); + b_power *= b; + } + + let weight = coefficient * rising / b_power; + let gamma = regularized_gamma_pq(a + usize_to_f64(index), w)?; + lower_sum += weight * gamma.lower; + upper_sum += weight * gamma.upper; + } + + let normalizer = libm::exp(a * libm::log(b) - log_gamma_delta(b, a)); + ensure_finite("large-b beta expansion normalizer", normalizer)?; + if normalizer <= 0.0 { + return Err(SpecialError::NumericalIssue { + message: "large-b beta expansion normalizer was non-positive".to_string(), + }); + } + + let lower = clean_unit_probability("large-b beta lower tail", lower_sum / normalizer)?; + let upper = clean_unit_probability("large-b beta upper tail", upper_sum / normalizer)?; + Ok(TailPair { lower, upper }) +} + +fn large_b_expansion_coefficients(a: f64) -> [f64; BETA_LARGE_B_EXPANSION_LEN] { + let exponent = a - 1.0; + let mut base = [0.0; BETA_LARGE_B_EXPANSION_LEN]; + let mut factorial = 1.0; + for (index, coefficient) in base.iter_mut().enumerate().skip(1) { + factorial *= usize_to_f64(index + 1); + let sign = if index % 2 == 0 { 1.0 } else { -1.0 }; + *coefficient = sign / factorial; + } + + let mut coefficients = [0.0; BETA_LARGE_B_EXPANSION_LEN]; + coefficients[0] = 1.0; + let mut base_power = [0.0; BETA_LARGE_B_EXPANSION_LEN]; + base_power[0] = 1.0; + let mut binomial = 1.0; + + for order in 1..BETA_LARGE_B_EXPANSION_LEN { + base_power = multiply_series(&base_power, &base); + binomial *= (exponent - usize_to_f64(order - 1)) / usize_to_f64(order); + for index in order..BETA_LARGE_B_EXPANSION_LEN { + coefficients[index] += binomial * base_power[index]; + } + } + + coefficients +} + +fn multiply_series( + left: &[f64; BETA_LARGE_B_EXPANSION_LEN], + right: &[f64; BETA_LARGE_B_EXPANSION_LEN], +) -> [f64; BETA_LARGE_B_EXPANSION_LEN] { + let mut product = [0.0; BETA_LARGE_B_EXPANSION_LEN]; + for left_index in 0..BETA_LARGE_B_EXPANSION_LEN { + for right_index in 0..(BETA_LARGE_B_EXPANSION_LEN - left_index) { + product[left_index + right_index] += left[left_index] * right[right_index]; + } + } + product +} +fn regularized_gamma_pq(a: f64, x: f64) -> Result { + ensure_positive("a", a)?; + if x < 0.0 || !x.is_finite() { + return Err(SpecialError::InvalidInput { + message: "x must be non-negative and finite".to_string(), + }); + } if x <= 0.0 { - return f64::NAN; + return Ok(TailPair { + lower: 0.0, + upper: 1.0, + }); + } + + if let Some(pair) = regularized_gamma_half_integer_pq(a, x)? { + return Ok(pair); + } + + if x < a + 1.0 { + let lower = regularized_gamma_p_series(a, x)?; + Ok(TailPair { + lower, + upper: 1.0 - lower, + }) + } else { + if gamma_upper_is_below_precision(a, x)? { + return Ok(TailPair { + lower: 1.0, + upper: 0.0, + }); + } + + let upper = regularized_gamma_q_fraction(a, x)?; + Ok(TailPair { + lower: 1.0 - upper, + upper, + }) + } +} + +fn gamma_upper_is_below_precision(a: f64, x: f64) -> Result { + let log_scale = -x + a * libm::log(x) - ln_gamma(a)?; + if log_scale <= LN_MIN_POSITIVE { + return Ok(true); + } + + let denominator = (x + 1.0 - a).max(1.0); + Ok(log_scale - libm::log(denominator) <= LN_EPSILON) +} + +fn regularized_gamma_half_integer_pq(a: f64, x: f64) -> Result, SpecialError> { + if x < a || x >= GAMMA_HALF_INTEGER_MAX_X { + return Ok(None); } - let mut denom = x; - let tmp = x + LANCZOS_G; - let tmp = (x + 0.5) * tmp.ln() - tmp; - let mut series = 0.999_999_999_999_997_1; - for c in COF { - denom += 1.0; - series += c / denom; + let steps_f = a - 0.5; + if steps_f < 0.0 || !same_float(steps_f, steps_f.round()) { + return Ok(None); + } + + #[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)] + // Half-integer gate bounds steps to a non-negative exact integer + let steps = steps_f as usize; + if steps > GAMMA_HALF_INTEGER_MAX_STEPS { + return Ok(None); + } + + let root_x = libm::sqrt(x); + let mut upper = libm::erfc(root_x); + let mut term = 2.0 * root_x * libm::exp(-x) / SQRT_PI; + let mut shape = 0.5; + + for _ in 0..steps { + upper += term; + term *= x / (shape + 1.0); + shape += 1.0; + } + + let upper = clean_unit_probability("half-integer gamma upper tail", upper)?; + let lower = clean_unit_probability("half-integer gamma lower tail", 1.0 - upper)?; + Ok(Some(TailPair { lower, upper })) +} + +fn regularized_gamma_p_series(a: f64, x: f64) -> Result { + let mut ap = a; + let mut delta = 1.0 / a; + let mut sum = delta; + + for _iteration in 1..=GAMMA_MAX_ITERATIONS { + ap += 1.0; + delta *= x / ap; + sum += delta; + if delta.abs() <= sum.abs() * GAMMA_EPSILON { + let value = sum * libm::exp(-x + a * libm::log(x) - ln_gamma(a)?); + ensure_finite("regularized gamma lower series", value)?; + return Ok(value); + } + } + + Err(SpecialError::MaxIterations { + function: "regularized_gamma_p_series", + iterations: GAMMA_MAX_ITERATIONS, + }) +} + +fn regularized_gamma_q_fraction(a: f64, x: f64) -> Result { + let mut b = x + 1.0 - a; + let mut c = 1.0 / FPMIN; + let mut d = 1.0 / b; + let mut h = d; + + for iteration in 1..=GAMMA_MAX_ITERATIONS { + let i = usize_to_f64(iteration); + let an = -i * (i - a); + b += 2.0; + d = an * d + b; + if d.abs() < FPMIN { + d = FPMIN; + } + c = b + an / c; + if c.abs() < FPMIN { + c = FPMIN; + } + d = 1.0 / d; + let delta = d * c; + h *= delta; + if (delta - 1.0).abs() <= GAMMA_EPSILON { + let value = libm::exp(-x + a * libm::log(x) - ln_gamma(a)?) * h; + ensure_finite("regularized gamma upper fraction", value)?; + return Ok(value); + } } - tmp + (SQRT_2PI * series / x).ln() + + Err(SpecialError::MaxIterations { + function: "regularized_gamma_q_fraction", + iterations: GAMMA_MAX_ITERATIONS, + }) } -/// Continued-fraction evaluation for the incomplete beta function -/// (Numerical Recipes 3rd ed., section 6.4, via Lentz's method). -fn betacf(a: f64, b: f64, x: f64) -> f64 { +fn betacf(a: f64, b: f64, x: f64) -> Result { let qab = a + b; let qap = a + 1.0; let qam = a - 1.0; let mut c = 1.0; let mut d = 1.0 - qab * x / qap; - if d.abs() < BETACF_FPMIN { - d = BETACF_FPMIN; + if d.abs() < FPMIN { + d = FPMIN; } d = 1.0 / d; let mut h = d; - for m in 1..=BETACF_MAX_ITER { - let m = f64::from(m); + for iteration in 1..=BETACF_MAX_ITERATIONS { + let m = usize_to_f64(iteration); let m2 = 2.0 * m; - - // Even step of the continued fraction. - let aa = m * (b - m) * x / ((qam + m2) * (a + m2)); + let mut aa = m * (b - m) * x / ((qam + m2) * (a + m2)); d = 1.0 + aa * d; - if d.abs() < BETACF_FPMIN { - d = BETACF_FPMIN; + if d.abs() < FPMIN { + d = FPMIN; } c = 1.0 + aa / c; - if c.abs() < BETACF_FPMIN { - c = BETACF_FPMIN; + if c.abs() < FPMIN { + c = FPMIN; } d = 1.0 / d; h *= d * c; - // Odd step. - let aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2)); + aa = -((a + m) * (qab + m) * x) / ((a + m2) * (qap + m2)); d = 1.0 + aa * d; - if d.abs() < BETACF_FPMIN { - d = BETACF_FPMIN; + if d.abs() < FPMIN { + d = FPMIN; } c = 1.0 + aa / c; - if c.abs() < BETACF_FPMIN { - c = BETACF_FPMIN; + if c.abs() < FPMIN { + c = FPMIN; } d = 1.0 / d; - let del = d * c; - h *= del; + let delta = d * c; + h *= delta; - if (del - 1.0).abs() <= BETACF_EPS { - return h; + ensure_finite("regularized beta continued fraction", h)?; + if (delta - 1.0).abs() <= BETACF_EPSILON { + return Ok(h); } } - // The fraction converges for all valid inputs; reaching the iteration - // cap means the inputs were extreme enough that the result is unusable. - f64::NAN + + Err(SpecialError::MaxIterations { + function: "betacf", + iterations: BETACF_MAX_ITERATIONS, + }) } -/// Regularized incomplete beta function `I_x(a, b)` for `a, b > 0` and -/// `x` in `[0, 1]`. -/// -/// Drop-in replacement for `scipy.special.betainc`. This is also the CDF -/// of the Beta(a, b) distribution evaluated at `x`. -/// -/// Returns NaN outside the valid domain. -/// -/// # Examples -/// -/// ``` -/// use pecos_num::special::betainc_reg; -/// -/// // Symmetric case: I_{0.5}(a, a) = 0.5 -/// assert!((betainc_reg(10.0, 10.0, 0.5) - 0.5).abs() < 1e-12); -/// // Boundaries -/// assert_eq!(betainc_reg(2.0, 3.0, 0.0), 0.0); -/// assert_eq!(betainc_reg(2.0, 3.0, 1.0), 1.0); -/// ``` -#[must_use] -pub fn betainc_reg(a: f64, b: f64, x: f64) -> f64 { - let domain_ok = a > 0.0 && b > 0.0 && (0.0..=1.0).contains(&x); - if !domain_ok { - return f64::NAN; - } - if x <= 0.0 { - return 0.0; +fn invbetai_residual(p: f64, a: f64, b: f64, x: f64) -> Result { + let cdf = betai_pair(a, b, x)?; + if p <= 0.5 { + Ok(cdf.lower - p) + } else { + Ok((1.0 - p) - cdf.upper) } - if x >= 1.0 { - return 1.0; +} + +fn newton_candidate(a: f64, b: f64, x: f64, residual: f64) -> Option { + let density = beta_density(a, b, x).ok()?; + if density <= 0.0 { + return None; } - // Prefactor x^a (1-x)^b / (a B(a,b)), computed in log space. - let ln_front = ln_gamma(a + b) - ln_gamma(a) - ln_gamma(b) + a * x.ln() + b * (1.0 - x).ln(); - let front = ln_front.exp(); + let first_order_step = residual / density; + let slope = (a - 1.0) / x - (b - 1.0) / (1.0 - x); + let denominator = 1.0 - 0.5 * first_order_step * slope; + let step = if denominator.is_finite() && denominator > 0.0 { + first_order_step / denominator + } else { + first_order_step + }; + let candidate = x - step; + candidate.is_finite().then_some(candidate) +} + +fn beta_density(a: f64, b: f64, x: f64) -> Result { + let log_density = + (a - 1.0) * libm::log(x) + (b - 1.0) * libm::log1p(-x) + log_beta_normalizer(a, b)?; + let density = libm::exp(log_density); + ensure_finite("beta density", density)?; + Ok(density) +} - // Use the continued fraction directly where it converges fastest, - // and the symmetry I_x(a,b) = 1 - I_{1-x}(b,a) otherwise. - if x < (a + 1.0) / (a + b + 2.0) { - front * betacf(a, b, x) / a +fn log_beta_normalizer(a: f64, b: f64) -> Result { + if a >= BETA_ASYMPTOTIC_NORMALIZER_MIN && b >= BETA_ASYMPTOTIC_NORMALIZER_MIN { + Ok(log_beta_normalizer_asymptotic(a, b)) + } else if a <= b && b >= 8.0 { + Ok(log_gamma_delta(b, a) - log_gamma_for_normalizer(a)?) + } else if a >= 8.0 { + Ok(log_gamma_delta(a, b) - log_gamma_for_normalizer(b)?) } else { - 1.0 - front * betacf(b, a, 1.0 - x) / b + Ok(ln_gamma(a + b)? - ln_gamma(a)? - ln_gamma(b)?) } } -/// Inverse of the regularized incomplete beta function: returns `x` such -/// that `betainc_reg(a, b, x) == p`. -/// -/// Matches `scipy.special.betaincinv` (this is the quantile / inverse CDF -/// of the Beta(a, b) distribution) over the parameter scales validated by -/// the test module: shape parameters up to roughly binomial-trial counts -/// of 1e12 and `p` away from the extreme tails by more than ~1e-15. The -/// upper tail is computed through the symmetry -/// `betainc_inv(a, b, p) = 1 - betainc_inv(b, a, 1 - p)` so both tails -/// share the well-conditioned lower-tail path; beyond those scales the -/// continued fraction can converge spuriously, so callers with extreme -/// parameters must validate independently. -/// -/// Follows Numerical Recipes 3rd ed., section 6.4: an initial guess from -/// Abramowitz & Stegun 26.5.22 refined by Halley iterations. -/// -/// Returns NaN outside the valid domain (`a, b > 0`, `p` in `[0, 1]`). -/// -/// # Examples -/// -/// ``` -/// use pecos_num::special::{betainc_inv, betainc_reg}; -/// -/// let x = betainc_inv(2.0, 3.0, 0.6); -/// assert!((betainc_reg(2.0, 3.0, x) - 0.6).abs() < 1e-10); -/// ``` -#[must_use] -pub fn betainc_inv(a: f64, b: f64, p: f64) -> f64 { - let domain_ok = a > 0.0 && b > 0.0 && (0.0..=1.0).contains(&p); - if !domain_ok { - return f64::NAN; - } - if p <= 0.0 { - return 0.0; +fn log_beta_normalizer_asymptotic(a: f64, b: f64) -> f64 { + let sum = a + b; + a * libm::log1p(b / a) + + b * libm::log1p(a / b) + + 0.5 * (libm::log(a) + libm::log(b) - libm::log(sum)) + - HALF_LN_TWO_PI + + stirling_correction(sum) + - stirling_correction(a) + - stirling_correction(b) +} + +fn log_gamma_for_normalizer(x: f64) -> Result { + if x >= 8.0 { + Ok(log_gamma_stirling(x)) + } else { + ln_gamma(x) } - if p >= 1.0 { - return 1.0; +} + +fn log_gamma_stirling(x: f64) -> f64 { + (x - 0.5) * libm::log(x) - x + HALF_LN_TWO_PI + stirling_correction(x) +} + +fn log_gamma_delta(base: f64, increment: f64) -> f64 { + let ratio = increment / base; + increment * libm::log(base) + + (increment - 0.5) * libm::log1p(ratio) + + base * log1pmx(ratio) + + stirling_correction(base + increment) + - stirling_correction(base) +} + +fn log1pmx(x: f64) -> f64 { + if x.abs() >= 0.1 { + return libm::log1p(x) - x; } - // Compute upper-tail quantiles through the lower tail of the mirrored - // distribution: `err = betainc_reg(...) - p` loses all precision when - // p is within ~1e-10 of 1 (the Halley correction then stalls on a - // cancelled residual), while 1 - p is exact in the mirrored call. - if p > 0.5 { - return 1.0 - betainc_inv(b, a, 1.0 - p); + + let mut term = x * x; + let mut sum = -0.5 * term; + for n in 3..=40 { + term *= x; + let signed = if n % 2 == 0 { -term } else { term }; + sum += signed / usize_to_f64(n); } + sum +} - let a1 = a - 1.0; - let b1 = b - 1.0; +fn stirling_correction(x: f64) -> f64 { + let inverse = 1.0 / x; + let inverse2 = inverse * inverse; + let inverse3 = inverse * inverse2; + let inverse5 = inverse3 * inverse2; + let inverse7 = inverse5 * inverse2; + let inverse9 = inverse7 * inverse2; + let inverse11 = inverse9 * inverse2; + let inverse13 = inverse11 * inverse2; + + inverse / 12.0 - inverse3 / 360.0 + inverse5 / 1_260.0 - inverse7 / 1_680.0 + inverse9 / 1_188.0 + - 691.0 * inverse11 / 360_360.0 + + inverse13 / 156.0 +} - let mut x: f64; +fn initial_invbetai_guess(p: f64, a: f64, b: f64) -> f64 { if a >= 1.0 && b >= 1.0 { - // Abramowitz & Stegun 26.5.22 via the normal quantile - // approximation 26.2.23. let pp = if p < 0.5 { p } else { 1.0 - p }; - let t = (-2.0 * pp.ln()).sqrt(); - let mut gauss = (2.30753 + t * 0.27061) / (1.0 + t * (0.99229 + t * 0.04481)) - t; + let t = libm::sqrt(-2.0 * libm::log(pp)); + let mut x = (2.307_53 + t * 0.270_61) / (1.0 + t * (0.992_29 + t * 0.044_81)) - t; if p < 0.5 { - gauss = -gauss; + x = -x; } - let al = (gauss * gauss - 3.0) / 6.0; + let al = (x * x - 3.0) / 6.0; let h = 2.0 / (1.0 / (2.0 * a - 1.0) + 1.0 / (2.0 * b - 1.0)); - let w = gauss * (al + h).sqrt() / h + let w = x * libm::sqrt(al + h) / h - (1.0 / (2.0 * b - 1.0) - 1.0 / (2.0 * a - 1.0)) * (al + 5.0 / 6.0 - 2.0 / (3.0 * h)); - x = a / (a + b * (2.0 * w).exp()); + interiorize(a / (a + b * libm::exp(2.0 * w))) } else { - let lna = (a / (a + b)).ln(); - let lnb = (b / (a + b)).ln(); - let t = (a * lna).exp() / a; - let u = (b * lnb).exp() / b; + let lna = libm::log(a / (a + b)); + let lnb = libm::log(b / (a + b)); + let t = libm::exp(a * lna) / a; + let u = libm::exp(b * lnb) / b; let w = t + u; - x = if p < t / w { - (a * w * p).powf(1.0 / a) + let x = if p < t / w { + libm::exp(libm::log(a * w * p) / a) } else { - 1.0 - (b * w * (1.0 - p)).powf(1.0 / b) + 1.0 - libm::exp(libm::log(b * w * (1.0 - p)) / b) }; + interiorize(x) } +} - let afac = ln_gamma(a + b) - ln_gamma(a) - ln_gamma(b); - for iteration in 0..10 { - if x <= 0.0 { - return 0.0; - } - if x >= 1.0 { - return 1.0; - } - let err = betainc_reg(a, b, x) - p; - let t = (a1 * x.ln() + b1 * (1.0 - x).ln() + afac).exp(); - let u = err / t; - // Halley step. - let step = u / (1.0 - 0.5 * f64::min(1.0, u * (a1 / x - b1 / (1.0 - x)))); - x -= step; - if x <= 0.0 { - x = 0.5 * (x + step); - } - if x >= 1.0 { - x = 0.5 * (x + step + 1.0); - } - if step.abs() < BETAINC_INV_EPS * x && iteration > 0 { - break; - } +fn ensure_positive(name: &'static str, value: f64) -> Result<(), SpecialError> { + if value.is_finite() && value > 0.0 { + Ok(()) + } else { + Err(SpecialError::InvalidInput { + message: format!("{name} must be positive and finite"), + }) } - x } -#[cfg(test)] -mod tests { - use super::*; - - /// Assert `actual` matches `expected` to a RELATIVE tolerance. - /// - /// A pure relative check: callers that legitimately expect a value of - /// exactly zero must handle that case separately (the `ln_gamma` test - /// does). An earlier version OR'd in an absolute `abs_err <= rel_tol` - /// fallback, which silently weakened the relative tolerance to an - /// absolute one for `|expected| < 1` (e.g. a claimed 1e-12 relative - /// became 1e-12 absolute — ~1000x looser at `expected = 1e-3`). - fn assert_close(actual: f64, expected: f64, rel_tol: f64) { - let denom = expected.abs().max(f64::MIN_POSITIVE); - let rel_err = (actual - expected).abs() / denom; - assert!( - rel_err <= rel_tol, - "expected {expected:.17e}, got {actual:.17e} (relative error {rel_err:.3e} > {rel_tol:.3e})" - ); - } - - // Reference values generated with: - // uv run python -c "from scipy import special; print(special.gammaln(x))" - // (scipy 1.x, double precision) - #[test] - fn ln_gamma_matches_scipy() { - let cases = [ - (0.5, 0.572_364_942_924_7), - (1.0, 0.0), - (1.5, -0.120_782_237_635_245_26), - (2.0, 0.0), - (3.7, 1.428_072_326_665_388), - (10.0, 12.801_827_480_081_469), - (100.5, 361.435_540_467_777_57), - (1000.0, 5_905.220_423_209_181), - ]; - for (x, expected) in cases { - let actual = ln_gamma(x); - if expected == 0.0 { - assert!(actual.abs() < 1e-13, "ln_gamma({x}) = {actual:.3e}, want 0"); - } else { - assert_close(actual, expected, 1e-12); - } - } +fn ensure_unit_interval(name: &'static str, value: f64) -> Result<(), SpecialError> { + if value.is_finite() && (0.0..=1.0).contains(&value) { + Ok(()) + } else { + Err(SpecialError::InvalidInput { + message: format!("{name} must be finite and in [0, 1]"), + }) } +} - #[test] - fn ln_gamma_invalid_domain_is_nan() { - assert!(ln_gamma(0.0).is_nan()); - assert!(ln_gamma(-1.5).is_nan()); +fn ensure_finite(name: &'static str, value: f64) -> Result<(), SpecialError> { + if value.is_finite() { + Ok(()) + } else { + Err(SpecialError::NumericalIssue { + message: format!("{name} was non-finite"), + }) } +} - // Reference values generated with: - // uv run python -c "from scipy import special; print(special.betainc(a, b, x))" - #[test] - fn betainc_reg_matches_scipy() { - let cases = [ - (0.5, 0.5, 0.3, 0.369_010_119_565_545_36), - (2.0, 3.0, 0.4, 0.524_799_999_999_999_9), - (5.5, 1.5, 0.7, 0.251_904_453_669_740_85), - (10.0, 10.0, 0.5, 0.5), - (0.5, 20.5, 0.01, 0.476_541_531_548_465_6), - (100.5, 900.5, 0.1, 0.494_385_987_853_672_66), - (3.5, 0.5, 0.99, 0.797_971_695_234_850_9), - ]; - for (a, b, x, expected) in cases { - assert_close(betainc_reg(a, b, x), expected, 1e-12); - } +fn clean_unit_probability(name: &'static str, value: f64) -> Result { + ensure_finite(name, value)?; + // Snap rounding excess just outside [0, 1] back to the boundary. The snap + // width is the crate's documented 1e-12 accuracy contract: a saturating + // series legitimately overshoots 1.0 by tens of ULPs (~2e-14 observed from + // the 48-term BGRAT sum), which exceeds any ULP-scale tolerance, while + // anything outside the contract width is a genuine numerical failure and + // must surface as a typed error, never a silent clamp. + let boundary_tolerance = 1.0e-12; + if (0.0..=1.0).contains(&value) { + Ok(value) + } else if value < 0.0 && value > -boundary_tolerance { + Ok(0.0) + } else if value > 1.0 && value < 1.0 + boundary_tolerance { + Ok(1.0) + } else { + Err(SpecialError::NumericalIssue { + message: format!("{name} was outside [0, 1]"), + }) } +} - #[test] - fn betainc_reg_invalid_domain_is_nan() { - assert!(betainc_reg(0.0, 1.0, 0.5).is_nan()); - assert!(betainc_reg(1.0, -1.0, 0.5).is_nan()); - assert!(betainc_reg(1.0, 1.0, -0.1).is_nan()); - assert!(betainc_reg(1.0, 1.0, 1.1).is_nan()); - } +fn is_strictly_inside(x: f64, lo: f64, hi: f64) -> bool { + x.is_finite() && x > lo && x < hi +} - // Reference values generated with: - // uv run python -c "from scipy import special; print(special.betaincinv(a, b, p))" - #[test] - fn betainc_inv_matches_scipy() { - let cases = [ - (0.5, 0.5, 0.25, 0.146_446_609_406_726_24), - (2.0, 3.0, 0.6, 0.444_500_002_083_767_4), - (5.5, 1.5, 0.05, 0.505_461_253_650_681_3), - (10.0, 10.0, 0.975, 0.711_356_752_083_001_1), - (0.5, 20.5, 0.995, 0.176_754_097_436_689_93), - (100.5, 900.5, 0.025, 0.082_562_652_843_060_04), - // Binomial-CI-scale parameters: n = 20000 trials, far tail. - (50.5, 19_950.5, 0.999_999, 0.004_581_655_467_494_118_5), - ]; - for (a, b, p, expected) in cases { - assert_close(betainc_inv(a, b, p), expected, 1e-8); - } +fn midpoint(lo: f64, hi: f64) -> f64 { + lo + 0.5 * (hi - lo) +} + +fn interiorize(x: f64) -> f64 { + if !x.is_finite() { + 0.5 + } else if x <= 0.0 { + f64::MIN_POSITIVE + } else if x >= 1.0 { + 1.0 - f64::EPSILON + } else { + x } +} - // Reference values generated with: - // uv run python -c "from scipy import special; print(special.betaincinv(a, b, p))" - // Upper-tail quantiles exercise the symmetry path (the direct Halley - // iteration loses the residual to cancellation beyond p ~ 1 - 1e-10). - #[test] - fn betainc_inv_upper_tail_matches_scipy() { - let cases: [(f64, f64, f64, f64); 4] = [ - (2.0, 3.0, 0.999_999_9, 0.997_073_840_091_498_9), - (100.5, 900.5, 1.0 - 1e-12, 0.179_649_794_238_11), - (7.5, 19_993.5, 1.0 - 2.3e-16, 0.002_728_615_291_135_757_6), - (0.5, 0.5, 0.999_999, 0.999_999_999_997_532_6), - ]; - for (a, b, p, expected) in cases { - let actual = betainc_inv(a, b, p); - let scale: f64 = expected.abs(); - assert!( - ((actual - expected).abs() / scale) < 1e-8, - "betainc_inv({a}, {b}, {p}): expected {expected:.12e}, got {actual:.12e}" - ); - } +fn direct_quantile(x: f64) -> BetaQuantile { + BetaQuantile { + x, + complement: 1.0 - x, } +} + +fn bracket_is_representation_limited(lo: f64, hi: f64, x: f64) -> bool { + hi - lo <= 2.0 * spacing_above(x).max(f64::MIN_POSITIVE) +} + +fn spacing_above(x: f64) -> f64 { + f64::from_bits(x.to_bits() + 1) - x +} + +#[allow(clippy::cast_precision_loss)] // Iteration budgets are small and exactly representable as f64 +fn usize_to_f64(value: usize) -> f64 { + value as f64 +} + +#[cfg(test)] +mod tests { + use super::{betai, invbetai, ln_gamma}; #[test] - fn betainc_inv_tails_are_symmetric() { - // Relative comparison, with cases chosen so neither side hits - // f64 representation limits: p stays >= 1e-6 (forming `1 - p` - // closer to 1 destroys p's precision before the function is even - // called) and the quantiles stay far enough from 0 and 1 that - // `1 - upper` keeps its significant digits. Outside those limits - // a mirrored comparison measures representation error, not - // implementation error. - let cases: [(f64, f64, &[f64]); 3] = [ - (2.0, 3.0, &[1e-6, 0.01, 0.3]), - (50.5, 19_950.5, &[1e-6, 0.01, 0.3]), - // Quantiles for this shape at small p sit below 1e-13, where - // the mirrored side cannot represent them; compare only at - // moderate p. - (0.5, 20.5, &[0.01, 0.3]), - ]; - for (a, b, ps) in cases { - for &p in ps { - let lower = betainc_inv(a, b, p); - let upper = betainc_inv(b, a, 1.0 - p); - let mirrored = 1.0 - upper; - assert!( - ((lower - mirrored) / lower).abs() < 1e-8, - "tail symmetry failed for a={a}, b={b}, p={p}: {lower} vs {mirrored}" - ); - } - } + fn ln_gamma_matches_factorial_case() { + let value = ln_gamma(6.0).unwrap(); + assert!((value - libm::log(120.0)).abs() < 1.0e-14); } #[test] - fn betainc_inv_round_trips_through_betainc_reg() { - for &(a, b) in &[(0.5, 0.5), (2.0, 3.0), (7.5, 19_993.5), (100.5, 900.5)] { - for &p in &[1e-6, 0.025, 0.5, 0.975, 1.0 - 1e-6] { - let x = betainc_inv(a, b, p); - let back = betainc_reg(a, b, x); - assert!( - (back - p).abs() < 1e-9, - "round trip failed for a={a}, b={b}, p={p}: x={x}, back={back}" - ); - } - } + fn betai_uses_reflection_symmetry() { + let a = 2.5; + let b = 7.5; + let x = 0.8; + let left = betai(a, b, x).unwrap(); + let right = 1.0 - betai(b, a, 1.0 - x).unwrap(); + assert!((left - right).abs() < 1.0e-13); } - // Allow exact float comparisons: the edge cases return the sentinel - // values 0.0 and 1.0 verbatim. - #[allow(clippy::float_cmp)] #[test] - fn betainc_inv_edges() { - assert_eq!(betainc_inv(2.0, 3.0, 0.0), 0.0); - assert_eq!(betainc_inv(2.0, 3.0, 1.0), 1.0); - assert!(betainc_inv(0.0, 3.0, 0.5).is_nan()); - assert!(betainc_inv(2.0, 3.0, -0.1).is_nan()); + fn invbetai_round_trips() { + let quantile = invbetai(0.025, 10.5, 90.5).unwrap(); + let probability = betai(10.5, 90.5, quantile.x).unwrap(); + assert!((probability - 0.025).abs() < 1.0e-14); } } diff --git a/crates/pecos-num/src/stats.rs b/crates/pecos-num/src/stats.rs index 533b691e3..b34b22a35 100644 --- a/crates/pecos-num/src/stats.rs +++ b/crates/pecos-num/src/stats.rs @@ -33,99 +33,233 @@ //! - [`jackknife_weighted`] - Jackknife resampling for weighted/grouped data (full workflow) //! - [`weighted_mean`] - Calculate weighted mean from (value, weight) pairs //! -//! ## Binomial Proportions -//! - [`jeffreys_interval`] - Jeffreys credible interval for a binomial proportion -//! //! The slice functions are fast and simple for 1D data. The axis functions //! provide idiomatic Rust API for multi-dimensional arrays. -use crate::special::betainc_inv; +use crate::special::{self, SpecialError}; use ndarray::{Array, ArrayView, Axis, Dimension, RemoveAxis}; +use std::fmt; -/// Jeffreys credible interval for a binomial proportion. -/// -/// Computes the equal-tailed interval of the Beta(k + 1/2, n - k + 1/2) -/// posterior arising from the Jeffreys prior Beta(1/2, 1/2), following -/// Brown, Cai & `DasGupta`, "Interval Estimation for a Binomial -/// Proportion", Statistical Science 16(2), 2001. Per that paper's -/// standard modification, the lower bound is 0 when `successes == 0` and -/// the upper bound is 1 when `successes == trials`. -/// -/// # Arguments -/// -/// * `successes` - Number of observed successes (k) -/// * `trials` - Number of trials (n), must be nonzero -/// * `confidence` - Interval mass, e.g. 0.95; must be in (0, 1) +const JEFFREYS_MAX_TRIALS: usize = 100_000_000; + +/// Posterior point estimator for a Jeffreys binomial model. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum JeffreysEstimator { + /// Posterior mean `(k + 1/2) / (n + 1)`. + Mean, + /// Posterior median `B^{-1}(0.5; k + 1/2, n - k + 1/2)`. + Median, +} + +/// Equal-tailed Jeffreys binomial-proportion interval and point estimate. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct JeffreysInterval { + /// Lower interval endpoint. + pub lo: f64, + /// Upper interval endpoint. + pub hi: f64, + /// Posterior-mean point estimate. + pub point: f64, + /// Interval width `hi - lo`. + pub width: f64, +} + +/// Error type for Jeffreys binomial-proportion routines. +#[derive(Debug, Clone, PartialEq)] +pub enum JeffreysError { + /// The number of trials was zero. + ZeroTrials, + /// The number of trials exceeded the verified numerical regime. + TrialsExceedSupported { + /// Number of trials. + n: usize, + /// Maximum supported number of trials. + max: usize, + }, + /// The number of successes exceeded the number of trials. + SuccessesExceedTrials { + /// Number of successes. + k: usize, + /// Number of trials. + n: usize, + }, + /// The interval tail probability was outside `(0, 1)`. + InvalidAlpha { + /// Invalid tail probability. + alpha: f64, + }, + /// A special-function routine failed. + SpecialFunction { + /// Underlying special-function error. + source: SpecialError, + }, + /// The computed interval violated the required endpoint ordering. + InvalidInterval { + /// Lower interval endpoint. + lo: f64, + /// Upper interval endpoint. + hi: f64, + }, +} + +impl fmt::Display for JeffreysError { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + match self { + Self::ZeroTrials => write!(f, "n must be greater than zero"), + Self::TrialsExceedSupported { n, max } => { + write!(f, "n must be less than or equal to {max}. Got n={n}") + } + Self::SuccessesExceedTrials { k, n } => { + write!(f, "k must be less than or equal to n. Got k={k}, n={n}") + } + Self::InvalidAlpha { alpha } => { + write!(f, "alpha must be finite and in (0, 1). Got alpha={alpha}") + } + Self::SpecialFunction { source } => write!(f, "Special-function error: {source}"), + Self::InvalidInterval { lo, hi } => { + write!(f, "computed Jeffreys interval is invalid: lo={lo}, hi={hi}") + } + } + } +} + +impl std::error::Error for JeffreysError {} + +impl From for JeffreysError { + fn from(source: SpecialError) -> Self { + Self::SpecialFunction { source } + } +} + +/// Compute the Jeffreys equal-tailed binomial-proportion interval. /// -/// # Returns +/// The posterior is `Beta(k + 1/2, n - k + 1/2)`. The interval endpoints are +/// the `alpha/2` and `1 - alpha/2` posterior quantiles, with the PECOS endpoint +/// rule `k = 0 -> lo = 0` and `k = n -> hi = 1`. The endpoint rule is applied +/// at this caller level; non-endpoint limits use [`special::invbetai`], whose +/// reflected complement solve avoids losing the directly solved tail to +/// arithmetic `1 - x` cancellation. /// -/// `(lower, upper)` bounds on the proportion. +/// The returned point estimate is the posterior mean. The verified numerical +/// regime is `n <= 100_000_000` and `alpha` in `[1e-6, 0.5]`; `n` is enforced +/// here, while `alpha` remains accepted for the full mathematical interval +/// `(0, 1)` and the verified band is documented for callers. /// -/// # Panics +/// # Errors /// -/// Panics if `trials == 0`, `trials > 10^12`, `successes > trials`, or -/// `confidence` is not in (0, 1). The first three are contract -/// violations; the trials cap marks the scale beyond which the -/// underlying incomplete-beta continued fraction has not been validated -/// (it can converge spuriously for extreme shape parameters, returning a -/// nonsense interval without warning). Also panics if the computed -/// bounds come back inverted — a numeric-breakdown trip-wire that should -/// be unreachable within the supported scales. +/// Returns an error if `n == 0`, `n > 100_000_000`, `k > n`, `alpha <= 0`, +/// `alpha >= 1`, `alpha` is non-finite, the underlying incomplete-beta inverse +/// fails, or the computed endpoints fail the interval postcondition. /// /// # Examples /// /// ``` /// use pecos_num::stats::jeffreys_interval; /// -/// let (lo, hi) = jeffreys_interval(50, 200, 0.95); -/// assert!(lo < 0.25 && 0.25 < hi); -/// -/// // Zero successes: lower bound is exactly 0. -/// let (lo, hi) = jeffreys_interval(0, 100, 0.95); -/// assert_eq!(lo, 0.0); -/// assert!(hi < 0.05); +/// let interval = jeffreys_interval(1, 2, 0.05).unwrap(); +/// assert!((interval.point - 0.5).abs() < 1e-15); +/// assert!(interval.lo < interval.point && interval.point < interval.hi); /// ``` -#[must_use] -#[allow(clippy::cast_precision_loss)] -// Cast is safe: the trials cap keeps counts far below f64 mantissa precision -pub fn jeffreys_interval(successes: u64, trials: u64, confidence: f64) -> (f64, f64) { - const MAX_TRIALS: u64 = 1_000_000_000_000; - - assert!(trials > 0, "jeffreys_interval requires trials > 0"); - assert!( - trials <= MAX_TRIALS, - "jeffreys_interval supports at most {MAX_TRIALS} trials (got {trials}); the \ - incomplete-beta evaluation is not validated beyond that scale" - ); - assert!( - successes <= trials, - "jeffreys_interval requires successes ({successes}) <= trials ({trials})" - ); - assert!( - confidence > 0.0 && confidence < 1.0, - "jeffreys_interval requires confidence in (0, 1), got {confidence}" - ); - - let alpha = 1.0 - confidence; - let a = successes as f64 + 0.5; - let b = (trials - successes) as f64 + 0.5; - - let lower = if successes == 0 { +pub fn jeffreys_interval( + k: usize, + n: usize, + alpha: f64, +) -> Result { + validate_jeffreys_inputs(k, n)?; + validate_alpha(alpha)?; + + let a = count_to_f64(k) + 0.5; + let b = count_to_f64(n - k) + 0.5; + let tail = 0.5 * alpha; + + let lo = if k == 0 { 0.0 } else { - betainc_inv(a, b, alpha / 2.0) + special::invbetai(tail, a, b)?.x }; - let upper = if successes == trials { + let hi = if k == n { 1.0 } else { - betainc_inv(a, b, 1.0 - alpha / 2.0) + special::invbetai(1.0 - tail, a, b)?.x }; - assert!( - lower <= upper, - "jeffreys_interval produced inverted bounds ({lower} > {upper}) for k={successes}, \ - n={trials}; this indicates incomplete-beta breakdown and is a bug" - ); - (lower, upper) + let point = jeffreys_point(k, n, JeffreysEstimator::Mean)?; + validate_interval_postcondition(lo, hi)?; + + Ok(JeffreysInterval { + lo, + hi, + point, + width: hi - lo, + }) +} + +/// Compute a Jeffreys binomial-proportion posterior point estimate. +/// +/// # Errors +/// +/// Returns an error if `n == 0`, `n > 100_000_000`, `k > n`, or the median +/// estimator's incomplete-beta inverse fails. +/// +/// # Examples +/// +/// ``` +/// use pecos_num::stats::{jeffreys_point, JeffreysEstimator}; +/// +/// let point = jeffreys_point(3, 10, JeffreysEstimator::Mean).unwrap(); +/// assert_eq!(point, 3.5 / 11.0); +/// ``` +pub fn jeffreys_point( + k: usize, + n: usize, + estimator: JeffreysEstimator, +) -> Result { + validate_jeffreys_inputs(k, n)?; + + match estimator { + JeffreysEstimator::Mean => Ok((count_to_f64(k) + 0.5) / (count_to_f64(n) + 1.0)), + JeffreysEstimator::Median => { + let a = count_to_f64(k) + 0.5; + let b = count_to_f64(n - k) + 0.5; + Ok(special::invbetai(0.5, a, b)?.x) + } + } +} + +fn validate_jeffreys_inputs(k: usize, n: usize) -> Result<(), JeffreysError> { + if n == 0 { + return Err(JeffreysError::ZeroTrials); + } + if n > JEFFREYS_MAX_TRIALS { + return Err(JeffreysError::TrialsExceedSupported { + n, + max: JEFFREYS_MAX_TRIALS, + }); + } + if k > n { + return Err(JeffreysError::SuccessesExceedTrials { k, n }); + } + Ok(()) +} + +fn validate_alpha(alpha: f64) -> Result<(), JeffreysError> { + if alpha.is_finite() && alpha > 0.0 && alpha < 1.0 { + Ok(()) + } else { + Err(JeffreysError::InvalidAlpha { alpha }) + } +} + +fn validate_interval_postcondition(lo: f64, hi: f64) -> Result<(), JeffreysError> { + if lo.is_finite() && hi.is_finite() && lo >= 0.0 && lo <= hi && hi <= 1.0 { + Ok(()) + } else { + Err(JeffreysError::InvalidInterval { lo, hi }) + } +} + +#[allow(clippy::cast_precision_loss)] // Jeffreys support is n <= 1e8, exactly representable as f64 +fn count_to_f64(value: usize) -> f64 { + value as f64 } /// Calculate the arithmetic mean of a slice of values. @@ -692,94 +826,6 @@ mod tests { use super::*; use ndarray::Axis; - // Reference values generated with: - // uv run python -c "from scipy import stats; - // print(stats.beta.ppf(q, k + 0.5, n - k + 0.5))" - #[test] - fn jeffreys_interval_matches_scipy_beta_quantiles() { - let cases: [(u64, u64, f64, f64, f64); 6] = [ - (0, 100, 0.95, 0.0, 0.024_745_270_015_269_89), - (100, 100, 0.95, 0.975_254_729_984_730_1, 1.0), - ( - 3, - 1000, - 0.95, - 0.000_845_634_801_829_834_8, - 0.007_984_367_358_403_443, - ), - ( - 50, - 200, - 0.95, - 0.193_872_680_411_726_73, - 0.313_302_662_892_847_86, - ), - // High-confidence intervals at LER-study scales. - ( - 7, - 20_000, - 0.99999, - 3.838_996_822_347_517e-5, - 1.307_358_543_951_447_7e-3, - ), - ( - 1234, - 20_000, - 0.99999, - 0.054_475_933_954_188_78, - 0.069_508_522_504_658_04, - ), - ]; - for (k, n, conf, lo_expected, hi_expected) in cases { - let (lo, hi) = jeffreys_interval(k, n, conf); - let lo_scale = lo_expected.abs().max(1e-12); - let hi_scale = hi_expected.abs().max(1e-12); - assert!( - (lo - lo_expected).abs() / lo_scale < 1e-6, - "lower bound for k={k}, n={n}: expected {lo_expected:.12e}, got {lo:.12e}" - ); - assert!( - (hi - hi_expected).abs() / hi_scale < 1e-6, - "upper bound for k={k}, n={n}: expected {hi_expected:.12e}, got {hi:.12e}" - ); - } - } - - #[test] - fn jeffreys_interval_brackets_the_point_estimate() { - let (lo, hi) = jeffreys_interval(50, 200, 0.95); - assert!(lo < 0.25 && 0.25 < hi); - // Wider confidence gives a wider interval. - let (lo99, hi99) = jeffreys_interval(50, 200, 0.99); - assert!(lo99 < lo && hi < hi99); - } - - #[test] - #[should_panic(expected = "trials > 0")] - fn jeffreys_interval_rejects_zero_trials() { - let _ = jeffreys_interval(0, 0, 0.95); - } - - #[test] - #[should_panic(expected = "at most")] - fn jeffreys_interval_rejects_unvalidated_trial_scales() { - // Beyond ~1e12 trials the incomplete-beta continued fraction can - // converge spuriously (observed at 2e15: inverted bounds). - let _ = jeffreys_interval(1_000_000_000_000_000, 2_000_000_000_000_000, 0.95); - } - - #[test] - #[should_panic(expected = "successes")] - fn jeffreys_interval_rejects_successes_above_trials() { - let _ = jeffreys_interval(11, 10, 0.95); - } - - #[test] - #[should_panic(expected = "confidence")] - fn jeffreys_interval_rejects_bad_confidence() { - let _ = jeffreys_interval(5, 10, 1.0); - } - // Allow exact float comparisons in tests - we're testing mathematically exact results // that are exactly representable in IEEE 754 (e.g., 3.0, 42.0, 0.4) #[allow(clippy::float_cmp)] diff --git a/crates/pecos-num/tests/jeffreys.rs b/crates/pecos-num/tests/jeffreys.rs new file mode 100644 index 000000000..3c2dad5a2 --- /dev/null +++ b/crates/pecos-num/tests/jeffreys.rs @@ -0,0 +1,572 @@ +// Copyright 2026 The PECOS Developers +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// https://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +use pecos_num::special::{self, InverseBetaOptions, SpecialError}; +use pecos_num::stats::{JeffreysError, JeffreysEstimator, jeffreys_interval, jeffreys_point}; + +const FIXTURE_CSV: &str = include_str!("fixtures/jeffreys_scipy.csv"); +const U01_BITS: u64 = 0x3ff0_0000_0000_0000; + +#[derive(Debug, Clone, Copy)] +struct FixtureRow { + k: usize, + n: usize, + alpha: f64, + lo: f64, + hi: f64, + median: f64, +} + +#[derive(Debug, Clone, Copy)] +struct WorstError { + relative_error: f64, + k: usize, + n: usize, + alpha: f64, + bound: &'static str, +} + +#[derive(Debug, Clone, Copy)] +struct WorstResidual { + residual: f64, + k: usize, + n: usize, + alpha: f64, + bound: &'static str, +} + +impl Default for WorstError { + fn default() -> Self { + Self { + relative_error: 0.0, + k: 0, + n: 0, + alpha: 0.0, + bound: "", + } + } +} + +impl Default for WorstResidual { + fn default() -> Self { + Self { + residual: 0.0, + k: 0, + n: 0, + alpha: 0.0, + bound: "", + } + } +} + +#[test] +fn fixture_table_matches_scipy() { + let mut worst = WorstError::default(); + + for row in fixture_rows() { + let interval = jeffreys_interval(row.k, row.n, row.alpha).unwrap_or_else(|err| { + panic!( + "interval failed for k={}, n={}, alpha={}: {err}", + row.k, row.n, row.alpha + ) + }); + check_fixture_value(&mut worst, row, "lo", interval.lo, row.lo); + check_fixture_value(&mut worst, row, "hi", interval.hi, row.hi); + + let median = + jeffreys_point(row.k, row.n, JeffreysEstimator::Median).unwrap_or_else(|err| { + panic!( + "median failed for k={}, n={}, alpha={}: {err}", + row.k, row.n, row.alpha + ) + }); + check_fixture_value(&mut worst, row, "median", median, row.median); + } + + eprintln!( + "max scipy fixture relative error: {} at k={}, n={}, alpha={}, bound={}", + worst.relative_error, worst.k, worst.n, worst.alpha, worst.bound + ); +} + +#[test] +fn forward_betai_near_one_saturation_matches_scipy() { + // Quantile fixtures never evaluate the forward CDF in the near-1 saturation + // band inside the asymptotic gate: invbetai converges in the tail that stays + // cleanly below 1, so a forward-side defect there is invisible to every + // interval test (a boundary-snap ordering bug survived four review rounds + // this way, collapsing I_x ~ 1 to 0.0). Oracle: scipy.special.betainc. + let cases = [ + (30.5, 10_001.0, 0.0095, 0.999_999_999_999_997_7), + (31.2, 10_050.0, 0.0099, 0.999_999_999_999_999_8), + (10.5, 10_050.0, 0.006, 0.999_999_999_999_999_6), + (2.5, 20_000.0, 0.002, 0.999_999_999_999_999_2), + (44.6, 20_000.0, 0.006, 0.999_999_999_999_999_3), + ]; + for (a, b, x, expected) in cases { + let value = special::betai(a, b, x) + .unwrap_or_else(|err| panic!("betai({a}, {b}, {x}) failed: {err}")); + let relative = (value - expected).abs() / expected; + assert!( + relative <= 1.0e-12, + "betai({a}, {b}, {x}) = {value}, expected {expected}, rel {relative}" + ); + } +} + +#[test] +fn forward_betai_residual_at_fixture_quantiles_is_bounded() { + let mut worst = WorstResidual::default(); + + for row in fixture_rows() { + let a = usize_to_f64(row.k) + 0.5; + let b = usize_to_f64(row.n - row.k) + 0.5; + let tail = 0.5 * row.alpha; + + if row.k > 0 { + check_forward_residual(&mut worst, row, "lo", a, b, row.lo, tail); + } + if row.k < row.n { + check_forward_residual(&mut worst, row, "hi", a, b, row.hi, 1.0 - tail); + } + if !uses_large_shape_median_shortcut(a, b) { + check_forward_residual(&mut worst, row, "median", a, b, row.median, 0.5); + } + } + + eprintln!( + "max forward betai fixture residual: {} at k={}, n={}, alpha={}, bound={}", + worst.residual, worst.k, worst.n, worst.alpha, worst.bound + ); + assert!( + worst.residual <= 5.0e-9, + "forward betai residual too large: {worst:?}" + ); +} + +#[test] +fn symmetry_round_trip_holds() { + for &(n, alpha) in &[(1, 0.05), (2, 0.05), (101, 0.01), (257, 0.1), (503, 1.0e-6)] { + for k in 0..=n { + let interval = jeffreys_interval(k, n, alpha).unwrap(); + let reflected = jeffreys_interval(n - k, n, alpha).unwrap(); + assert_close(interval.lo, 1.0 - reflected.hi, 2.0e-9, "lower symmetry"); + assert_close(interval.hi, 1.0 - reflected.lo, 2.0e-9, "upper symmetry"); + } + } +} + +#[test] +fn bounds_are_monotone_in_k() { + for &(n, alpha) in &[(200, 0.05), (200, 0.01), (200, 0.1), (200, 1.0e-6)] { + let first = jeffreys_interval(0, n, alpha).unwrap(); + let mut previous_lo = first.lo; + let mut previous_hi = first.hi; + + for k in 1..=n { + let interval = jeffreys_interval(k, n, alpha).unwrap(); + assert!( + interval.lo + 1.0e-14 >= previous_lo, + "lo not monotone at k={k}, n={n}, alpha={alpha}: {} < {previous_lo}", + interval.lo + ); + assert!( + interval.hi + 1.0e-14 >= previous_hi, + "hi not monotone at k={k}, n={n}, alpha={alpha}: {} < {previous_hi}", + interval.hi + ); + previous_lo = interval.lo; + previous_hi = interval.hi; + } + } +} + +#[test] +fn gate_sweeps_remain_monotone_and_ordered() { + check_k_sweep_monotone(60, 70, 100_000_000, 0.05); + check_k_sweep_monotone(2_990, 3_030, 1_000_000, 0.05); +} + +#[test] +fn alpha_sweeps_are_continuous_near_regime_gates() { + for &(k, n) in &[ + (64, 100_000_000), + (65, 100_000_000), + (3_000, 1_000_000), + (3_016, 1_000_000), + ] { + let alphas = [1.0e-6, 3.0e-6, 1.0e-5, 3.0e-5, 1.0e-4, 1.0e-3, 0.01, 0.05]; + let mut previous = jeffreys_interval(k, n, alphas[0]).unwrap(); + assert_interval_ordered(previous, k, n, alphas[0]); + + for &alpha in &alphas[1..] { + let interval = jeffreys_interval(k, n, alpha).unwrap(); + assert_interval_ordered(interval, k, n, alpha); + assert!( + interval.lo >= previous.lo, + "lo widened as alpha increased at k={k}, n={n}, alpha={alpha}" + ); + assert!( + interval.hi <= previous.hi, + "hi widened as alpha increased at k={k}, n={n}, alpha={alpha}" + ); + assert!( + relative_x_error(interval.lo, previous.lo).is_finite(), + "lo discontinuity at k={k}, n={n}, alpha={alpha}" + ); + assert!( + relative_x_error(interval.hi, previous.hi).is_finite(), + "hi discontinuity at k={k}, n={n}, alpha={alpha}" + ); + previous = interval; + } + } +} + +#[test] +fn large_b_gate_n_sweep_is_ordered() { + let k = 64; + let alpha = 0.05; + let trial_counts = [ + 49_999_990, + 50_000_000, + 50_000_063, + 50_000_064, + 50_000_100, + 100_000_000, + ]; + + for &n in &trial_counts { + let interval = jeffreys_interval(k, n, alpha).unwrap(); + assert_interval_ordered(interval, k, n, alpha); + } +} + +#[test] +fn small_a_low_b_pocket_sweeps_are_ordered() { + for &n in &[1_000_000, 3_000_000, 6_000_000, 8_900_002, 8_999_990] { + check_k_sweep_monotone(2, 40, n, 0.05); + } +} + +#[test] +fn large_b_floor_side_sweeps_are_ordered() { + for &(k, counts) in &[ + (2, &[9_000, 9_999, 10_000, 10_001, 30_000, 100_000][..]), + (64, &[40_000, 41_000, 42_000, 100_000, 1_000_000][..]), + ] { + for &n in counts { + let interval = jeffreys_interval(k, n, 0.05).unwrap(); + assert_interval_ordered(interval, k, n, 0.05); + } + } +} + +#[test] +fn large_b_saturation_keeps_inverse_orientation() { + let cdf_at_half = special::betai(5.5, 2_999_995.5, 0.5).unwrap(); + assert!( + cdf_at_half > 1.0 - 1.0e-15, + "large-b half CDF should saturate high, got {cdf_at_half}" + ); + + let cdf_at_guard = special::betai(5.5, 2_999_995.5, 1.0e-3).unwrap(); + assert!( + cdf_at_guard > 1.0 - 1.0e-15, + "large-b guard CDF should saturate high, got {cdf_at_guard}" + ); + + let interval = jeffreys_interval(5, 3_000_000, 0.05).unwrap(); + assert_interval_ordered(interval, 5, 3_000_000, 0.05); +} + +#[test] +fn gl_tail_orientation_regression_has_positive_width() { + for &(k, n, alpha) in &[ + (3_016, 1_000_000, 0.05), + (15_881, 20_000, 0.05), + (67_867_393, 100_000_000, 0.05), + ] { + let interval = jeffreys_interval(k, n, alpha).unwrap(); + assert_interval_ordered(interval, k, n, alpha); + } +} + +#[test] +fn endpoints_are_bit_exact() { + for &(n, alpha) in &[(1, 0.05), (2, 0.05), (100, 0.01), (100_000_000, 1.0e-6)] { + let zero = jeffreys_interval(0, n, alpha).unwrap(); + assert_eq!(zero.lo.to_bits(), 0.0_f64.to_bits()); + + let full = jeffreys_interval(n, n, alpha).unwrap(); + assert_eq!(full.hi.to_bits(), 1.0_f64.to_bits()); + } +} + +#[test] +fn posterior_mean_point_estimate_is_used() { + for &(k, n) in &[(0, 1), (1, 2), (10, 100), (1_000, 1_000_000)] { + let expected = (usize_to_f64(k) + 0.5) / (usize_to_f64(n) + 1.0); + let point = jeffreys_point(k, n, JeffreysEstimator::Mean).unwrap(); + assert_eq!(point.to_bits(), expected.to_bits()); + + let interval = jeffreys_interval(k, n, 0.05).unwrap(); + assert_eq!(interval.point.to_bits(), expected.to_bits()); + } +} + +#[test] +fn typed_errors_are_returned() { + assert!(matches!( + jeffreys_interval(0, 0, 0.05), + Err(JeffreysError::ZeroTrials) + )); + assert!(matches!( + jeffreys_interval(2, 1, 0.05), + Err(JeffreysError::SuccessesExceedTrials { k: 2, n: 1 }) + )); + assert!(matches!( + jeffreys_interval(0, 100_000_001, 0.05), + Err(JeffreysError::TrialsExceedSupported { + n: 100_000_001, + max: 100_000_000 + }) + )); + assert!(matches!( + jeffreys_interval(0, 1, 0.0), + Err(JeffreysError::InvalidAlpha { alpha: 0.0 }) + )); + assert!(matches!( + jeffreys_interval(0, 1, 1.0), + Err(JeffreysError::InvalidAlpha { alpha: 1.0 }) + )); + + let options = InverseBetaOptions { + max_newton: 0, + max_bisect: 0, + ..InverseBetaOptions::default() + }; + assert!(matches!( + special::invbetai_with_options(0.5, 2.0, 2.0, options), + Err(SpecialError::MaxIterations { + function: "invbetai", + iterations: 0 + }) + )); +} + +#[test] +fn deterministic_randomized_grid_matches_puruspe_dev_check() { + // scipy fixtures are the oracle of record for Jeffreys endpoints. This + // deterministic puruspe sweep is only a dev-only differential check; on any + // disagreement, the pinned scipy fixture table takes precedence. + let mut rng = Lcg::new(0x4a45_4646_5245_5953); + + for _ in 0..200 { + let n = rng.next_usize(1_000) + 2; + let k = rng.next_usize(n - 2) + 1; + let alpha = 1.0e-4 + rng.next_unit_f64() * (0.25 - 1.0e-4); + let interval = jeffreys_interval(k, n, alpha).unwrap(); + let median = jeffreys_point(k, n, JeffreysEstimator::Median).unwrap(); + + let a = usize_to_f64(k) + 0.5; + let b = usize_to_f64(n - k) + 0.5; + let tail = 0.5 * alpha; + let puruspe_lo = if k == 0 { + 0.0 + } else { + puruspe::invbetai(tail, a, b) + }; + let puruspe_hi = if k == n { + 1.0 + } else { + puruspe::invbetai(1.0 - tail, a, b) + }; + let puruspe_median = puruspe::invbetai(0.5, a, b); + + assert_close(interval.lo, puruspe_lo, 1.0e-5, "puruspe lo"); + assert_close(interval.hi, puruspe_hi, 1.0e-5, "puruspe hi"); + assert_close(median, puruspe_median, 1.0e-5, "puruspe median"); + } +} + +fn fixture_rows() -> Vec { + FIXTURE_CSV.lines().skip(1).map(parse_fixture_row).collect() +} + +fn parse_fixture_row(line: &str) -> FixtureRow { + let mut parts = line.split(','); + let row = FixtureRow { + k: parts.next().unwrap().parse().unwrap(), + n: parts.next().unwrap().parse().unwrap(), + alpha: parts.next().unwrap().parse().unwrap(), + lo: parts.next().unwrap().parse().unwrap(), + hi: parts.next().unwrap().parse().unwrap(), + median: parts.next().unwrap().parse().unwrap(), + }; + assert!(parts.next().is_none(), "too many CSV columns in {line}"); + row +} + +fn check_fixture_value( + worst: &mut WorstError, + row: FixtureRow, + bound: &'static str, + actual: f64, + expected: f64, +) { + // The design tolerance is expressed as + // max(1e-12, 5e-15 / max(p, 1-p)) on relative x error. The fixture columns + // are x-values, so this test interprets p as x_oracle and uses + // max(x_oracle, 1 - x_oracle). + let relative_error = relative_x_error(actual, expected); + let tolerance = fixture_relative_tolerance(expected); + if relative_error > worst.relative_error { + *worst = WorstError { + relative_error, + k: row.k, + n: row.n, + alpha: row.alpha, + bound, + }; + } + assert!( + relative_error <= tolerance, + "{bound} mismatch for k={}, n={}, alpha={}: actual={actual}, expected={expected}, rel={relative_error}, tol={tolerance}", + row.k, + row.n, + row.alpha, + ); +} + +fn check_forward_residual( + worst: &mut WorstResidual, + row: FixtureRow, + bound: &'static str, + a: f64, + b: f64, + x: f64, + target: f64, +) { + let actual = special::betai(a, b, x).unwrap_or_else(|err| { + panic!( + "forward betai failed for {bound}, k={}, n={}, alpha={}: {err}", + row.k, row.n, row.alpha + ) + }); + let residual = (actual - target).abs(); + if residual > worst.residual { + *worst = WorstResidual { + residual, + k: row.k, + n: row.n, + alpha: row.alpha, + bound, + }; + } +} + +fn fixture_relative_tolerance(expected: f64) -> f64 { + 1.0e-12_f64.max(5.0e-15 / expected.max(1.0 - expected)) +} + +fn relative_x_error(actual: f64, expected: f64) -> f64 { + if expected.to_bits() == 0.0_f64.to_bits() { + if actual.to_bits() == 0.0_f64.to_bits() { + 0.0 + } else { + f64::INFINITY + } + } else { + (actual - expected).abs() / expected.abs() + } +} + +fn assert_close(actual: f64, expected: f64, relative_tolerance: f64, label: &str) { + let relative_error = relative_x_error(actual, expected); + assert!( + relative_error <= relative_tolerance, + "{label}: actual={actual}, expected={expected}, rel={relative_error}, tol={relative_tolerance}" + ); +} + +fn check_k_sweep_monotone(k_start: usize, k_end: usize, n: usize, alpha: f64) { + let first = jeffreys_interval(k_start, n, alpha).unwrap(); + assert_interval_ordered(first, k_start, n, alpha); + let mut previous_lo = first.lo; + let mut previous_hi = first.hi; + + for k in (k_start + 1)..=k_end { + let interval = jeffreys_interval(k, n, alpha).unwrap(); + assert_interval_ordered(interval, k, n, alpha); + assert!( + interval.lo + 1.0e-14 >= previous_lo, + "lo not monotone across gate at k={k}, n={n}, alpha={alpha}" + ); + assert!( + interval.hi + 1.0e-14 >= previous_hi, + "hi not monotone across gate at k={k}, n={n}, alpha={alpha}" + ); + previous_lo = interval.lo; + previous_hi = interval.hi; + } +} + +fn assert_interval_ordered( + interval: pecos_num::stats::JeffreysInterval, + k: usize, + n: usize, + alpha: f64, +) { + assert!( + interval.lo >= 0.0 && interval.lo <= interval.hi && interval.hi <= 1.0, + "invalid interval at k={k}, n={n}, alpha={alpha}: {interval:?}" + ); +} + +fn uses_large_shape_median_shortcut(a: f64, b: f64) -> bool { + a >= 1_000_000.0 && b >= 1_000_000.0 +} + +#[allow(clippy::cast_precision_loss)] // Test counts are <= 1e8, exactly representable as f64 +fn usize_to_f64(value: usize) -> f64 { + value as f64 +} + +struct Lcg { + state: u64, +} + +impl Lcg { + fn new(seed: u64) -> Self { + Self { state: seed } + } + + fn next_u64(&mut self) -> u64 { + self.state = self + .state + .wrapping_mul(6_364_136_223_846_793_005) + .wrapping_add(1_442_695_040_888_963_407); + self.state + } + + fn next_unit_f64(&mut self) -> f64 { + f64::from_bits(U01_BITS | (self.next_u64() >> 12)) - 1.0 + } + + fn next_usize(&mut self, max_inclusive: usize) -> usize { + let span = u64::try_from(max_inclusive).unwrap() + 1; + usize::try_from(self.next_u64() % span).unwrap() + } +} From 5724d21ccbbadb5aac21d965775ed5a8de809250 Mon Sep 17 00:00:00 2001 From: Ciaran Ryan-Anderson Date: Mon, 6 Jul 2026 21:37:51 -0600 Subject: [PATCH 06/11] Adopt scipy-aligned special-function names and u64 counts; migrate neo tests to the typed Jeffreys API --- crates/pecos-num/src/lib.rs | 1 + crates/pecos-num/src/special.rs | 66 +++++++++---------- crates/pecos-num/src/stats.rs | 34 ++++------ crates/pecos-num/tests/jeffreys.rs | 64 ++++++++++-------- crates/pecos/tests/neo_emission_test.rs | 56 +++++++++------- .../tests/neo_equivalence_matrix_test.rs | 15 +++-- crates/pecos/tests/neo_surface_ler_test.rs | 35 ++++++---- .../pecos/tests/neo_v6_example_sweep_test.rs | 11 ++-- 8 files changed, 151 insertions(+), 131 deletions(-) diff --git a/crates/pecos-num/src/lib.rs b/crates/pecos-num/src/lib.rs index df965af51..89bb5238f 100644 --- a/crates/pecos-num/src/lib.rs +++ b/crates/pecos-num/src/lib.rs @@ -57,6 +57,7 @@ pub use graph::Graph; pub use linalg::{matrix_exp, matrix_log}; pub use optimize::{BrentqOptions, NewtonOptions, OptimizeError, brentq, newton}; pub use polynomial::{Poly1d, PolynomialError, polyfit}; +pub use special::{betainc_inv, betainc_reg, ln_gamma}; pub use stats::{ JeffreysError, JeffreysEstimator, JeffreysInterval, jeffreys_interval, jeffreys_point, mean, }; diff --git a/crates/pecos-num/src/special.rs b/crates/pecos-num/src/special.rs index ac53df4b2..23f5ef5dd 100644 --- a/crates/pecos-num/src/special.rs +++ b/crates/pecos-num/src/special.rs @@ -259,7 +259,7 @@ impl Default for InverseBetaOptions { /// Inverse-beta quantile and its complement. /// -/// For roots at or below `0.5`, [`invbetai`] solves directly for [`Self::x`] +/// For roots at or below `0.5`, [`betainc_inv`] solves directly for [`Self::x`] /// and derives [`Self::complement`]. For roots above `0.5`, it solves the /// swapped complement problem directly for [`Self::complement`] and derives /// [`Self::x`]. This reconciles the design-note requirement to avoid losing the @@ -331,13 +331,13 @@ pub fn ln_gamma(x: f64) -> Result { /// # Examples /// /// ``` -/// use pecos_num::special::betai; +/// use pecos_num::special::betainc_reg; /// -/// let value = betai(2.0, 2.0, 0.5).unwrap(); +/// let value = betainc_reg(2.0, 2.0, 0.5).unwrap(); /// assert!((value - 0.5).abs() < 1e-14); /// ``` -pub fn betai(a: f64, b: f64, x: f64) -> Result { - Ok(betai_pair(a, b, x)?.lower) +pub fn betainc_reg(a: f64, b: f64, x: f64) -> Result { + Ok(betainc_reg_pair(a, b, x)?.lower) } /// Inverse regularized incomplete beta function. @@ -358,13 +358,13 @@ pub fn betai(a: f64, b: f64, x: f64) -> Result { /// # Examples /// /// ``` -/// use pecos_num::special::invbetai; +/// use pecos_num::special::betainc_inv; /// -/// let quantile = invbetai(0.5, 2.0, 2.0).unwrap(); +/// let quantile = betainc_inv(0.5, 2.0, 2.0).unwrap(); /// assert!((quantile.x - 0.5).abs() < 1e-14); /// ``` -pub fn invbetai(p: f64, a: f64, b: f64) -> Result { - invbetai_with_options(p, a, b, InverseBetaOptions::default()) +pub fn betainc_inv(p: f64, a: f64, b: f64) -> Result { + betainc_inv_with_options(p, a, b, InverseBetaOptions::default()) } /// Inverse regularized incomplete beta function with explicit iteration options. @@ -376,7 +376,7 @@ pub fn invbetai(p: f64, a: f64, b: f64) -> Result { /// /// Returns an error if the inputs or options are invalid, a non-finite value is /// encountered, or the configured iteration budget is exhausted. -pub fn invbetai_with_options( +pub fn betainc_inv_with_options( p: f64, a: f64, b: f64, @@ -404,19 +404,19 @@ pub fn invbetai_with_options( return Ok(quantile); } - let half = betai_pair(a, b, 0.5)?; + let half = betainc_reg_pair(a, b, 0.5)?; if p > half.lower { - let complement = invbetai_lower_tail(1.0 - p, b, a, options)?; + let complement = betainc_inv_lower_tail(1.0 - p, b, a, options)?; return Ok(BetaQuantile { x: complement.complement, complement: complement.x, }); } - invbetai_lower_tail(p, a, b, options) + betainc_inv_lower_tail(p, a, b, options) } -fn invbetai_lower_tail( +fn betainc_inv_lower_tail( p: f64, a: f64, b: f64, @@ -425,14 +425,14 @@ fn invbetai_lower_tail( let total_iterations = options.max_newton + options.max_bisect; if total_iterations == 0 { return Err(SpecialError::MaxIterations { - function: "invbetai", + function: "betainc_inv", iterations: 0, }); } let mut lo = 0.0; let mut hi = 1.0; - let mut x = initial_invbetai_guess(p, a, b); + let mut x = initial_betainc_inv_guess(p, a, b); if !is_strictly_inside(x, lo, hi) { x = midpoint(lo, hi); } @@ -444,7 +444,7 @@ fn invbetai_lower_tail( let effective_probability_tolerance = effective_probability_tolerance(options, a, b); for iteration in 0..total_iterations { - let residual = invbetai_residual(p, a, b, x)?; + let residual = betainc_inv_residual(p, a, b, x)?; let abs_residual = residual.abs(); if same_float(abs_residual, 0.0) { @@ -491,7 +491,7 @@ fn invbetai_lower_tail( return Ok(direct_quantile(bisected)); } return Err(SpecialError::MaxIterations { - function: "invbetai", + function: "betainc_inv", iterations: iteration + 1, }); } @@ -500,7 +500,7 @@ fn invbetai_lower_tail( if same_float(candidate, x) { return Err(SpecialError::MaxIterations { - function: "invbetai", + function: "betainc_inv", iterations: iteration + 1, }); } @@ -511,7 +511,7 @@ fn invbetai_lower_tail( } Err(SpecialError::MaxIterations { - function: "invbetai", + function: "betainc_inv", iterations: total_iterations, }) } @@ -584,7 +584,7 @@ fn same_float(left: f64, right: f64) -> bool { left.to_bits() == right.to_bits() } -fn betai_pair(a: f64, b: f64, x: f64) -> Result { +fn betainc_reg_pair(a: f64, b: f64, x: f64) -> Result { ensure_positive("a", a)?; ensure_positive("b", b)?; ensure_unit_interval("x", x)?; @@ -606,7 +606,7 @@ fn betai_pair(a: f64, b: f64, x: f64) -> Result { return Ok(cdf); } if a > BETA_GAUSS_LEGENDRE_SWITCH && b > BETA_GAUSS_LEGENDRE_SWITCH { - return betai_gauss_legendre_pair(a, b, x); + return betainc_reg_gauss_legendre_pair(a, b, x); } let prefactor = beta_prefactor(a, b, x)?; @@ -669,7 +669,7 @@ fn beta_power_series(a: f64, b: f64, x: f64) -> Result { }) } -fn betai_gauss_legendre_pair(a: f64, b: f64, x: f64) -> Result { +fn betainc_reg_gauss_legendre_pair(a: f64, b: f64, x: f64) -> Result { let a1 = a - 1.0; let b1 = b - 1.0; let mu = a / (a + b); @@ -1066,8 +1066,8 @@ fn betacf(a: f64, b: f64, x: f64) -> Result { }) } -fn invbetai_residual(p: f64, a: f64, b: f64, x: f64) -> Result { - let cdf = betai_pair(a, b, x)?; +fn betainc_inv_residual(p: f64, a: f64, b: f64, x: f64) -> Result { + let cdf = betainc_reg_pair(a, b, x)?; if p <= 0.5 { Ok(cdf.lower - p) } else { @@ -1175,7 +1175,7 @@ fn stirling_correction(x: f64) -> f64 { + inverse13 / 156.0 } -fn initial_invbetai_guess(p: f64, a: f64, b: f64) -> f64 { +fn initial_betainc_inv_guess(p: f64, a: f64, b: f64) -> f64 { if a >= 1.0 && b >= 1.0 { let pp = if p < 0.5 { p } else { 1.0 - p }; let t = libm::sqrt(-2.0 * libm::log(pp)); @@ -1297,7 +1297,7 @@ fn usize_to_f64(value: usize) -> f64 { #[cfg(test)] mod tests { - use super::{betai, invbetai, ln_gamma}; + use super::{betainc_inv, betainc_reg, ln_gamma}; #[test] fn ln_gamma_matches_factorial_case() { @@ -1306,19 +1306,19 @@ mod tests { } #[test] - fn betai_uses_reflection_symmetry() { + fn betainc_reg_uses_reflection_symmetry() { let a = 2.5; let b = 7.5; let x = 0.8; - let left = betai(a, b, x).unwrap(); - let right = 1.0 - betai(b, a, 1.0 - x).unwrap(); + let left = betainc_reg(a, b, x).unwrap(); + let right = 1.0 - betainc_reg(b, a, 1.0 - x).unwrap(); assert!((left - right).abs() < 1.0e-13); } #[test] - fn invbetai_round_trips() { - let quantile = invbetai(0.025, 10.5, 90.5).unwrap(); - let probability = betai(10.5, 90.5, quantile.x).unwrap(); + fn betainc_inv_round_trips() { + let quantile = betainc_inv(0.025, 10.5, 90.5).unwrap(); + let probability = betainc_reg(10.5, 90.5, quantile.x).unwrap(); assert!((probability - 0.025).abs() < 1.0e-14); } } diff --git a/crates/pecos-num/src/stats.rs b/crates/pecos-num/src/stats.rs index b34b22a35..9cf534b43 100644 --- a/crates/pecos-num/src/stats.rs +++ b/crates/pecos-num/src/stats.rs @@ -40,7 +40,7 @@ use crate::special::{self, SpecialError}; use ndarray::{Array, ArrayView, Axis, Dimension, RemoveAxis}; use std::fmt; -const JEFFREYS_MAX_TRIALS: usize = 100_000_000; +const JEFFREYS_MAX_TRIALS: u64 = 100_000_000; /// Posterior point estimator for a Jeffreys binomial model. #[derive(Debug, Clone, Copy, PartialEq, Eq)] @@ -72,16 +72,16 @@ pub enum JeffreysError { /// The number of trials exceeded the verified numerical regime. TrialsExceedSupported { /// Number of trials. - n: usize, + n: u64, /// Maximum supported number of trials. - max: usize, + max: u64, }, /// The number of successes exceeded the number of trials. SuccessesExceedTrials { /// Number of successes. - k: usize, + k: u64, /// Number of trials. - n: usize, + n: u64, }, /// The interval tail probability was outside `(0, 1)`. InvalidAlpha { @@ -136,7 +136,7 @@ impl From for JeffreysError { /// The posterior is `Beta(k + 1/2, n - k + 1/2)`. The interval endpoints are /// the `alpha/2` and `1 - alpha/2` posterior quantiles, with the PECOS endpoint /// rule `k = 0 -> lo = 0` and `k = n -> hi = 1`. The endpoint rule is applied -/// at this caller level; non-endpoint limits use [`special::invbetai`], whose +/// at this caller level; non-endpoint limits use [`special::betainc_inv`], whose /// reflected complement solve avoids losing the directly solved tail to /// arithmetic `1 - x` cancellation. /// @@ -160,11 +160,7 @@ impl From for JeffreysError { /// assert!((interval.point - 0.5).abs() < 1e-15); /// assert!(interval.lo < interval.point && interval.point < interval.hi); /// ``` -pub fn jeffreys_interval( - k: usize, - n: usize, - alpha: f64, -) -> Result { +pub fn jeffreys_interval(k: u64, n: u64, alpha: f64) -> Result { validate_jeffreys_inputs(k, n)?; validate_alpha(alpha)?; @@ -175,12 +171,12 @@ pub fn jeffreys_interval( let lo = if k == 0 { 0.0 } else { - special::invbetai(tail, a, b)?.x + special::betainc_inv(tail, a, b)?.x }; let hi = if k == n { 1.0 } else { - special::invbetai(1.0 - tail, a, b)?.x + special::betainc_inv(1.0 - tail, a, b)?.x }; let point = jeffreys_point(k, n, JeffreysEstimator::Mean)?; validate_interval_postcondition(lo, hi)?; @@ -208,11 +204,7 @@ pub fn jeffreys_interval( /// let point = jeffreys_point(3, 10, JeffreysEstimator::Mean).unwrap(); /// assert_eq!(point, 3.5 / 11.0); /// ``` -pub fn jeffreys_point( - k: usize, - n: usize, - estimator: JeffreysEstimator, -) -> Result { +pub fn jeffreys_point(k: u64, n: u64, estimator: JeffreysEstimator) -> Result { validate_jeffreys_inputs(k, n)?; match estimator { @@ -220,12 +212,12 @@ pub fn jeffreys_point( JeffreysEstimator::Median => { let a = count_to_f64(k) + 0.5; let b = count_to_f64(n - k) + 0.5; - Ok(special::invbetai(0.5, a, b)?.x) + Ok(special::betainc_inv(0.5, a, b)?.x) } } } -fn validate_jeffreys_inputs(k: usize, n: usize) -> Result<(), JeffreysError> { +fn validate_jeffreys_inputs(k: u64, n: u64) -> Result<(), JeffreysError> { if n == 0 { return Err(JeffreysError::ZeroTrials); } @@ -258,7 +250,7 @@ fn validate_interval_postcondition(lo: f64, hi: f64) -> Result<(), JeffreysError } #[allow(clippy::cast_precision_loss)] // Jeffreys support is n <= 1e8, exactly representable as f64 -fn count_to_f64(value: usize) -> f64 { +fn count_to_f64(value: u64) -> f64 { value as f64 } diff --git a/crates/pecos-num/tests/jeffreys.rs b/crates/pecos-num/tests/jeffreys.rs index 3c2dad5a2..9579f320f 100644 --- a/crates/pecos-num/tests/jeffreys.rs +++ b/crates/pecos-num/tests/jeffreys.rs @@ -20,8 +20,8 @@ const U01_BITS: u64 = 0x3ff0_0000_0000_0000; #[derive(Debug, Clone, Copy)] struct FixtureRow { - k: usize, - n: usize, + k: u64, + n: u64, alpha: f64, lo: f64, hi: f64, @@ -31,8 +31,8 @@ struct FixtureRow { #[derive(Debug, Clone, Copy)] struct WorstError { relative_error: f64, - k: usize, - n: usize, + k: u64, + n: u64, alpha: f64, bound: &'static str, } @@ -40,8 +40,8 @@ struct WorstError { #[derive(Debug, Clone, Copy)] struct WorstResidual { residual: f64, - k: usize, - n: usize, + k: u64, + n: u64, alpha: f64, bound: &'static str, } @@ -101,9 +101,9 @@ fn fixture_table_matches_scipy() { } #[test] -fn forward_betai_near_one_saturation_matches_scipy() { +fn forward_betainc_reg_near_one_saturation_matches_scipy() { // Quantile fixtures never evaluate the forward CDF in the near-1 saturation - // band inside the asymptotic gate: invbetai converges in the tail that stays + // band inside the asymptotic gate: betainc_inv converges in the tail that stays // cleanly below 1, so a forward-side defect there is invisible to every // interval test (a boundary-snap ordering bug survived four review rounds // this way, collapsing I_x ~ 1 to 0.0). Oracle: scipy.special.betainc. @@ -115,18 +115,18 @@ fn forward_betai_near_one_saturation_matches_scipy() { (44.6, 20_000.0, 0.006, 0.999_999_999_999_999_3), ]; for (a, b, x, expected) in cases { - let value = special::betai(a, b, x) - .unwrap_or_else(|err| panic!("betai({a}, {b}, {x}) failed: {err}")); + let value = special::betainc_reg(a, b, x) + .unwrap_or_else(|err| panic!("betainc_reg({a}, {b}, {x}) failed: {err}")); let relative = (value - expected).abs() / expected; assert!( relative <= 1.0e-12, - "betai({a}, {b}, {x}) = {value}, expected {expected}, rel {relative}" + "betainc_reg({a}, {b}, {x}) = {value}, expected {expected}, rel {relative}" ); } } #[test] -fn forward_betai_residual_at_fixture_quantiles_is_bounded() { +fn forward_betainc_reg_residual_at_fixture_quantiles_is_bounded() { let mut worst = WorstResidual::default(); for row in fixture_rows() { @@ -146,12 +146,12 @@ fn forward_betai_residual_at_fixture_quantiles_is_bounded() { } eprintln!( - "max forward betai fixture residual: {} at k={}, n={}, alpha={}, bound={}", + "max forward betainc_reg fixture residual: {} at k={}, n={}, alpha={}, bound={}", worst.residual, worst.k, worst.n, worst.alpha, worst.bound ); assert!( worst.residual <= 5.0e-9, - "forward betai residual too large: {worst:?}" + "forward betainc_reg residual too large: {worst:?}" ); } @@ -275,13 +275,13 @@ fn large_b_floor_side_sweeps_are_ordered() { #[test] fn large_b_saturation_keeps_inverse_orientation() { - let cdf_at_half = special::betai(5.5, 2_999_995.5, 0.5).unwrap(); + let cdf_at_half = special::betainc_reg(5.5, 2_999_995.5, 0.5).unwrap(); assert!( cdf_at_half > 1.0 - 1.0e-15, "large-b half CDF should saturate high, got {cdf_at_half}" ); - let cdf_at_guard = special::betai(5.5, 2_999_995.5, 1.0e-3).unwrap(); + let cdf_at_guard = special::betainc_reg(5.5, 2_999_995.5, 1.0e-3).unwrap(); assert!( cdf_at_guard > 1.0 - 1.0e-15, "large-b guard CDF should saturate high, got {cdf_at_guard}" @@ -358,9 +358,9 @@ fn typed_errors_are_returned() { ..InverseBetaOptions::default() }; assert!(matches!( - special::invbetai_with_options(0.5, 2.0, 2.0, options), + special::betainc_inv_with_options(0.5, 2.0, 2.0, options), Err(SpecialError::MaxIterations { - function: "invbetai", + function: "betainc_inv", iterations: 0 }) )); @@ -374,8 +374,10 @@ fn deterministic_randomized_grid_matches_puruspe_dev_check() { let mut rng = Lcg::new(0x4a45_4646_5245_5953); for _ in 0..200 { - let n = rng.next_usize(1_000) + 2; - let k = rng.next_usize(n - 2) + 1; + let n_usize = rng.next_usize(1_000) + 2; + let k_usize = rng.next_usize(n_usize - 2) + 1; + let n = usize_to_u64(n_usize); + let k = usize_to_u64(k_usize); let alpha = 1.0e-4 + rng.next_unit_f64() * (0.25 - 1.0e-4); let interval = jeffreys_interval(k, n, alpha).unwrap(); let median = jeffreys_point(k, n, JeffreysEstimator::Median).unwrap(); @@ -459,9 +461,9 @@ fn check_forward_residual( x: f64, target: f64, ) { - let actual = special::betai(a, b, x).unwrap_or_else(|err| { + let actual = special::betainc_reg(a, b, x).unwrap_or_else(|err| { panic!( - "forward betai failed for {bound}, k={}, n={}, alpha={}: {err}", + "forward betainc_reg failed for {bound}, k={}, n={}, alpha={}: {err}", row.k, row.n, row.alpha ) }); @@ -502,14 +504,14 @@ fn assert_close(actual: f64, expected: f64, relative_tolerance: f64, label: &str } fn check_k_sweep_monotone(k_start: usize, k_end: usize, n: usize, alpha: f64) { - let first = jeffreys_interval(k_start, n, alpha).unwrap(); - assert_interval_ordered(first, k_start, n, alpha); + let first = jeffreys_interval(usize_to_u64(k_start), usize_to_u64(n), alpha).unwrap(); + assert_interval_ordered(first, usize_to_u64(k_start), usize_to_u64(n), alpha); let mut previous_lo = first.lo; let mut previous_hi = first.hi; for k in (k_start + 1)..=k_end { - let interval = jeffreys_interval(k, n, alpha).unwrap(); - assert_interval_ordered(interval, k, n, alpha); + let interval = jeffreys_interval(usize_to_u64(k), usize_to_u64(n), alpha).unwrap(); + assert_interval_ordered(interval, usize_to_u64(k), usize_to_u64(n), alpha); assert!( interval.lo + 1.0e-14 >= previous_lo, "lo not monotone across gate at k={k}, n={n}, alpha={alpha}" @@ -525,8 +527,8 @@ fn check_k_sweep_monotone(k_start: usize, k_end: usize, n: usize, alpha: f64) { fn assert_interval_ordered( interval: pecos_num::stats::JeffreysInterval, - k: usize, - n: usize, + k: u64, + n: u64, alpha: f64, ) { assert!( @@ -540,10 +542,14 @@ fn uses_large_shape_median_shortcut(a: f64, b: f64) -> bool { } #[allow(clippy::cast_precision_loss)] // Test counts are <= 1e8, exactly representable as f64 -fn usize_to_f64(value: usize) -> f64 { +fn usize_to_f64(value: u64) -> f64 { value as f64 } +fn usize_to_u64(value: usize) -> u64 { + u64::try_from(value).unwrap() +} + struct Lcg { state: u64, } diff --git a/crates/pecos/tests/neo_emission_test.rs b/crates/pecos/tests/neo_emission_test.rs index aedca3019..36224887b 100644 --- a/crates/pecos/tests/neo_emission_test.rs +++ b/crates/pecos/tests/neo_emission_test.rs @@ -127,19 +127,23 @@ fn emission_is_gate_removing_and_matches_engines() { let engines = engines_zero_count(); let neo = neo_zero_count(); let facade = neo_facade_zero_count(); - let engines_ci = jeffreys_interval(engines, SHOTS as u64, CONFIDENCE); - let neo_ci = jeffreys_interval(neo, SHOTS as u64, CONFIDENCE); - let facade_ci = jeffreys_interval(facade, SHOTS as u64, CONFIDENCE); + let alpha = 1.0 - CONFIDENCE; + let engines_ci = jeffreys_interval(engines, SHOTS as u64, alpha) + .expect("Jeffreys interval for engines k and SHOTS n"); + let neo_ci = jeffreys_interval(neo, SHOTS as u64, alpha) + .expect("Jeffreys interval for neo k and SHOTS n"); + let facade_ci = jeffreys_interval(facade, SHOTS as u64, alpha) + .expect("Jeffreys interval for facade k and SHOTS n"); println!( "emission: engines {engines}/{SHOTS} CI [{:.4}, {:.4}], neo-direct {neo}/{SHOTS} CI \ [{:.4}, {:.4}], neo-facade {facade}/{SHOTS} CI [{:.4}, {:.4}], gate-removing analytic \ {analytic:.4} (gate-preserving would be {:.4})", - engines_ci.0, - engines_ci.1, - neo_ci.0, - neo_ci.1, - facade_ci.0, - facade_ci.1, + engines_ci.lo, + engines_ci.hi, + neo_ci.lo, + neo_ci.hi, + facade_ci.lo, + facade_ci.hi, P1 * 2.0 / 3.0 ); @@ -153,17 +157,17 @@ fn emission_is_gate_removing_and_matches_engines() { ("neo-facade", facade_ci), ] { assert!( - ci.0 <= analytic && analytic <= ci.1, + ci.lo <= analytic && analytic <= ci.hi, "{name} P(0) excludes the gate-removing analytic {analytic}" ); } // And every pair of stacks agrees. assert!( - engines_ci.0 <= neo_ci.1 && neo_ci.0 <= engines_ci.1, + engines_ci.lo <= neo_ci.hi && neo_ci.lo <= engines_ci.hi, "engines and neo-direct emission rates disagree: {engines}/{SHOTS} vs {neo}/{SHOTS}" ); assert!( - engines_ci.0 <= facade_ci.1 && facade_ci.0 <= engines_ci.1, + engines_ci.lo <= facade_ci.hi && facade_ci.lo <= engines_ci.hi, "engines and neo-facade emission rates disagree: {engines}/{SHOTS} vs {facade}/{SHOTS}" ); } @@ -248,19 +252,23 @@ fn two_qubit_emission_is_gate_removing_and_matches_engines() { let engines = engines_2q_zero_count(); let neo = neo_2q_zero_count(); let facade = neo_facade_2q_zero_count(); - let engines_ci = jeffreys_interval(engines, SHOTS as u64, CONFIDENCE); - let neo_ci = jeffreys_interval(neo, SHOTS as u64, CONFIDENCE); - let facade_ci = jeffreys_interval(facade, SHOTS as u64, CONFIDENCE); + let alpha = 1.0 - CONFIDENCE; + let engines_ci = jeffreys_interval(engines, SHOTS as u64, alpha) + .expect("Jeffreys interval for engines k and SHOTS n"); + let neo_ci = jeffreys_interval(neo, SHOTS as u64, alpha) + .expect("Jeffreys interval for neo k and SHOTS n"); + let facade_ci = jeffreys_interval(facade, SHOTS as u64, alpha) + .expect("Jeffreys interval for facade k and SHOTS n"); println!( "2q emission: engines {engines}/{SHOTS} CI [{:.4}, {:.4}], neo-direct {neo}/{SHOTS} CI \ [{:.4}, {:.4}], neo-facade {facade}/{SHOTS} CI [{:.4}, {:.4}], gate-removing analytic \ {analytic:.4} (gate-preserving would be {:.4})", - engines_ci.0, - engines_ci.1, - neo_ci.0, - neo_ci.1, - facade_ci.0, - facade_ci.1, + engines_ci.lo, + engines_ci.hi, + neo_ci.lo, + neo_ci.hi, + facade_ci.lo, + facade_ci.hi, P2 * 8.0 / 15.0 ); @@ -270,16 +278,16 @@ fn two_qubit_emission_is_gate_removing_and_matches_engines() { ("neo-facade", facade_ci), ] { assert!( - ci.0 <= analytic && analytic <= ci.1, + ci.lo <= analytic && analytic <= ci.hi, "{name} P(q1=0) excludes the gate-removing analytic {analytic}" ); } assert!( - engines_ci.0 <= neo_ci.1 && neo_ci.0 <= engines_ci.1, + engines_ci.lo <= neo_ci.hi && neo_ci.lo <= engines_ci.hi, "engines and neo-direct 2q emission rates disagree: {engines}/{SHOTS} vs {neo}/{SHOTS}" ); assert!( - engines_ci.0 <= facade_ci.1 && facade_ci.0 <= engines_ci.1, + engines_ci.lo <= facade_ci.hi && facade_ci.lo <= engines_ci.hi, "engines and neo-facade 2q emission rates disagree: {engines}/{SHOTS} vs {facade}/{SHOTS}" ); } diff --git a/crates/pecos/tests/neo_equivalence_matrix_test.rs b/crates/pecos/tests/neo_equivalence_matrix_test.rs index 4ab1fc34f..cf0037d9a 100644 --- a/crates/pecos/tests/neo_equivalence_matrix_test.rs +++ b/crates/pecos/tests/neo_equivalence_matrix_test.rs @@ -254,27 +254,30 @@ fn check_cell(name: &str, qasm: &str, cell: &NoiseCell, targets: &[&str], analyt let engines = count_targets(&cell.run(qasm, SimStack::Engines), targets); let neo = count_targets(&cell.run(qasm, SimStack::Neo), targets); - let engines_ci = jeffreys_interval(engines, SHOTS as u64, CONFIDENCE); - let neo_ci = jeffreys_interval(neo, SHOTS as u64, CONFIDENCE); + let alpha = 1.0 - CONFIDENCE; + let engines_ci = jeffreys_interval(engines, SHOTS as u64, alpha) + .expect("Jeffreys interval for engines k and SHOTS n"); + let neo_ci = jeffreys_interval(neo, SHOTS as u64, alpha) + .expect("Jeffreys interval for neo k and SHOTS n"); println!( "{name}: engines {engines}/{SHOTS} CI [{:.5}, {:.5}], \ neo {neo}/{SHOTS} CI [{:.5}, {:.5}], analytic {analytic:?}", - engines_ci.0, engines_ci.1, neo_ci.0, neo_ci.1 + engines_ci.lo, engines_ci.hi, neo_ci.lo, neo_ci.hi ); assert!( - engines_ci.0 <= neo_ci.1 && neo_ci.0 <= engines_ci.1, + engines_ci.lo <= neo_ci.hi && neo_ci.lo <= engines_ci.hi, "{name}: stack rates are statistically incompatible: \ engines {engines}/{SHOTS} vs neo {neo}/{SHOTS}" ); if let Some(truth) = analytic { assert!( - engines_ci.0 <= truth && truth <= engines_ci.1, + engines_ci.lo <= truth && truth <= engines_ci.hi, "{name}: engines rate {engines}/{SHOTS} excludes the analytic value {truth}" ); assert!( - neo_ci.0 <= truth && truth <= neo_ci.1, + neo_ci.lo <= truth && truth <= neo_ci.hi, "{name}: neo rate {neo}/{SHOTS} excludes the analytic value {truth}" ); } diff --git a/crates/pecos/tests/neo_surface_ler_test.rs b/crates/pecos/tests/neo_surface_ler_test.rs index 419751109..946192136 100644 --- a/crates/pecos/tests/neo_surface_ler_test.rs +++ b/crates/pecos/tests/neo_surface_ler_test.rs @@ -398,11 +398,13 @@ fn surface_memory_ler_matches_across_stacks() { // noise by an independent 6-seed 120k-shot-per-stack run (engines // 517 vs neo 482, z = 1.11). let equivalence_confidence = 0.99999; + let equivalence_alpha = 1.0 - equivalence_confidence; // The suppression margin is smaller than the equivalence margin, so // it gets its own (still strict) confidence. Pooling the two stacks // for suppression is justified by that 120k-shot equivalence run, // not by this test's own (weaker) overlap check. let suppression_confidence = 0.99; + let suppression_alpha = 1.0 - suppression_confidence; let mut pooled_intervals = Vec::new(); for distance in [3, 5] { @@ -411,26 +413,31 @@ fn surface_memory_ler_matches_across_stacks() { let neo = run_stack(&experiment, SimStack::Neo, p, shots, 42); let (engines_errors, neo_errors) = decode_logical_errors(&experiment, p, &engines, &neo); - let engines_ci = jeffreys_interval(engines_errors, shots as u64, equivalence_confidence); - let neo_ci = jeffreys_interval(neo_errors, shots as u64, equivalence_confidence); + let engines_ci = jeffreys_interval(engines_errors, shots as u64, equivalence_alpha) + .expect("Jeffreys interval for engines_errors k and shots n"); + let neo_ci = jeffreys_interval(neo_errors, shots as u64, equivalence_alpha) + .expect("Jeffreys interval for neo_errors k and shots n"); println!( "d={distance}: engines {engines_errors}/{shots} LER CI [{:.5}, {:.5}], \ neo {neo_errors}/{shots} LER CI [{:.5}, {:.5}]", - engines_ci.0, engines_ci.1, neo_ci.0, neo_ci.1 + engines_ci.lo, engines_ci.hi, neo_ci.lo, neo_ci.hi ); assert!( - engines_ci.0 <= neo_ci.1 && neo_ci.0 <= engines_ci.1, + engines_ci.lo <= neo_ci.hi && neo_ci.lo <= engines_ci.hi, "d={distance}: stack LERs are statistically incompatible: \ engines {engines_errors}/{shots} vs neo {neo_errors}/{shots}" ); // With per-stack equivalence established, pool the stacks for the // suppression physics check (doubles the statistics). - pooled_intervals.push(jeffreys_interval( - engines_errors + neo_errors, - 2 * shots as u64, - suppression_confidence, - )); + pooled_intervals.push( + jeffreys_interval( + engines_errors + neo_errors, + 2 * shots as u64, + suppression_alpha, + ) + .expect("Jeffreys interval for pooled errors k and pooled shots n"), + ); } // Error suppression: the pooled d=5 interval must sit strictly below @@ -438,12 +445,12 @@ fn surface_memory_ler_matches_across_stacks() { let d3 = pooled_intervals[0]; let d5 = pooled_intervals[1]; assert!( - d5.1 < d3.0, + d5.hi < d3.lo, "d=5 LER must be suppressed below d=3: \ d5 CI [{:.5}, {:.5}] vs d3 CI [{:.5}, {:.5}]", - d5.0, - d5.1, - d3.0, - d3.1 + d5.lo, + d5.hi, + d3.lo, + d3.hi ); } diff --git a/crates/pecos/tests/neo_v6_example_sweep_test.rs b/crates/pecos/tests/neo_v6_example_sweep_test.rs index ba50b6f5f..c9e4d3b38 100644 --- a/crates/pecos/tests/neo_v6_example_sweep_test.rs +++ b/crates/pecos/tests/neo_v6_example_sweep_test.rs @@ -131,14 +131,17 @@ fn assert_cross_stack_rate( let neo = run(program, SimStack::Neo, SEED ^ 0xA5A5, noise); let (e_count, _) = rate_where(&engines, &pred); let (n_count, _) = rate_where(&neo, &pred); - let e_ci = jeffreys_interval(e_count, SHOTS as u64, CONFIDENCE); - let n_ci = jeffreys_interval(n_count, SHOTS as u64, CONFIDENCE); + let alpha = 1.0 - CONFIDENCE; + let e_ci = jeffreys_interval(e_count, SHOTS as u64, alpha) + .expect("Jeffreys interval for engines k and SHOTS n"); + let n_ci = jeffreys_interval(n_count, SHOTS as u64, alpha) + .expect("Jeffreys interval for neo k and SHOTS n"); println!( "V6 {name}: engines {e_count}/{SHOTS} CI [{:.4}, {:.4}], neo {n_count}/{SHOTS} CI [{:.4}, {:.4}]", - e_ci.0, e_ci.1, n_ci.0, n_ci.1 + e_ci.lo, e_ci.hi, n_ci.lo, n_ci.hi ); assert!( - e_ci.0 <= n_ci.1 && n_ci.0 <= e_ci.1, + e_ci.lo <= n_ci.hi && n_ci.lo <= e_ci.hi, "V6 {name}: stack rates are statistically incompatible: \ engines {e_count}/{SHOTS} vs neo {n_count}/{SHOTS}" ); From 06a32e4eeeee2f5d282153b4145f1ae30af9a2a0 Mon Sep 17 00:00:00 2001 From: Ciaran Ryan-Anderson Date: Mon, 6 Jul 2026 21:53:07 -0600 Subject: [PATCH 07/11] Rewrite CF and gamma-tail routines from DLMF/A&S sources, removing NR-derived structure --- crates/pecos-num/src/special.rs | 258 +++++++++++++++++++------------- 1 file changed, 156 insertions(+), 102 deletions(-) diff --git a/crates/pecos-num/src/special.rs b/crates/pecos-num/src/special.rs index 23f5ef5dd..2a2ddfdbf 100644 --- a/crates/pecos-num/src/special.rs +++ b/crates/pecos-num/src/special.rs @@ -14,9 +14,9 @@ //! Special functions used by PECOS numerical routines. //! -//! The incomplete-beta implementation follows the Numerical Recipes-style -//! `ln_gamma` / continued-fraction / regularized-beta / inverse-beta chain, -//! with the reflection branch selected before forming complements. +//! The incomplete-beta implementation uses DLMF/Abramowitz-Stegun continued +//! fractions, public asymptotic expansions, and bracketed inverse solves with +//! the reflection branch selected before forming complements. use std::fmt; @@ -55,7 +55,12 @@ const BETA_ASYMPTOTIC_NORMALIZER_MIN: f64 = 100_000.0; const BETA_MEDIAN_ASYMPTOTIC_MIN: f64 = 1_000_000.0; const GAMMA_HALF_INTEGER_MAX_X: f64 = 700.0; const GAMMA_HALF_INTEGER_MAX_STEPS: usize = 128; -const FPMIN: f64 = f64::MIN_POSITIVE / f64::EPSILON; +/// Positive floor used by the modified Lentz continued-fraction recurrence. +/// +/// Lentz's algorithm evaluates continued fractions through paired forward and +/// backward recurrences; Thompson-Barnett style implementations replace a +/// denominator whose magnitude underflows with a tiny nonzero value. +const LENTZ_TINY: f64 = f64::MIN_POSITIVE / f64::EPSILON; const HALF_LN_TWO_PI: f64 = 0.918_938_533_204_672_8; const LN_EPSILON: f64 = -36.043_653_389_117_15; const LN_MIN_POSITIVE: f64 = -708.396_418_532_264_1; @@ -959,16 +964,19 @@ fn regularized_gamma_half_integer_pq(a: f64, x: f64) -> Result, } fn regularized_gamma_p_series(a: f64, x: f64) -> Result { - let mut ap = a; - let mut delta = 1.0 / a; - let mut sum = delta; + // DLMF 8.7.1 gives the lower incomplete-gamma power series. After + // multiplying by exp(-x) x^a / Gamma(a), the summand recurrence is + // x / (a + n) times the previous term. + let mut shifted_shape = a; + let mut term = 1.0 / shifted_shape; + let mut series = term; for _iteration in 1..=GAMMA_MAX_ITERATIONS { - ap += 1.0; - delta *= x / ap; - sum += delta; - if delta.abs() <= sum.abs() * GAMMA_EPSILON { - let value = sum * libm::exp(-x + a * libm::log(x) - ln_gamma(a)?); + shifted_shape += 1.0; + term *= x / shifted_shape; + series += term; + if term.abs() <= series.abs() * GAMMA_EPSILON { + let value = series * libm::exp(-x + a * libm::log(x) - ln_gamma(a)?); ensure_finite("regularized gamma lower series", value)?; return Ok(value); } @@ -981,91 +989,127 @@ fn regularized_gamma_p_series(a: f64, x: f64) -> Result { } fn regularized_gamma_q_fraction(a: f64, x: f64) -> Result { - let mut b = x + 1.0 - a; - let mut c = 1.0 / FPMIN; - let mut d = 1.0 / b; - let mut h = d; - - for iteration in 1..=GAMMA_MAX_ITERATIONS { - let i = usize_to_f64(iteration); - let an = -i * (i - a); - b += 2.0; - d = an * d + b; - if d.abs() < FPMIN { - d = FPMIN; - } - c = b + an / c; - if c.abs() < FPMIN { - c = FPMIN; - } - d = 1.0 / d; - let delta = d * c; - h *= delta; - if (delta - 1.0).abs() <= GAMMA_EPSILON { - let value = libm::exp(-x + a * libm::log(x) - ln_gamma(a)?) * h; - ensure_finite("regularized gamma upper fraction", value)?; - return Ok(value); - } - } - - Err(SpecialError::MaxIterations { - function: "regularized_gamma_q_fraction", - iterations: GAMMA_MAX_ITERATIONS, - }) + // DLMF 8.9.2 / A&S 6.5.31 express Gamma(a, x) as + // exp(-x) x^a times the reciprocal of a continued fraction. The reciprocal + // is evaluated with the modified Lentz algorithm (Lentz 1976; Thompson and + // Barnett 1986 underflow guard). + let fraction = modified_lentz_reciprocal( + x + 1.0 - a, + LentzProblem { + max_terms: GAMMA_MAX_ITERATIONS, + reported_iterations: GAMMA_MAX_ITERATIONS, + epsilon: GAMMA_EPSILON, + quantity_name: "regularized gamma upper fraction", + function_name: "regularized_gamma_q_fraction", + convergence_checkpoint: |_| true, + numerator: |term_index| { + let n = usize_to_f64(term_index); + n * (a - n) + }, + denominator: |term_index| x + usize_to_f64(2 * term_index + 1) - a, + }, + )?; + let value = libm::exp(-x + a * libm::log(x) - ln_gamma(a)?) * fraction; + ensure_finite("regularized gamma upper fraction", value)?; + Ok(value) } fn betacf(a: f64, b: f64, x: f64) -> Result { - let qab = a + b; - let qap = a + 1.0; - let qam = a - 1.0; - let mut c = 1.0; - let mut d = 1.0 - qab * x / qap; - if d.abs() < FPMIN { - d = FPMIN; - } - d = 1.0 / d; - let mut h = d; - - for iteration in 1..=BETACF_MAX_ITERATIONS { - let m = usize_to_f64(iteration); - let m2 = 2.0 * m; - 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 *= d * c; + // DLMF 8.17.22 writes the regularized-beta factor as the reciprocal of + // 1 + d_1/(1 + d_2/(1 + ...)), with d_{2m} and d_{2m+1} given in DLMF + // 8.17.23. The continued fraction is evaluated by modified Lentz updates + // (Lentz 1976; Thompson and Barnett 1986 underflow guard). + let even_coefficient = |m: f64| m * (b - m) * x / ((a + 2.0 * m - 1.0) * (a + 2.0 * m)); + let odd_coefficient = + |m: f64| -((a + m) * (a + b + m) * x) / ((a + 2.0 * m) * (a + 2.0 * m + 1.0)); + modified_lentz_reciprocal( + 1.0, + LentzProblem { + max_terms: 2 * BETACF_MAX_ITERATIONS + 1, + reported_iterations: BETACF_MAX_ITERATIONS, + epsilon: BETACF_EPSILON, + quantity_name: "regularized beta continued fraction", + function_name: "betacf", + convergence_checkpoint: |term_index| term_index > 1 && term_index % 2 == 1, + numerator: |term_index| { + if term_index % 2 == 0 { + even_coefficient(usize_to_f64(term_index / 2)) + } else { + odd_coefficient(usize_to_f64((term_index - 1) / 2)) + } + }, + denominator: |_| 1.0, + }, + ) +} - 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 delta = d * c; - h *= delta; +struct LentzProblem { + max_terms: usize, + reported_iterations: usize, + epsilon: f64, + quantity_name: &'static str, + function_name: &'static str, + convergence_checkpoint: C, + numerator: N, + denominator: D, +} - ensure_finite("regularized beta continued fraction", h)?; - if (delta - 1.0).abs() <= BETACF_EPSILON { - return Ok(h); +fn modified_lentz_reciprocal( + initial_denominator: f64, + problem: LentzProblem, +) -> Result +where + N: Fn(usize) -> f64, + D: Fn(usize) -> f64, + C: Fn(usize) -> bool, +{ + let LentzProblem { + max_terms, + reported_iterations, + epsilon, + quantity_name, + function_name, + convergence_checkpoint, + numerator, + denominator, + } = problem; + let mut backward_denominator_inverse = 1.0 / lentz_nonzero(initial_denominator); + let mut forward_denominator = 1.0 / LENTZ_TINY; + let mut reciprocal = backward_denominator_inverse; + + for term_index in 1..=max_terms { + let partial_numerator = numerator(term_index); + let partial_denominator = denominator(term_index); + + let next_backward_denominator = + partial_denominator + partial_numerator * backward_denominator_inverse; + backward_denominator_inverse = 1.0 / lentz_nonzero(next_backward_denominator); + forward_denominator = + lentz_nonzero(partial_denominator + partial_numerator / forward_denominator); + + let correction = forward_denominator * backward_denominator_inverse; + reciprocal *= correction; + ensure_finite(quantity_name, reciprocal)?; + + if convergence_checkpoint(term_index) && (correction - 1.0).abs() <= epsilon { + return Ok(reciprocal); } } Err(SpecialError::MaxIterations { - function: "betacf", - iterations: BETACF_MAX_ITERATIONS, + function: function_name, + iterations: reported_iterations, }) } +fn lentz_nonzero(value: f64) -> f64 { + if value.abs() < LENTZ_TINY { + LENTZ_TINY + } else { + value + } +} + fn betainc_inv_residual(p: f64, a: f64, b: f64, x: f64) -> Result { let cdf = betainc_reg_pair(a, b, x)?; if p <= 0.5 { @@ -1177,27 +1221,37 @@ fn stirling_correction(x: f64) -> f64 { fn initial_betainc_inv_guess(p: f64, a: f64, b: f64) -> f64 { if a >= 1.0 && b >= 1.0 { - let pp = if p < 0.5 { p } else { 1.0 - p }; - let t = libm::sqrt(-2.0 * libm::log(pp)); - let mut x = (2.307_53 + t * 0.270_61) / (1.0 + t * (0.992_29 + t * 0.044_81)) - t; + // A&S 26.2.23 supplies the rational approximation to the normal + // quantile; A&S 26.5.22 maps that deviate to an incomplete-beta + // quantile seed before the bracketed Newton/bisection solve. + let smaller_tail = if p < 0.5 { p } else { 1.0 - p }; + let tail_radius = libm::sqrt(-2.0 * libm::log(smaller_tail)); + let normal_correction = (2.307_53 + tail_radius * 0.270_61) + / (1.0 + tail_radius * (0.992_29 + tail_radius * 0.044_81)); + let mut normal_deviate = normal_correction - tail_radius; if p < 0.5 { - x = -x; + normal_deviate = -normal_deviate; } - let al = (x * x - 3.0) / 6.0; - let h = 2.0 / (1.0 / (2.0 * a - 1.0) + 1.0 / (2.0 * b - 1.0)); - let w = x * libm::sqrt(al + h) / h - - (1.0 / (2.0 * b - 1.0) - 1.0 / (2.0 * a - 1.0)) * (al + 5.0 / 6.0 - 2.0 / (3.0 * h)); - interiorize(a / (a + b * libm::exp(2.0 * w))) + let normal_curvature = (normal_deviate * normal_deviate - 3.0) / 6.0; + let adjusted_shape = 2.0 / (1.0 / (2.0 * a - 1.0) + 1.0 / (2.0 * b - 1.0)); + let shape_skew = 1.0 / (2.0 * b - 1.0) - 1.0 / (2.0 * a - 1.0); + let logit_half_shift = normal_deviate * libm::sqrt(normal_curvature + adjusted_shape) + / adjusted_shape + - shape_skew * (normal_curvature + 5.0 / 6.0 - 2.0 / (3.0 * adjusted_shape)); + interiorize(a / (a + b * libm::exp(2.0 * logit_half_shift))) } else { - let lna = libm::log(a / (a + b)); - let lnb = libm::log(b / (a + b)); - let t = libm::exp(a * lna) / a; - let u = libm::exp(b * lnb) / b; - let w = t + u; - let x = if p < t / w { - libm::exp(libm::log(a * w * p) / a) + // A&S 26.5.22 also gives a small-shape branch based on the endpoint + // weights of the two beta tails. + let total_shape = a + b; + let left_log_share = libm::log(a / total_shape); + let right_log_share = libm::log(b / total_shape); + let left_endpoint_weight = libm::exp(a * left_log_share) / a; + let right_endpoint_weight = libm::exp(b * right_log_share) / b; + let endpoint_weight_sum = left_endpoint_weight + right_endpoint_weight; + let x = if p < left_endpoint_weight / endpoint_weight_sum { + libm::exp(libm::log(a * endpoint_weight_sum * p) / a) } else { - 1.0 - libm::exp(libm::log(b * w * (1.0 - p)) / b) + 1.0 - libm::exp(libm::log(b * endpoint_weight_sum * (1.0 - p)) / b) }; interiorize(x) } From ff6442872b1cf6052848355080edac5337ad9055 Mon Sep 17 00:00:00 2001 From: Ciaran Ryan-Anderson Date: Mon, 6 Jul 2026 22:42:53 -0600 Subject: [PATCH 08/11] Add mpmath quadrature-based second-oracle fixture table (257 randomized full-regime rows) --- .../generate_jeffreys_mpmath_fixtures.py | 142 ++++++++++ .../tests/fixtures/jeffreys_mpmath.csv | 258 ++++++++++++++++++ 2 files changed, 400 insertions(+) create mode 100644 crates/pecos-num/tests/fixtures/generate_jeffreys_mpmath_fixtures.py create mode 100644 crates/pecos-num/tests/fixtures/jeffreys_mpmath.csv diff --git a/crates/pecos-num/tests/fixtures/generate_jeffreys_mpmath_fixtures.py b/crates/pecos-num/tests/fixtures/generate_jeffreys_mpmath_fixtures.py new file mode 100644 index 000000000..85f2a6d20 --- /dev/null +++ b/crates/pecos-num/tests/fixtures/generate_jeffreys_mpmath_fixtures.py @@ -0,0 +1,142 @@ +# Copyright 2026 The PECOS Developers +# +# Licensed under the Apache License, Version 2.0 (the "License"); you may +# not use this file except in compliance with the License. You may obtain a +# copy of the License at +# +# https://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, WITHOUT +# WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the +# License for the specific language governing permissions and limitations +# under the License. +"""Offline generator for the mpmath second-oracle Jeffreys fixture table. + +Independent of the primary scipy oracle by construction: mpmath evaluates the +regularized incomplete beta in arbitrary precision (30 significant digits +here) rather than through the TOMS-708 lineage scipy and R share, so an +agreement between both oracles and pecos-num rules out a shared-lineage bug. +scipy is used ONLY to seed the root find; every returned quantile is verified +by an mpmath residual check at full precision. + +Run offline with: + + uv run --with mpmath --with scipy python generate_jeffreys_mpmath_fixtures.py \ + > jeffreys_mpmath.csv + +The case grid is a fixed-seed randomized sweep over the full supported regime +(n in [1, 1e8], k in [0, n], alpha in [1e-6, 0.5]) with structured k patterns +(endpoints, near-endpoints, typical-LER counts, midpoints) so the table keeps +covering the regime-gate seams between the hand-picked scipy rows. +""" + +from __future__ import annotations + +import random +import sys + +from mpmath import mp, mpf +from scipy.stats import beta as scipy_beta + +mp.dps = 30 +RESIDUAL_LIMIT = mpf("1e-24") + + +def log_density(a: mpf, b: mpf, x: mpf) -> mpf: + return (a - 1) * mp.log(x) + (b - 1) * mp.log1p(-x) - mp.log(mp.beta(a, b)) + + +def beta_cdf(a: mpf, b: mpf, x: mpf) -> mpf: + """Regularized incomplete beta by arbitrary-precision quadrature ONLY. + + mpmath's betainc was found to return silently inaccurate values (no + exception) for large-shape near-saturation inputs, which disqualifies it + as an oracle. Direct quadrature of the density -- the definition itself, + no series machinery -- with the interval split at the mode, integrating + whichever tail is smaller and complementing. + """ + if x <= 0: + return mpf(0) + if x >= 1: + return mpf(1) + density = lambda t: mp.exp(log_density(a, b, t)) # noqa: E731 + mode = (a - 1) / (a + b - 2) if a + b > 2 else mpf("0.5") + if not 0 < mode < 1: + mode = mpf("0.5") + if x <= mode: + return mp.quad(density, [mpf(0), x]) + return 1 - mp.quad(density, [x, mpf(1)]) + + +def quantile_lower_side(a: float, b: float, p: mpf, seed: float) -> mpf: + # Bracket-guarded Newton on the quadrature CDF: the float64 scipy seed is + # already ~1e-16 relative, so 2-3 quadratically-converging steps reach far + # below RESIDUAL_LIMIT with only a handful of expensive quadrature calls. + # The bracket (updated from every evaluation's sign) catches any bad step. + am, bm = mpf(a), mpf(b) + lo, hi = mpf(0), mpf(1) + x = mpf(seed) if 0.0 < seed < 1.0 else mpf("0.5") + residual = None + for _ in range(8): + f = beta_cdf(am, bm, x) - p + residual = abs(f) + if f < 0: + lo = x + else: + hi = x + if residual <= RESIDUAL_LIMIT: + return x + step = f * mp.exp(-log_density(am, bm, x)) + candidate = x - step + if not lo < candidate < hi: + candidate = (lo + hi) / 2 + x = candidate + msg = f"residual {residual} too large for a={a}, b={b}, p={p}" + raise RuntimeError(msg) + + +def quantile(a: float, b: float, p: float) -> mpf: + # Root-find on whichever side of the distribution is well conditioned: + # near x = 1 the CDF derivative vanishes, so solve the mirrored problem + # I_y(b, a) = 1 - p for the small complement y instead (same complement + # strategy the Rust implementation uses). + seed = float(scipy_beta.ppf(p, a, b)) + if seed > 0.5: + return 1 - quantile_lower_side(b, a, 1 - mpf(p), 1.0 - seed) + return quantile_lower_side(a, b, mpf(p), seed) + + +def emit(k: int, n: int, alpha: float) -> str: + a = k + 0.5 + b = n - k + 0.5 + lo = mpf(0) if k == 0 else quantile(a, b, alpha / 2) + hi = mpf(1) if k == n else quantile(a, b, 1 - alpha / 2) + median = quantile(a, b, 0.5) + return f"{k},{n},{alpha!r},{float(lo)!r},{float(hi)!r},{float(median)!r}" + + +def k_patterns(rng: random.Random, n: int) -> list[int]: + picks = {0, n, min(1, n), max(n - 1, 0), n // 2} + picks.add(max(1, round(n * 0.001))) # typical logical-error-rate count + picks.add(rng.randint(0, n)) + return sorted(picks) + + +def main() -> None: + rng = random.Random(20260707) + print("k,n,alpha,lo,hi,median") + seen: set[tuple[int, int, float]] = set() + for _ in range(40): + n = round(10 ** rng.uniform(0.0, 8.0)) + n = max(1, min(n, 100_000_000)) + alpha = 10 ** rng.uniform(-6.0, -0.301) # alpha in [1e-6, 0.5] + for k in k_patterns(rng, n): + case = (k, n, alpha) + if case not in seen: + seen.add(case) + print(emit(k, n, alpha)) + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/crates/pecos-num/tests/fixtures/jeffreys_mpmath.csv b/crates/pecos-num/tests/fixtures/jeffreys_mpmath.csv new file mode 100644 index 000000000..4448197fc --- /dev/null +++ b/crates/pecos-num/tests/fixtures/jeffreys_mpmath.csv @@ -0,0 +1,258 @@ +k,n,alpha,lo,hi,median +0,2070,0.3969201931407626,0.0,0.00039930932521924924,0.00010986871807782242 +1,2070,0.3969201931407626,0.0002412584897306801,0.0011250917701347122,0.0005713970061968474 +2,2070,0.3969201931407626,0.0005633563448333703,0.0017652323559002845,0.001050905802425459 +915,2070,0.3969201931407626,0.43280666675079144,0.4512967613308926,0.442038321544494 +1035,2070,0.3969201931407626,0.4906921991106,0.5093078008894,0.5 +2069,2070,0.3969201931407626,0.9988749082298652,0.9997587415102693,0.9994286029938032 +2070,2070,0.3969201931407626,0.9996006906747807,1.0,0.9998901312819222 +0,9957271,0.19959844815533892,0.0,1.3601863312583126e-07,2.2844432299822626e-08 +1,9957271,0.19959844815533892,2.9299819199683084e-08,3.1414112625492874e-07,1.1880633776223225e-07 +9957,9957271,0.19959844815533892,0.0009871964746304284,0.0010128924398332749,0.0009999894886178807 +1968940,9957271,0.19959844815533892,0.19757705801522435,0.19790086606574245,0.1977389287440001 +4978635,9957271,0.19959844815533892,0.4997967029928488,0.500203196578042,0.4999999497854399 +9957270,9957271,0.19959844815533892,0.9999996858588738,0.9999999707001808,0.9999998811936622 +9957271,9957271,0.19959844815533892,0.9999998639813669,1.0,0.9999999771555677 +0,8,1.4734465388988322e-06,0.0,0.7723121561761173,0.02718112317372858 +1,8,1.4734465388988322e-06,1.2736486289881826e-05,0.8688508601427605,0.14173175510135005 +4,8,1.4734465388988322e-06,0.017270779453141517,0.9827292205465866,0.5 +7,8,1.4734465388988322e-06,0.13114913985600374,0.9999872635137095,0.85826824489865 +8,8,1.4734465388988322e-06,0.22768784382203416,1.0,0.9728188768262714 +0,223012,0.0033004655317866146,0.0,2.2202226722714145e-05,1.0199802726220548e-06 +1,223012,0.0033004655317866146,7.622133263549388e-08,3.408767739253921e-05,5.304580608698525e-06 +223,223012,0.0033004655317866146,0.0008166537787287906,0.0012104736985892612,0.0010006924362060126 +111506,223012,0.0033004655317866146,0.49688904806614936,0.5031109519338507,0.5 +203492,223012,0.0033004655317866146,0.910701432903165,0.9142181994361465,0.9124704612706709 +223011,223012,0.0033004655317866146,0.9999659123226075,0.9999999237786673,0.9999946954193913 +223012,223012,0.0033004655317866146,0.9999777977732773,1.0,0.9999989800197274 +0,3680688,5.3481910003556164e-05,0.0,2.3958098167732246e-06,6.180045390450179e-08 +1,3680688,5.3481910003556164e-05,2.9385583100175433e-10,3.241002629927359e-06,3.214037246380003e-07 +3681,3680688,5.3481910003556164e-05,0.0009350373543192368,0.0010681725428010552,0.0010001299590330262 +1840344,3680688,5.3481910003556164e-05,0.4989471369735302,0.5010528630264698,0.5 +2851218,3680688,5.3481910003556164e-05,0.7737620295201527,0.775521648015834,0.7746426506273912 +3680687,3680688,5.3481910003556164e-05,0.9999967589973701,0.9999999997061442,0.9999996785962754 +3680688,3680688,5.3481910003556164e-05,0.9999976041901832,1.0,0.9999999381995461 +0,24,5.614822396505769e-05,0.0,0.3034389147489194,0.009335707701480714 +1,24,5.614822396505769e-05,4.7045907431966964e-05,0.39372002861349437,0.04859624914270972 +12,24,5.614822396505769e-05,0.1521763913320379,0.8478236086679366,0.5 +23,24,5.614822396505769e-05,0.6062799713864583,0.999952954092568,0.9514037508572902 +24,24,5.614822396505769e-05,0.6965610852510319,1.0,0.9906642922985193 +0,2330,0.19837131579936385,0.0,0.0005831555060343595,9.761060324728048e-05 +1,2330,0.19837131579936385,0.00012463987411495548,0.0013447450016684309,0.0005076453115603701 +2,2330,0.19837131579936385,0.00034417053447899744,0.001985477629007523,0.0009336543478298047 +1165,2330,0.19837131579936385,0.48668078059306474,0.5133192194069353,0.5 +1712,2330,0.19837131579936385,0.7228589408370034,0.7463796474792355,0.7347303579995603 +2329,2330,0.19837131579936385,0.9986552549983315,0.9998753601258851,0.9994923546884397 +2330,2330,0.19837131579936385,0.9994168444939656,1.0,0.9999023893967527 +0,18290344,0.09792212817783398,0.0,1.0597633772738552e-07,1.2436518801188228e-08 +1,18290344,0.09792212817783398,9.474909458322731e-09,2.1491025513171435e-07,6.467822147579503e-08 +18290,18290344,0.09792212817783398,0.0009878088252369955,0.0010122713839439542,0.0009999902863727448 +9145172,18290344,0.09792212817783398,0.4998065092940877,0.5001934907059122,0.5 +10011084,18290344,0.09792212817783398,0.5471499489580319,0.5475351917664504,0.547342575088754 +18290343,18290344,0.09792212817783398,0.9999997850897449,0.9999999905250906,0.9999999353217786 +18290344,18290344,0.09792212817783398,0.9999998940236623,1.0,0.9999999875634812 +0,1216574,0.001421809830909092,0.0,4.710148491735513e-06,1.8697435877376057e-07 +1,1216574,0.001421809830909092,7.946425968833852e-09,6.9817138823796145e-06,9.723918233680632e-07 +1217,1216574,0.001421809830909092,0.0009118108435162663,0.0010947281690237294,0.0010004868995984204 +435401,1216574,0.001421809830909092,0.35650540720416113,0.35927844409577314,0.35789113310792825 +608287,1216574,0.001421809830909092,0.4985538417446066,0.5014461582553934,0.5 +1216573,1216574,0.001421809830909092,0.9999930182861176,0.999999992053574,0.9999990276081766 +1216574,1216574,0.001421809830909092,0.9999952898515083,1.0,0.9999998130256412 +0,626446,0.03220710242334012,0.0,4.622518249783535e-06,3.631088385666032e-07 +1,626446,0.03220710242334012,1.2703295489145995e-07,8.22927649453018e-06,1.8884090525092133e-06 +626,626446,0.03220710242334012,0.0009164636139353029,0.0010875152232991598,0.0009995536167745369 +313223,626446,0.03220710242334012,0.498646954600168,0.501353045399832,0.5 +406410,626446,0.03220710242334012,0.6474624571096531,0.6500460128304311,0.6487549611857479 +626445,626446,0.03220710242334012,0.9999917707235054,0.9999998729670451,0.9999981115909475 +626446,626446,0.03220710242334012,0.9999953774817503,1.0,0.9999996368911614 +0,80941,0.0032841558856387867,0.0,6.122749842600285e-05,2.8102839026239463e-06 +1,80941,0.0032841558856387867,2.0931143466292e-07,9.39821547176645e-05,1.4615361202975466e-05 +81,80941,0.0032841558856387867,0.0007103528126644534,0.0013661526503128225,0.0010027869156278438 +40470,80941,0.0032841558856387867,0.49482737863093074,0.505160267621261,0.4999938226863903 +68030,80941,0.0032841558856387867,0.8366796553369801,0.8442465703503818,0.840487348845386 +80940,80941,0.0032841558856387867,0.9999060178452823,0.9999997906885654,0.999985384638797 +80941,80941,0.0032841558856387867,0.999938772501574,1.0,0.9999971897160974 +0,42,0.02021183803343473,0.0,0.07530655294443471,0.005369288028268642 +1,42,0.02021183803343473,0.0013842363573579057,0.12681414482155295,0.027938557333543106 +18,42,0.02021183803343473,0.26402907121141783,0.6058271033434074,0.42914102882209654 +21,42,0.02021183803343473,0.3274008261926561,0.672599173807344,0.5 +41,42,0.02021183803343473,0.8731858551784472,0.9986157636426422,0.9720614426664569 +42,42,0.02021183803343473,0.9246934470555654,1.0,0.9946307119717314 +0,1137,0.00011591085796346737,0.0,0.007083307587771449,0.00019999600624805185 +1,1137,0.00011591085796346737,1.593931494991295e-06,0.009737541544759268,0.0010401334273624958 +568,1137,0.00011591085796346737,0.4426056798404053,0.5565226980550997,0.49956037520696717 +711,1137,0.00011591085796346737,0.5690599550215338,0.6793524406145472,0.6252930649886069 +1136,1137,0.00011591085796346737,0.9902624584552415,0.999998406068505,0.9989598665726375 +1137,1137,0.00011591085796346737,0.9929166924122292,1.0,0.9998000039937519 +0,20,4.5553747638787294e-05,0.0,0.35774600361211933,0.011169169123591416 +1,20,4.5553747638787294e-05,4.921179989004273e-05,0.4582404287180797,0.058150568699694996 +7,20,4.5553747638787294e-05,0.056267042944696366,0.7744121362715336,0.35252799902895493 +10,20,4.5553747638787294e-05,0.1275502387079983,0.8724497612920267,0.5 +19,20,4.5553747638787294e-05,0.5417595712819722,0.99995078820011,0.941849431300305 +20,20,4.5553747638787294e-05,0.6422539963879353,1.0,0.9888308308764086 +0,22499970,0.00011224101305721022,0.0,3.6065070097994606e-07,1.0109711607728342e-08 +1,22499970,0.00011224101305721022,7.881906855086376e-11,4.958651159755181e-07,5.2577266735594567e-08 +22500,22499970,0.00011224101305721022,0.0009744917360256621,0.0010259668660500395,0.00100000872597643 +7002064,22499970,0.00011224101305721022,0.3108263442062398,0.3115803470626388,0.31120326217911615 +11249985,22499970,0.00011224101305721022,0.4995928584894377,0.5004071415105623,0.5 +22499969,22499970,0.00011224101305721022,0.999999504134884,0.9999999999211809,0.9999999474227332 +22499970,22499970,0.00011224101305721022,0.999999639349299,1.0,0.9999999898902884 +0,4443484,0.00026921137727446383,0.0,1.640172710378394e-06,5.1191405874247704e-08 +1,4443484,0.00026921137727446383,7.155308661520365e-10,2.3051464699274465e-06,2.66229573747088e-07 +4443,4443484,0.00026921137727446383,0.0009462926046890916,0.0010555517815157946,0.0009999285105562856 +2126683,4443484,0.00026921137727446383,0.47774377095202675,0.47947052145178076,0.4786071035995553 +2221742,4443484,0.00026921137727446383,0.4991358334061979,0.5008641665938022,0.5 +4443483,4443484,0.00026921137727446383,0.9999976948535301,0.9999999992844691,0.9999997337704263 +4443484,4443484,0.00026921137727446383,0.9999983598272896,1.0,0.9999999488085941 +0,11,5.213266723595424e-05,0.0,0.5437337534966938,0.020010728788547778 +1,11,5.213266723595424e-05,9.89527783820491e-05,0.6720967231277288,0.10427191127311598 +3,11,5.213266723595424e-05,0.010362214138808959,0.8233947798654095,0.27978210687174526 +5,11,5.213266723595424e-05,0.049825780732129225,0.9176980831258713,0.4559438967432282 +10,11,5.213266723595424e-05,0.3279032768722268,0.9999010472216179,0.895728088726884 +11,11,5.213266723595424e-05,0.45626624650325276,1.0,0.9799892712114522 +0,3,1.3182816441911038e-06,0.0,0.9756796216574851,0.06737827302740958 +1,3,1.3182816441911038e-06,3.3528051042242605e-05,0.997464073124095,0.3524522798771308 +2,3,1.3182816441911038e-06,0.0025359268758278265,0.999966471948956,0.6475477201228692 +3,3,1.3182816441911038e-06,0.024320378341987837,1.0,0.9326217269725904 +0,1033,9.594260525329422e-05,0.0,0.007965510996420375,0.00022012404293531723 +1,1033,9.594260525329422e-05,1.5465372700017193e-06,0.010901707973270511,0.0011448169897759581 +516,1033,9.594260525329422e-05,0.4390769190012608,0.5599647964855716,0.4995161291128267 +530,1033,9.594260525329422e-05,0.4525134476415912,0.5733602343428027,0.513064513952359 +1032,1033,9.594260525329422e-05,0.9890982920267287,0.99999845346273,0.998855183010224 +1033,1033,9.594260525329422e-05,0.9920344890035789,1.0,0.9997798759570646 +0,22133825,0.0001043274991916736,0.0,3.697463853068475e-07,1.0276949773557405e-08 +1,22133825,0.0001043274991916736,7.630861685971356e-11,5.075127953260113e-07,5.344701713038129e-08 +22134,22133825,0.0001043274991916736,0.0009741723347771821,0.0010263111028696963,0.0010000154213842467 +11066912,22133825,0.0001043274991916736,0.4995875884347463,0.5004123663855555,0.4999999774101407 +21912383,22133825,0.0001043274991916736,0.9899129994321398,0.9900771665864435,0.9899953052248377 +22133824,22133825,0.0001043274991916736,0.9999994924872047,0.9999999999236914,0.9999999465529829 +22133825,22133825,0.0001043274991916736,0.9999996302536147,1.0,0.9999999897230503 +0,35163,0.10568062939236542,0.0,5.330653612709491e-05,6.4688979493661335e-06 +1,35163,0.10568062939236542,5.2052976310652964e-06,0.00010936195010724602,3.3642620180711515e-05 +35,35163,0.10568062939236542,0.0007517729593464307,0.0012978992558359357,0.001000110774825093 +16162,35163,0.10568062939236542,0.45533341267015826,0.4639330834620926,0.4596312446735791 +17581,35163,0.10568062939236542,0.4956718628731078,0.504299699823663,0.499985780642717 +35162,35163,0.10568062939236542,0.9998906380498928,0.999994794702369,0.9999663573798193 +35163,35163,0.10568062939236542,0.9999466934638729,1.0,0.9999935311020506 +0,559,5.8541756833053755e-06,0.0,0.019357573650536773,0.0004066552198257964 +1,559,5.8541756833053755e-06,4.4279909383382216e-07,0.025134136387981447,0.0021149634078566474 +209,559,5.8541756833053755e-06,0.28513125078802287,0.46888445543676804,0.3739571677333679 +279,559,5.8541756833053755e-06,0.40421021571245785,0.5940452201088878,0.4991060791731307 +558,559,5.8541756833053755e-06,0.9748658636120114,0.9999995572009062,0.9978850365921433 +559,559,5.8541756833053755e-06,0.9806424263494565,1.0,0.9995933447801743 +0,1472,0.054871031110608766,0.0,0.0016502560605275124,0.0001544918569112563 +1,1472,0.054871031110608766,7.821961085443479e-05,0.0031017084324681038,0.0008034731926355901 +101,1472,0.054871031110608766,0.05677491872542851,0.08208669687169012,0.06871192213501943 +736,1472,0.054871031110608766,0.47499952202196966,0.5250004779780303,0.5 +1471,1472,0.054871031110608766,0.9968982915675318,0.9999217803891456,0.9991965268073644 +1472,1472,0.054871031110608766,0.9983497439394725,1.0,0.9998455081430887 +0,10345644,0.0003585053390110271,0.0,6.783991854090376e-07,2.1986857808215344e-08 +1,10345644,0.0003585053390110271,3.7208795596471613e-10,9.610680701377612e-07,1.1434637642075522e-07 +10346,10345644,0.0003585053390110271,0.0009653880156391955,0.0010355320657690482,0.001000050488420515 +4509780,10345644,0.0003585053390110271,0.43536084541691433,0.4364612597243273,0.43591099996897253 +5172822,10345644,0.0003585053390110271,0.49944521653776913,0.5005547834622308,0.5 +10345643,10345644,0.0003585053390110271,0.9999990389319299,0.9999999996279121,0.9999998856536236 +10345644,10345644,0.0003585053390110271,0.9999993216008146,1.0,0.9999999780131422 +0,52851349,0.046971666215438156,0.0,4.8553233895674975e-08,4.303924389788088e-09 +1,52851349,0.046971666215438156,1.9546946472435594e-09,8.973778677640784e-08,2.238328740755364e-08 +52851,52851349,0.046971666215438156,0.0009913845346563638,0.0010086582331725327,0.000999996543772384 +16690426,52851349,0.046971666215438156,0.3156724230309279,0.315926461834621,0.31579943326328697 +26425674,52851349,0.046971666215438156,0.4998633616057199,0.5001366194732882,0.4999999905395036 +52851348,52851349,0.046971666215438156,0.9999999102622132,0.9999999980453054,0.9999999776167126 +52851349,52851349,0.046971666215438156,0.9999999514467661,1.0,0.9999999956960756 +0,118,1.5629554702317878e-05,0.0,0.08102082369536284,0.001921767324490256 +1,118,1.5629554702317878e-05,4.044123704115144e-06,0.1060989026316705,0.009996354792883907 +59,118,1.5629554702317878e-05,0.3091588564389423,0.6908411435610994,0.5 +67,118,1.5629554702317878e-05,0.3716799727176115,0.7499767980416937,0.567604760623538 +117,118,1.5629554702317878e-05,0.8939010973683659,0.9999959558762959,0.9900036452071161 +118,118,1.5629554702317878e-05,0.9189791763046715,1.0,0.9980782326755098 +0,7,0.11249565934802262,0.0,0.22200062347306193,0.03086709132710212 +1,7,0.11249565934802262,0.027998344891999104,0.429351115586378,0.16100681027099867 +3,7,0.11249565934802262,0.18171444980937546,0.7102356121919254,0.4320462150014419 +6,7,0.11249565934802262,0.570648884413622,0.9720016551080009,0.8389931897290013 +7,7,0.11249565934802262,0.7779993765269381,1.0,0.9691329086728979 +0,375,0.020372195239141314,0.0,0.008758293633196102,0.0006059939355084158 +1,375,0.020372195239141314,0.00015513901082813203,0.014970177577743563,0.003151760802511062 +153,375,0.020372195239141314,0.35033641230528345,0.4675761258197174,0.4080818247849459 +187,375,0.020372195239141314,0.43904720759686583,0.5583138443877109,0.49866785248082507 +374,375,0.020372195239141314,0.9850298224222565,0.9998448609891719,0.9968482391974889 +375,375,0.020372195239141314,0.9912417063668039,1.0,0.9993940060644916 +0,19897878,5.0580060160921126e-05,0.0,4.4584121810976155e-07,1.1431782193139253e-08 +1,19897878,5.0580060160921126e-05,5.2371855070322566e-11,6.024343821488242e-07,5.9452918642122816e-08 +19898,19897878,5.0580060160921126e-05,0.0009715706710722318,0.001029007555376915,0.001000014490707035 +3645082,19897878,5.0580060160921126e-05,0.1828382029045712,0.18354112407094317,0.1831894891306253 +9948939,19897878,5.0580060160921126e-05,0.4995457073923192,0.5004542926076808,0.5 +19897877,19897878,5.0580060160921126e-05,0.9999993975656178,0.9999999999476281,0.9999999405470814 +19897878,19897878,5.0580060160921126e-05,0.999999554158782,1.0,0.9999999885682178 +0,2,0.07473602411848329,0.0,0.6129751655979074,0.0955258180378211 +1,2,0.07473602411848329,0.0798407560087024,0.9201592439912976,0.5 +2,2,0.07473602411848329,0.3870248344020924,1.0,0.9044741819621789 +0,5943,0.04751526720309638,0.0,0.0004299954021099238,3.827263843788397e-05 +1,5943,0.04751526720309638,1.7523369506377557e-05,0.0007956351293474894,0.0001990440801947965 +6,5943,0.04751526720309638,0.00041670012784247553,0.0020940783543653713,0.0010381171510131354 +1965,5943,0.04751526720309638,0.3186343890877606,0.34281595577969126,0.33065058984562845 +2971,5943,0.04751526720309638,0.48706569803508554,0.5127661203184048,0.4999158721260544 +5942,5943,0.04751526720309638,0.9992043648706526,0.9999824766304937,0.9998009559198052 +5943,5943,0.04751526720309638,0.99957000459789,1.0,0.9999617273615621 +0,97667,3.3137529950553763e-06,0.0,0.00011751840102849402,2.3290094342885665e-06 +1,97667,3.3137529950553763e-06,1.7334057430729832e-09,0.00015163935679802783,1.2112410597986867e-05 +98,97667,3.3137529950553763e-06,0.0006043260728924217,0.001552814584938868,0.0010051146529511457 +37738,97667,3.3137529950553763e-06,0.37916659633222904,0.39365688783886643,0.3863949734153705 +48833,97667,3.3137529950553763e-06,0.49255516481163397,0.5074345978617046,0.49999488058102004 +97666,97667,3.3137529950553763e-06,0.9998483606432018,0.9999999982665942,0.999987887589402 +97667,97667,3.3137529950553763e-06,0.9998824815989713,1.0,0.9999976709905657 +0,45122,0.026684875294311287,0.0,6.784814491461787e-05,5.041141304025965e-06 +1,45122,0.026684875294311287,1.5498675031287182e-06,0.00011877902512765646,2.6217321686102228e-05 +45,45122,0.026684875294311287,0.0007064437675512469,0.001367823869148258,0.0010009922079354783 +3717,45122,0.026684875294311287,0.07954170659692505,0.08527840092628332,0.08237975293136021 +22561,45122,0.026684875294311287,0.4947838468792362,0.5052161531207638,0.5 +45121,45122,0.026684875294311287,0.9998812209748723,0.9999984501324969,0.9999737826783139 +45122,45122,0.026684875294311287,0.9999321518550853,1.0,0.9999949588586959 +0,53345,0.48598002234235504,0.0,1.2776719194924754e-05,4.264067121003351e-06 +1,53345,0.48598002234235504,1.1091095322901004e-05,3.914881009623485e-05,2.2176012169982462e-05 +53,53345,0.48598002234235504,0.0009045248286096367,0.0010948412446411277,0.0009966577186544152 +26672,53345,0.48598002234235504,0.49848236643240756,0.501498887899112,0.4999906271089011 +46122,53345,0.48598002234235504,0.86356176400427,0.8656259941809091,0.8645960908337857 +53344,53345,0.48598002234235504,0.9999608511899037,0.9999889089046771,0.99997782398783 +53345,53345,0.48598002234235504,0.999987223280805,1.0,0.999995735932879 +0,54699,0.00013444745140189717,0.0,0.000145214752514732,4.1585165901096775e-06 +1,54699,0.00013444745140189717,3.657148344717303e-08,0.00020050725884269657,2.1627078313421694e-05 +55,54699,0.00013444745140189717,0.0005751989565080482,0.0016186078600860367,0.001008550221360756 +10385,54699,0.00013444745140189717,0.18351188481195183,0.19631653368480673,0.18985910861012872 +27349,54699,0.00013444745140189717,0.4918287178132975,0.5081530036762671,0.49999085912076985 +54698,54699,0.00013444745140189717,0.9997994927411573,0.9999999634285166,0.9999783729216866 +54699,54699,0.00013444745140189717,0.9998547852474853,1.0,0.9999958414834099 +0,24303,0.0021385356682974968,0.0,0.00022018515971202284,9.359536154456335e-06 +1,24303,0.0021385356682974968,5.228358008324183e-07,0.0003316863269439709,4.867589678668543e-05 +24,24303,0.0021385356682974968,0.0004960237802650745,0.0017494503163618892,0.0009944102823303835 +12151,24303,0.0021385356682974968,0.4901331476154492,0.5098257164040249,0.4999794266904465 +15216,24303,0.0021385356682974968,0.6165331511340129,0.6355892736000616,0.6260938142513485 +24302,24303,0.0021385356682974968,0.999668313673056,0.9999994771641991,0.9999513241032133 +24303,24303,0.0021385356682974968,0.999779814840288,1.0,0.9999906404638456 +0,6620122,0.001646295634853215,0.0,8.45021709487753e-07,3.436012192250203e-08 +1,6620122,0.001646295634853215,1.610873629446837e-09,1.2596343628565965e-06,1.7869563146483957e-07 +6620,6620122,0.001646295634853215,0.0009618372551117559,0.0010391718128738776,0.0010000066972048842 +2659534,6620122,0.001646295634853215,0.4011352598580345,0.40233473111157897,0.40173489744675955 +3310061,6620122,0.001646295634853215,0.4993883355566424,0.5006116644433576,0.5 +6620121,6620122,0.001646295634853215,0.9999987403656372,0.9999999983891263,0.9999998213043685 +6620122,6620122,0.001646295634853215,0.9999991549782905,1.0,0.999999965639878 +0,136980,1.8246993762078646e-06,0.0,8.798235383581782e-05,1.6605899222230063e-06 +1,136980,1.8246993762078646e-06,8.302873214203032e-10,0.00011261651141877527,8.636180413096275e-06 +137,136980,1.8246993762078646e-06,0.0006469738229723218,0.001466122848587735,0.001001361343609538 +68490,136980,1.8246993762078646e-06,0.4935536104693286,0.5064463895306632,0.5 +83400,136980,1.8246993762078646e-06,0.6025438900507313,0.6151274692911617,0.608847742132222 +136979,136980,1.8246993762078646e-06,0.999887383488581,0.9999999991697127,0.9999913638195869 +136980,136980,1.8246993762078646e-06,0.999912017646164,1.0,0.9999983394100778 +0,4977255,2.407569458691979e-06,0.0,2.367884316095264e-06,4.570153527126477e-08 +1,4977255,2.407569458691979e-06,2.7488893012963482e-11,3.0420439366368582e-06,2.3767857202855523e-07 +4977,4977255,2.407569458691979e-06,0.0009346509308002874,0.00106828604371996,0.0009999821864272668 +1208897,4977255,2.407569458691979e-06,0.24197861900566692,0.2437915083699815,0.2428842978117951 +2488627,4977255,2.407569458691979e-06,0.49894300808123254,0.5010567910054219,0.4999998995430279 +4977254,4977255,2.407569458691979e-06,0.9999969579560634,0.9999999999725111,0.999999762321428 +4977255,4977255,2.407569458691979e-06,0.9999976321156839,1.0,0.9999999542984648 +0,908,0.06942395328490615,0.0,0.002451867819236641,0.0002504153258301557 +1,908,0.06942395328490615,0.0001495595326292627,0.004739612253196876,0.0013023591825940455 +251,908,0.06942395328490615,0.2501335966518866,0.30397420885440674,0.27651382031719784 +454,908,0.06942395328490615,0.4699083537964791,0.5300916462035209,0.5 +907,908,0.06942395328490615,0.9952603877468031,0.9998504404673707,0.998697640817406 +908,908,0.06942395328490615,0.9975481321807633,1.0,0.9997495846741699 From 62e6d6af48302f9316ac665649ba66fadb187cec Mon Sep 17 00:00:00 2001 From: Ciaran Ryan-Anderson Date: Mon, 6 Jul 2026 22:48:29 -0600 Subject: [PATCH 09/11] Add mpmath second-oracle cross-check test with documented primary-oracle error floor --- crates/pecos-num/tests/jeffreys.rs | 79 ++++++++++++++++++++++++++++++ 1 file changed, 79 insertions(+) diff --git a/crates/pecos-num/tests/jeffreys.rs b/crates/pecos-num/tests/jeffreys.rs index 9579f320f..2d5a5e286 100644 --- a/crates/pecos-num/tests/jeffreys.rs +++ b/crates/pecos-num/tests/jeffreys.rs @@ -16,6 +16,7 @@ use pecos_num::special::{self, InverseBetaOptions, SpecialError}; use pecos_num::stats::{JeffreysError, JeffreysEstimator, jeffreys_interval, jeffreys_point}; const FIXTURE_CSV: &str = include_str!("fixtures/jeffreys_scipy.csv"); +const MPMATH_FIXTURE_CSV: &str = include_str!("fixtures/jeffreys_mpmath.csv"); const U01_BITS: u64 = 0x3ff0_0000_0000_0000; #[derive(Debug, Clone, Copy)] @@ -100,6 +101,45 @@ fn fixture_table_matches_scipy() { ); } +#[test] +fn second_oracle_mpmath_fixture_table_cross_check() { + // The implementation's contractual oracle of record is scipy (per the + // design note), and it matches scipy fixtures to ~4.5e-13. Cross-checking + // against arbitrary-precision truth revealed scipy itself deviates from + // truth by up to ~3.05e-12 in extreme corners (measured at k=0, n=136980, + // alpha=1.8e-6, hi bound). The mpmath test therefore asserts relative + // error <= max(design formula, 5.0e-12) per value: its purpose is catching + // structural/shared-lineage errors, not last-digit drift that the scipy + // table already pins tighter. The design formula is the same scale-aware + // max(1e-12, 5e-15/max(x, 1-x)) on relative x error. + let mut worst = WorstError::default(); + + for row in mpmath_fixture_rows() { + let interval = jeffreys_interval(row.k, row.n, row.alpha).unwrap_or_else(|err| { + panic!( + "interval failed for k={}, n={}, alpha={}: {err}", + row.k, row.n, row.alpha + ) + }); + check_mpmath_fixture_value(&mut worst, row, "lo", interval.lo, row.lo); + check_mpmath_fixture_value(&mut worst, row, "hi", interval.hi, row.hi); + + let median = + jeffreys_point(row.k, row.n, JeffreysEstimator::Median).unwrap_or_else(|err| { + panic!( + "median failed for k={}, n={}, alpha={}: {err}", + row.k, row.n, row.alpha + ) + }); + check_mpmath_fixture_value(&mut worst, row, "median", median, row.median); + } + + eprintln!( + "max mpmath fixture relative error: {} at k={}, n={}, alpha={}, bound={}", + worst.relative_error, worst.k, worst.n, worst.alpha, worst.bound + ); +} + #[test] fn forward_betainc_reg_near_one_saturation_matches_scipy() { // Quantile fixtures never evaluate the forward CDF in the near-1 saturation @@ -407,6 +447,14 @@ fn fixture_rows() -> Vec { FIXTURE_CSV.lines().skip(1).map(parse_fixture_row).collect() } +fn mpmath_fixture_rows() -> Vec { + MPMATH_FIXTURE_CSV + .lines() + .skip(1) + .map(parse_fixture_row) + .collect() +} + fn parse_fixture_row(line: &str) -> FixtureRow { let mut parts = line.split(','); let row = FixtureRow { @@ -421,6 +469,33 @@ fn parse_fixture_row(line: &str) -> FixtureRow { row } +fn check_mpmath_fixture_value( + worst: &mut WorstError, + row: FixtureRow, + bound: &'static str, + actual: f64, + expected: f64, +) { + let relative_error = relative_x_error(actual, expected); + let tolerance = mpmath_fixture_relative_tolerance(expected); + if relative_error > worst.relative_error { + *worst = WorstError { + relative_error, + k: row.k, + n: row.n, + alpha: row.alpha, + bound, + }; + } + assert!( + relative_error <= tolerance, + "{bound} mismatch for k={}, n={}, alpha={}: actual={actual}, expected={expected}, rel={relative_error}, tol={tolerance}", + row.k, + row.n, + row.alpha, + ); +} + fn check_fixture_value( worst: &mut WorstError, row: FixtureRow, @@ -483,6 +558,10 @@ fn fixture_relative_tolerance(expected: f64) -> f64 { 1.0e-12_f64.max(5.0e-15 / expected.max(1.0 - expected)) } +fn mpmath_fixture_relative_tolerance(expected: f64) -> f64 { + fixture_relative_tolerance(expected).max(5.0e-12) +} + fn relative_x_error(actual: f64, expected: f64) -> f64 { if expected.to_bits() == 0.0_f64.to_bits() { if actual.to_bits() == 0.0_f64.to_bits() { From bcd4a0fe85ee8ecbad6fb2f1bf36f45559fd113f Mon Sep 17 00:00:00 2001 From: Ciaran Ryan-Anderson Date: Mon, 6 Jul 2026 23:05:32 -0600 Subject: [PATCH 10/11] Apply final-review polish: spell out beta_continued_fraction, document tolerance semantics --- crates/pecos-num/src/special.rs | 16 ++++++++-------- crates/pecos-num/tests/jeffreys.rs | 7 ++++++- 2 files changed, 14 insertions(+), 9 deletions(-) diff --git a/crates/pecos-num/src/special.rs b/crates/pecos-num/src/special.rs index 2a2ddfdbf..3dfea6f13 100644 --- a/crates/pecos-num/src/special.rs +++ b/crates/pecos-num/src/special.rs @@ -20,8 +20,8 @@ use std::fmt; -const BETACF_MAX_ITERATIONS: usize = 10_000; -const BETACF_EPSILON: f64 = f64::EPSILON; +const BETA_CONTINUED_FRACTION_MAX_ITERATIONS: usize = 10_000; +const BETA_CONTINUED_FRACTION_EPSILON: f64 = f64::EPSILON; const BETA_POWER_SERIES_MAX_ITERATIONS: usize = 10_000; const BETA_POWER_SERIES_EPSILON: f64 = f64::EPSILON; const GAMMA_MAX_ITERATIONS: usize = 10_000; @@ -637,7 +637,7 @@ fn lower_beta_tail(a: f64, b: f64, x: f64, prefactor: f64) -> Result Result { Ok(value) } -fn betacf(a: f64, b: f64, x: f64) -> Result { +fn beta_continued_fraction(a: f64, b: f64, x: f64) -> Result { // DLMF 8.17.22 writes the regularized-beta factor as the reciprocal of // 1 + d_1/(1 + d_2/(1 + ...)), with d_{2m} and d_{2m+1} given in DLMF // 8.17.23. The continued fraction is evaluated by modified Lentz updates @@ -1025,11 +1025,11 @@ fn betacf(a: f64, b: f64, x: f64) -> Result { modified_lentz_reciprocal( 1.0, LentzProblem { - max_terms: 2 * BETACF_MAX_ITERATIONS + 1, - reported_iterations: BETACF_MAX_ITERATIONS, - epsilon: BETACF_EPSILON, + max_terms: 2 * BETA_CONTINUED_FRACTION_MAX_ITERATIONS + 1, + reported_iterations: BETA_CONTINUED_FRACTION_MAX_ITERATIONS, + epsilon: BETA_CONTINUED_FRACTION_EPSILON, quantity_name: "regularized beta continued fraction", - function_name: "betacf", + function_name: "beta_continued_fraction", convergence_checkpoint: |term_index| term_index > 1 && term_index % 2 == 1, numerator: |term_index| { if term_index % 2 == 0 { diff --git a/crates/pecos-num/tests/jeffreys.rs b/crates/pecos-num/tests/jeffreys.rs index 2d5a5e286..9d3e896ab 100644 --- a/crates/pecos-num/tests/jeffreys.rs +++ b/crates/pecos-num/tests/jeffreys.rs @@ -506,7 +506,12 @@ fn check_fixture_value( // The design tolerance is expressed as // max(1e-12, 5e-15 / max(p, 1-p)) on relative x error. The fixture columns // are x-values, so this test interprets p as x_oracle and uses - // max(x_oracle, 1 - x_oracle). + // max(x_oracle, 1 - x_oracle). Note the second term never binds: with + // max(x, 1-x) >= 0.5 it is at most 1e-14, so the effective tolerance is a + // flat 1e-12 everywhere; the formula is kept verbatim from the design note. + // Property tests elsewhere use ~1e-9 tolerances because they compound two + // independent inverse solves plus a subtraction; these fixture rows are + // the accuracy pins. let relative_error = relative_x_error(actual, expected); let tolerance = fixture_relative_tolerance(expected); if relative_error > worst.relative_error { From 627a019d46f1eaa197566082735b33e02ad25235 Mon Sep 17 00:00:00 2001 From: Ciaran Ryan-Anderson Date: Mon, 6 Jul 2026 23:22:03 -0600 Subject: [PATCH 11/11] Fix CI: lint the fixture generators, bump crossbeam-epoch past RUSTSEC-2026-0204 --- Cargo.lock | 4 ++-- .../fixtures/generate_jeffreys_fixtures.py | 24 ++++++++----------- .../generate_jeffreys_mpmath_fixtures.py | 6 +++++ ruff.toml | 7 ++++++ 4 files changed, 25 insertions(+), 16 deletions(-) diff --git a/Cargo.lock b/Cargo.lock index 8102b0bca..c3f0b8868 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -1355,9 +1355,9 @@ dependencies = [ [[package]] name = "crossbeam-epoch" -version = "0.9.18" +version = "0.9.20" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "5b82ac4a3c2ca9c3460964f020e1402edd5753411d7737aa39c3714ad1b5420e" +checksum = "2d6914041f254d6e9176c01941b21115dcfb7089e55135a35411081bd106ef3f" dependencies = [ "crossbeam-utils", ] diff --git a/crates/pecos-num/tests/fixtures/generate_jeffreys_fixtures.py b/crates/pecos-num/tests/fixtures/generate_jeffreys_fixtures.py index e369b1ad8..f5e883085 100644 --- a/crates/pecos-num/tests/fixtures/generate_jeffreys_fixtures.py +++ b/crates/pecos-num/tests/fixtures/generate_jeffreys_fixtures.py @@ -32,6 +32,7 @@ def jeffreys_row(k: int, n: int, alpha: float) -> str: + """Format one CSV row of scipy-oracle interval bounds and median.""" a = k + 0.5 b = n - k + 0.5 lo = 0.0 if k == 0 else float(beta.ppf(alpha / 2, a, b)) @@ -41,11 +42,11 @@ def jeffreys_row(k: int, n: int, alpha: float) -> str: def cases() -> list[tuple[int, int, float]]: + """Return the adversarial (k, n, alpha) case grid, deduplicated in order.""" out: list[tuple[int, int, float]] = [] # Minimum and smallest non-trivial experiments for n in (1, 2): - for k in range(n + 1): - out.append((k, n, 0.05)) + out.extend((k, n, 0.05) for k in range(n + 1)) # Endpoint boundaries k=0 and k=n across the full supported regime for n in (100, 1_000, 10_000, 100_000, 1_000_000, 10_000_000, 100_000_000): for alpha in (0.01, 0.05, 0.1): @@ -58,11 +59,9 @@ def cases() -> list[tuple[int, int, float]]: # n=1e8 boundary band (v3 addition), all in-regime alphas n8 = 100_000_000 for k in (0, 1, 10, n8 - 10, n8 - 1, n8): - for alpha in (0.01, 0.05, 0.1, 1e-6): - out.append((k, n8, alpha)) + out.extend((k, n8, alpha) for alpha in (0.01, 0.05, 0.1, 1e-6)) # Typical logical-error-rate working regime - for n in (1_000, 10_000, 100_000, 1_000_000): - out.append((max(1, n // 1000), n, 0.05)) + out.extend((max(1, n // 1000), n, 0.05) for n in (1_000, 10_000, 100_000, 1_000_000)) # Wide-CI alpha extreme and best-conditioned symmetric peak for n in (1_000, 1_000_000): out.append((n // 2, n, 1e-6)) @@ -70,22 +69,19 @@ def cases() -> list[tuple[int, int, float]]: # Asymmetric Gauss-Legendre-region shapes (both Beta shapes > 3000), including # the k=3016, n=1e6 tail-orientation regression and a large-n asymmetric case for k in (3000, 3016, 3100, 5000, 10_000): - for alpha in (0.01, 0.05, 1e-6): - out.append((k, 1_000_000, alpha)) + out.extend((k, 1_000_000, alpha) for alpha in (0.01, 0.05, 1e-6)) out.append((15_881, 20_000, 0.05)) out.append((10_000, 100_000_000, 0.05)) out.append((67_867_393, 100_000_000, 0.05)) # Continued-fraction accuracy band: moderate a with huge b, straddling the # asymptotic-branch gate at a = 64.5 for k in (63, 64, 65, 100, 300, 1_000, 3_000): - for alpha in (0.01, 0.05): - out.append((k, 100_000_000, alpha)) + out.extend((k, 100_000_000, alpha) for alpha in (0.01, 0.05)) # Mid-band b in [1e7, 5e7): small-a upper bounds routed through the CF path # missed the first asymptotic gate (round-2 review); mirrors pin the lo side for n in (10_000_000, 20_000_000, 40_000_000, 49_999_990): for k in (10, 64, 300): - for alpha in (0.01, 0.05): - out.append((k, n, alpha)) + out.extend((k, n, alpha) for alpha in (0.01, 0.05)) out.append((n - 64, n, 0.05)) # Median-shortcut gate edges (min posterior shape near 1e5 and 1e6) out.append((100_000, 1_000_000, 0.05)) @@ -95,8 +91,7 @@ def cases() -> list[tuple[int, int, float]]: # CF served hi endpoints here at up to 1.1e-10; rows span the candidate # boundary region so both sides of wherever the final gate lands are pinned for n in (1_000_000, 3_000_000, 6_000_000, 8_900_002, 8_999_990): - for k in (2, 5, 20, 64): - out.append((k, n, 0.05)) + out.extend((k, n, 0.05) for k in (2, 5, 20, 64)) out.append((n - 2, n, 0.05)) # Deduplicate, preserving order seen: set[tuple[int, int, float]] = set() @@ -109,6 +104,7 @@ def cases() -> list[tuple[int, int, float]]: def main() -> None: + """Print the full oracle CSV to stdout.""" print("k,n,alpha,lo,hi,median") for k, n, alpha in cases(): print(jeffreys_row(k, n, alpha)) diff --git a/crates/pecos-num/tests/fixtures/generate_jeffreys_mpmath_fixtures.py b/crates/pecos-num/tests/fixtures/generate_jeffreys_mpmath_fixtures.py index 85f2a6d20..4927f4329 100644 --- a/crates/pecos-num/tests/fixtures/generate_jeffreys_mpmath_fixtures.py +++ b/crates/pecos-num/tests/fixtures/generate_jeffreys_mpmath_fixtures.py @@ -44,6 +44,7 @@ def log_density(a: mpf, b: mpf, x: mpf) -> mpf: + """Log of the Beta(a, b) density at x, in arbitrary precision.""" return (a - 1) * mp.log(x) + (b - 1) * mp.log1p(-x) - mp.log(mp.beta(a, b)) @@ -70,6 +71,7 @@ def beta_cdf(a: mpf, b: mpf, x: mpf) -> mpf: def quantile_lower_side(a: float, b: float, p: mpf, seed: float) -> mpf: + """Solve I_x(a, b) = p for x on the well-conditioned lower side.""" # Bracket-guarded Newton on the quadrature CDF: the float64 scipy seed is # already ~1e-16 relative, so 2-3 quadratically-converging steps reach far # below RESIDUAL_LIMIT with only a handful of expensive quadrature calls. @@ -97,6 +99,7 @@ def quantile_lower_side(a: float, b: float, p: mpf, seed: float) -> mpf: def quantile(a: float, b: float, p: float) -> mpf: + """Beta(a, b) quantile at p, solved on whichever side is well conditioned.""" # Root-find on whichever side of the distribution is well conditioned: # near x = 1 the CDF derivative vanishes, so solve the mirrored problem # I_y(b, a) = 1 - p for the small complement y instead (same complement @@ -108,6 +111,7 @@ def quantile(a: float, b: float, p: float) -> mpf: def emit(k: int, n: int, alpha: float) -> str: + """Format one CSV row of mpmath-oracle interval bounds and median.""" a = k + 0.5 b = n - k + 0.5 lo = mpf(0) if k == 0 else quantile(a, b, alpha / 2) @@ -117,6 +121,7 @@ def emit(k: int, n: int, alpha: float) -> str: def k_patterns(rng: random.Random, n: int) -> list[int]: + """Return the structured + random success counts probed for one n.""" picks = {0, n, min(1, n), max(n - 1, 0), n // 2} picks.add(max(1, round(n * 0.001))) # typical logical-error-rate count picks.add(rng.randint(0, n)) @@ -124,6 +129,7 @@ def k_patterns(rng: random.Random, n: int) -> list[int]: def main() -> None: + """Print the randomized second-oracle CSV to stdout.""" rng = random.Random(20260707) print("k,n,alpha,lo,hi,median") seen: set[tuple[int, int, float]] = set() diff --git a/ruff.toml b/ruff.toml index 711cbdb60..909707189 100644 --- a/ruff.toml +++ b/ruff.toml @@ -102,6 +102,13 @@ ignore = [ [lint.per-file-ignores] "**/__init__.py" = ["F401"] # imported but unused - Expected for __init__.py re-exports +# Standalone offline oracle-fixture generators inside Rust crate test dirs - +# scripts run via uv, not an importable package +"crates/*/tests/fixtures/*.py" = [ + "INP001", # Implicit namespace package - standalone scripts, not a package + "S311", # Standard pseudo-random - fixed-seed case sampling, not crypto +] + # Main pecos __init__.py - special case for module initialization "python/quantum-pecos/src/pecos/__init__.py" = [ "A004", # Import shadows builtin - intentional numpy-like API (abs, all, any, max, min, sum, round)