Skip to content

fix: pass betainc the argument below one half in stats/base/dists/f/cdf - #15603

Open
Abhist17 wants to merge 2 commits into
stdlib-js:developfrom
Abhist17:fix/f-cdf-small-x
Open

Abhist17 wants to merge 2 commits into
stdlib-js:developfrom
Abhist17:fix/f-cdf-small-x

Conversation

@Abhist17

Copy link
Copy Markdown
Contributor

Description

What is the purpose of this pull request?

This pull request:

  • fixes stats/base/dists/f/cdf losing precision, down to returning 0, for small x.
var cdf = require( '@stdlib/stats/base/dists/f/cdf' );

cdf( 1.0e-10, 2.0, 3.0 );
// returns 1.0000000827237043e-10, exact is 9.999999999166667e-11 (8e-8 relative error)

cdf( 3.0e-19, 1.285, 4.949 );
// returns 0, exact is ~1e-12

The two branches hand betainc the argument that is close to one:

if ( d1 * x > d2 ) {
    return betainc( (d1*x)/(d2+(d1*x)), d1/2.0, d2/2.0, true, false ); // argument > 1/2
}
return betainc( d2/(d2+(d1*x)), d2/2.0, d1/2.0, true, true );          // argument > 1/2

For small x, d2/(d2+d1*x) is 1 - d1*x/d2 rounded, so every digit of d1*x/d2 below machine epsilon is lost before betainc runs, and the upper tail comes back as 0. This PR swaps the branches so each passes the argument below one half, which is exact to rounding. That is the same split as Boost's cdf( fisher_f_distribution ): ibetac( df2/2, df1/2, df2/(df2+v1x) ) when v1x > df2, and ibeta( df1/2, df2/2, v1x/(df2+v1x) ) otherwise. The change is to lib/main.js and lib/factory.js.

Checked against mpmath's regularized incomplete beta (40 digits) on 3,000 random (x, d1, d2) with x in [1e-17, 1e9] and d1, d2 in [0.1, 800], excluding results below 1e-290:

relative error > 1e-12 median error
develop 1,282 of 2,569 ~116,000 ULP
this PR 0 ~1 ULP

Related Issues

Does this pull request have any related issues?

This pull request has the following related issues:

Questions

Any questions for reviewers of this pull request?

No.

Other

Any other information relevant to this pull request? This may include screenshots, references, and/or implementation notes.

  • New fixture small_x.json (Julia, x = 10^U(-20, 0), d1, d2 in [0.5, 20.5]), added to runner.jl. The other fixture files are unchanged. On develop, 942 assertions in test.cdf.js fail, 200 of them because develop returns exactly 0.
  • The tolerance on the new fixture is 520 ULP. That is Julia's error, not stdlib's: against mpmath, the Julia values are off by up to 329 ULP (median 19), while this implementation is off by at most 20 ULP (median 2).
  • The existing tolerances drop from 58, 542 and 171 ULP to 10, 17 and 16.
  • Side note, not changed here: runner.jl assigns d2 twice in the both_large block and never sets d1, so both_large.json reuses d1 from the previous block.

Checklist

Please ensure the following tasks are completed before submitting this pull request.

AI Assistance

When authoring the changes proposed in this PR, did you use any kind of AI assistance?

  • Yes
  • No

If you answered "yes" above, how did you use AI assistance?

  • Code generation (e.g., when writing an implementation or fixing a bug)
  • Test/benchmark generation
  • Documentation (including examples)
  • Research and understanding

Disclosure

If you answered "yes" to using AI assistance, please provide a short disclosure indicating how you used AI assistance. This helps reviewers determine how much scrutiny to apply when reviewing your contribution. Example disclosures: "This PR was written primarily by Claude Code." or "I consulted ChatGPT to understand the codebase, but the proposed changes were fully authored manually by myself.".

This PR was written primarily by Claude Code. The bug was found by comparing every stats/base/dists function against SciPy on random parameters.


@stdlib-js/reviewers

…f/cdf`

Both branches handed `betainc` the argument close to one, so it was
rounded before the incomplete beta saw it. For small `x`,
`d2/(d2+d1*x)` rounds to `1` and the upper tail returns `0`, and cdf
values below about `1e-8` lost most of their digits. Swap the branches
so each uses the argument below one half, as Boost does, and add a
fixture for `x` down to `1e-20`.
@Abhist17
Abhist17 requested a review from a team September 27, 2026 22:51
@stdlib-bot stdlib-bot added Statistics Issue or pull request related to statistical functionality. Needs Review A pull request which needs code review. labels Sep 27, 2026
@stdlib-bot

stdlib-bot commented Sep 27, 2026 •

Copy link
Copy Markdown
Contributor

Coverage Report

Package Statements Branches Functions Lines
stats/base/dists/f/cdf $\\color{green}263/263$
$\\color{green}+0.00\\%$
$\\color{green}31/31$
$\\color{green}+0.00\\%$
$\\color{green}3/3$
$\\color{green}+0.00\\%$
$\\color{green}263/263$
$\\color{green}+0.00\\%$

The above coverage report was generated for the changes in this PR.

This branch has not been deployed

No deployments
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

Needs Review A pull request which needs code review. Statistics Issue or pull request related to statistical functionality.

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants