Conversation
`(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.
Coverage Report
The above coverage report was generated for the changes in this PR. |
| if ( x < 0.0 ) { | ||
| return NINF; | ||
| } | ||
| lt = k * ln( x / lambda ); |
There was a problem hiding this comment.
| 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 ) ); |
There was a problem hiding this comment.
| 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>
|
@0PrashantYadav0 thanks, applied both in 28f76a4. |
`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.
Description
This pull request:
(x/lambda)^kfor smallxinstats/base/dists/weibull/logcdf(JS, factory and C).For small
x,t = (x/lambda)^kunderflows to0and the function returns-Infinity, although the log-CDF is finite there:ln(1 - e^{-t}) = ln(t) + ln(1 - t/2 + ...). Onceln(t) = k*ln(x/lambda)is below-36, the correction is under1.2e-16, which is less than half an ULP ofln(t). So the function returnsk*ln(x/lambda)there, which never underflows, and keeps the existinglog1p/expm1paths above it.x = 0still gives-Infinity.Fixtures and tests:
runner.jlnow evaluates the textbook formulaln(1 - exp(-(x/lambda)^k))in 2048-bitBigFloatinstead ofDistributions.logcdf. It's written withexpm1so it stays exact whentis far below even 2048 bits.Distributionsis dropped fromREQUIRE, and the three fixtures are regenerated with the same input ranges.small_x.jsonwithxlog-spaced over[1e-300, 0.01]. It passes at 1 ULP in JS and native. Ondevelop, 552 of the 4,021 assertions intest.logcdf.jsfail, most of them returning-Infinity.developand 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 onlarge_shapeat results like-3.2e-279, i.e.-exp(-t)witht ≈ 640. A half-ULP rounding intchangesexp(-t)byttimes 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
None. Related to #15762 (
weibull/quantile).Questions
No.
Other
No.
Checklist
AI Assistance
If you answered "yes" above, how did you use AI assistance?
Disclosure
This PR was written primarily by Claude Code. I found the bug by evaluating every
stats/base/distslogcdf just above the lower end of its support.@stdlib-js/reviewers