Fix inverse CDF tail refinement and Gamma regression

[ChrisRackauckas-Claude]

Please ignore this PR until it has been reviewed by @ChrisRackauckas.

Summary

• retain the Chebyshev nodes that bracket each CDF tail cutoff, instead of selecting the first nodes inside the retained mass
• preserve the existing public sampling API while allowing the package's internal inverse-CDF path to consume an explicit probability grid
• replace the flaky Gamma sample-mean assertion with a deterministic conditional-CDF identity against Distributions.jl

Root cause

The statistical failure and the numerical defect were related but distinct.

For Gamma(2, 1) with 5,000 samples, the sample-mean standard error is sqrt(2 / 5000) = 0.02. A 200-seed scan produced an empirical standard deviation of 0.02159, and 4/200 means exceeded the test's fixed atol = 0.05. For example, seed 21 gives 1.9366639977, so the existing assertion can fail under ordinary Monte Carlo variation.

The deterministic CDF comparison also exposed a real refinement bug. The old lower boundary was 0.0547810463, which discards 0.00144679235 of the reference Gamma mass. The code selected the first Chebyshev node whose CDF was above the cutoff, rather than the preceding node that brackets it. The same off-by-one direction was present at the upper cutoff.

After retaining the bracketing nodes, the explicit 999-point conditional-CDF grid has a maximum error of 1.67e-14 with no custom tolerance.

An isolated git bisect run identified b059b431e821bfdc57ceb9d69d5db3233abb6f81 (Add inverse transform sampling using Chebyshev technology) as the first commit containing the failing behavior and assertion.

Investigation workflow

1. Reproduced the stochastic failure on untouched upstream/master at 2f6205d.
2. Checked the explicit seed-22 stream across Julia 1.10.0, 1.10.8, 1.10.10, 1.10.11, 1.11.9, and 1.12.6; all produced the same passing mean 1.9784477486. The previously recorded 2.052681 result therefore depended on different ambient RNG consumption, not a Julia patch-level RNG boundary.
3. Quantified mean variance at 1,000, 5,000, and 20,000 samples and compared it with the theoretical standard error.
4. Compared generated samples to the documented public Distributions.Gamma and cdf APIs, which revealed the dropped lower-tail mass.
5. Added the deterministic probability-grid regression first and observed the focused suite fail 15/16 before applying the boundary fix.
6. Applied the bracketing fix and reran focused, compatibility, Core, QA, and formatting checks.

Verification

• deterministic regression before fix: 15 passed, 1 failed
• focused inverse-transform suite after fix: 16/16
• focused suite repeated 10 times: 160/160
• SpecialFunctions 1.8.8 compatibility graph: 16/16
• current SpecialFunctions 2.8.0 graph: 16/16
• GROUP=Core Pkg.test(): 20,522/20,522
• GROUP=QA Pkg.test(): 17/17
• Runic check across all 43 Julia files: 43/43

[ChrisRackauckas-Claude]

Investigation scratchpad for commit d35f558:

• Started from untouched upstream/master@2f6205d; abandoned tolerance-widening commit b9a7c9b was not used.
• Reproduced the mean assertion failure deterministically with seed 21 (1.9366639977 versus 2.0 ± 0.05).
• Measured 200 seeds at n = 5000: empirical SD 0.0215893, theoretical SD 0.02, 4 failures at the old bound.
• Bisected the behavior to b059b43.
• Found the numerical cause independently via CDF identity: old refined support [0.0547810463, 19.9452189537] drops 0.00144679235 lower-tail mass.
• Added the 999-point deterministic regression before fixing it: 15/16 focused assertions passed, with the new identity failing as intended.
• After bracketing the cutoff nodes: maximum conditional-CDF error 1.67e-14; focused repeat, SpecialFunctions 1/2 graphs, Core, QA, and Runic all pass.

No public API was added or changed; the explicit-probability hook is package-internal and drives the same inverse-CDF implementation used by sampleInverseCDF.

[ChrisRackauckas-Claude]

Command-level verification appendix

Commit and base

• base: upstream/master@2f6205d15562dfa10ca7c5ea2a50b715ccd47102
• commit: d35f558e66b2c2ebf751470add3808a9680775f2
• required trailer present: Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
• post-push git fetch upstream && git log HEAD..upstream/master --oneline: no output
• post-push git status --short --branch: clean branch tracking its fork remote

Clean-base stochastic reproduction

julia --project=. --startup-file=no -e '
using PolyChaos, Random, Statistics, Test
Random.seed!(21)
samples = sampleMeasure(5000, GammaMeasure(2.0, 1.0); method = "inversecdf")
@test isapprox(mean(samples), 2.0; atol = 0.05)
'

Observed 1.9366639976972615; the assertion failed. The theoretical standard error is sqrt(2 / 5000) = 0.02.

Across 200 seeds:

n	empirical mean	empirical SD	theoretical SD	failures at ±0.05
1,000	2.0015875661	0.0455758415	0.0447213595	51
5,000	2.0037292935	0.0215893326	0.0200000000	4
20,000	2.0023447123	0.0100023346	0.0100000000	0

An independent git bisect run from ddefb884 through current master, using the seed-21 reproduction as the bad predicate and absence of the feature as good, reported b059b431e821bfdc57ceb9d69d5db3233abb6f81 as the first bad commit.

Deterministic regression before the fix

After adding the explicit 999-point probability-grid identity but before changing the cutoff indices:

julia --project=. --startup-file=no -e 'include("test/inverse_transform_sampling.jl")'

Result: 15 passed, 1 failed. At requested probability 0.001, the old algorithm returned a sample with conditional reference CDF 0.0024453456. The old refined domain was [0.0547810463, 19.9452189537]; the omitted lower mass was 0.00144679235.

Owner-public reference API evidence

The regression uses only Distributions.jl's documented public reference interface:

• Gamma(α, θ) is documented as shape/scale
• cdf(d, x) is documented as cumulative probability evaluation

Local ownership/export verification on the tested graph:

(gamma_exported = true, gamma_owner = Distributions,
 cdf_exported = true, cdf_owner = Distributions)

No dependency internal or re-export-only name is used.

After-fix commands and counts

# focused
julia --project=. --startup-file=no -e 'include("test/inverse_transform_sampling.jl")'
# 16/16

# focused, repeated in one process
julia --project=. --startup-file=no -e '
for iteration in 1:10
    include("test/inverse_transform_sampling.jl")
end
'
# 10 × 16/16

# package groups
GROUP=Core julia --project=. --startup-file=no -e 'using Pkg; Pkg.test()'
# 20,522/20,522

GROUP=QA julia --project=. --startup-file=no -e 'using Pkg; Pkg.test()'
# 17/17

# all repository Julia files
julia --project=/path/to/runic-env --startup-file=no \
    -e 'using Runic; exit(Runic.main(ARGS))' -- --check --verbose .
# 43/43

The focused suite also passed 16/16 in separate exact graphs with SpecialFunctions 1.8.8 and 2.8.0. The corrected 999-point grid's maximum conditional-CDF error is 1.6653345369e-14. git diff --check and the final staged diff check were clean; no tolerance, assertion, test runner, or warning policy was loosened.