Skip to content

fix: avoid underflow for small x in stats/base/dists/weibull/logcdf - #15797

Open
Abhist17 wants to merge 3 commits into
stdlib-js:developfrom
Abhist17:fix/weibull-logcdf-small-x
Open

Abhist17 wants to merge 3 commits into
stdlib-js:developfrom
Abhist17:fix/weibull-logcdf-small-x

Conversation

@Abhist17

@Abhist17 Abhist17 commented Oct 2, 2026

Copy link
Copy Markdown
Contributor

Description

What is the purpose of this pull request?

This pull request:

  • avoids the underflow of (x/lambda)^k for small x in stats/base/dists/weibull/logcdf (JS, factory and C).

For small x, t = (x/lambda)^k underflows to 0 and the function returns -Infinity, although the log-CDF is finite there:

var logcdf = require( '@stdlib/stats/base/dists/weibull/logcdf' );

logcdf( 1.0e-300, 2.5, 1.5 );
// develop: -Infinity
// this PR: -1727.9524825158046

ln(1 - e^{-t}) = ln(t) + ln(1 - t/2 + ...). Once ln(t) = k*ln(x/lambda) is below -36, the correction is under 1.2e-16, which is less than half an ULP of ln(t). So the function returns k*ln(x/lambda) there, which never underflows, and keeps the existing log1p/expm1 paths above it. x = 0 still gives -Infinity.

Fixtures and tests:

  • runner.jl now evaluates the textbook formula ln(1 - exp(-(x/lambda)^k)) in 2048-bit BigFloat instead of Distributions.logcdf. It's written with expm1 so it stays exact when t is far below even 2048 bits. Distributions is dropped from REQUIRE, and the three fixtures are regenerated with the same input ranges.
  • There's a new small_x.json with x log-spaced over [1e-300, 0.01]. It passes at 1 ULP in JS and native. On develop, 552 of the 4,021 assertions in test.logcdf.js fail, most of them returning -Infinity.
  • Against the exact values, develop and this PR have the same error on the regenerated fixtures: 2, 1 and 3,236 ULP. The tolerances go from 0 to those values. The 3,236 is on large_shape at results like -3.2e-279, i.e. -exp(-t) with t ≈ 640. A half-ULP rounding in t changes exp(-t) by t times as much, so that's the conditioning of the function. The old 0 only held because Distributions computes the same double-precision expression.

Related Issues

Does this pull request have any related issues?

None. Related to #15762 (weibull/quantile).

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.

No.

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. I found the bug by evaluating every stats/base/dists logcdf just above the lower end of its support.


@stdlib-js/reviewers

`(x/lambda)^k` underflows to `0` for small `x`, so the function returned
`-Infinity` although the log-CDF is `k*ln(x/lambda)` to working
precision there. Return `k*ln(x/lambda)` once it is below `-36`. The
fixtures are now generated from the textbook formula in extended
precision.
@Abhist17
Abhist17 requested a review from a team October 2, 2026 08:13
@stdlib-bot stdlib-bot added Statistics Issue or pull request related to statistical functionality. Needs Review A pull request which needs code review. labels Oct 2, 2026
@stdlib-bot

stdlib-bot commented Oct 2, 2026 •

Copy link
Copy Markdown
Contributor

Coverage Report

Package Statements Branches Functions Lines
stats/base/dists/weibull/logcdf $\\color{green}337/337$
$\\color{green}+0.00\\%$
$\\color{green}38/38$
$\\color{green}+0.00\\%$
$\\color{green}4/4$
$\\color{green}+0.00\\%$
$\\color{green}337/337$
$\\color{green}+0.00\\%$

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

if ( x < 0.0 ) {
return NINF;
}
lt = k * ln( x / lambda );

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
lt = k * ln( x / lambda );
r = x / lambda;
if ( r < FLOAT64_SMALLEST_NORMAL ) {
lt = k * ( ln( x ) - ln( lambda ) );
} else {
lt = k * ln( r );
}

x / lambda can still underflow: logcdf( 5.0e-324, 2.5, 5.0 ) returns -Infinity (exact ≈ -1865.12). The difference is only used when the ratio is subnormal, because ln( x ) - ln( lambda ) cancels near x = lambda (~100 ULP at k = 10).

Same in lib/factory.js:82 and src/main.c:55. Needs var r; and @stdlib/constants/float64/smallest-normal (plus the manifest.json dep).

gen( x, k, lambda, "both_large.json" );

# Small `x`:
x = exp10.( range( -300.0, stop = -2.0, length = 1000 ) );

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
x = exp10.( range( -300.0, stop = -2.0, length = 1000 ) );
x = exp10.( range( -323.0, stop = -2.0, length = 1000 ) );

This covers a subnormal x / lambda. With the change above it still passes at 1 ULP.

Co-authored-by: 0PrashantYadav0 <0PrashantYadav0@users.noreply.github.com>
@Abhist17

Abhist17 commented Oct 3, 2026

Copy link
Copy Markdown
Contributor Author

@0PrashantYadav0 thanks, applied both in 28f76a4. r = x/lambda with the ln(x) - ln(lambda) fallback only for a subnormal ratio, in JS, factory and C (plus the smallest-normal dependency). small_x.json now goes down to 1e-323 and still passes at 1 ULP in JS and native. logcdf( 5.0e-324, 2.5, 5.0 ) now returns -1865.1237745845383.

`log1p( -exp( -Infinity ) )` is `log1p( -0.0 )`, so the function
returns `-0.0` there, as it does on develop; the doctest flagged the
example once the file changed.

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.

3 participants