fix: evaluate stats/base/dists/normal/quantile via erfcinv to keep the lower tail

[Abhist17]

Description

What is the purpose of this pull request?

This pull request:

• fixes stats/base/dists/normal/quantile (JS, factory, and C) returning -Infinity, or losing digits, in the lower tail. lognormal/quantile, which calls it, inherits the fix.

var quantile = require( '@stdlib/stats/base/dists/normal/quantile' );
var lognormalQuantile = require( '@stdlib/stats/base/dists/lognormal/quantile' );

quantile( 1.0e-20, 0.0, 1.0 );          // -Infinity (expected ~-9.2623)
quantile( 1.0e-300, 0.0, 1.0 );         // -Infinity (expected ~-37.047)
lognormalQuantile( 1.0e-20, 0.0, 1.0 ); // 0         (expected ~9.49e-5)

The implementation evaluates mu + sigma*sqrt(2)*erfinv(2p-1). Forming 2p-1 throws away the low-order digits of small p, and once p < 2^-54 the result is exactly -1, so erfinv returns -Infinity. Before that point, the lower tail already loses accuracy (relative error ~1e-16/p in the probability).

The fix evaluates the equivalent mu - sigma*sqrt(2)*erfcinv(2p), which is what binomial/quantile and poisson/quantile already do for their normal approximations (-erfcinv( 2.0*p ) * SQRT2). 2p is exact, and erfcinv of a small argument is accurate. The upper tail isn't hurt: there, erfcinv receives an argument near 2, which it handles as 2 - y without cancellation. The method is documented in the JSDoc and C doc comments.

Related Issues

Does this pull request have any related issues?

No.

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.

• Fixtures: a new Julia runner evaluates quantile( Normal( mu, sigma ), p ) from Distributions.jl in 2048-bit BigFloat (an independent reference, not the formula used here) for 500 points in each tail. lower_tail.json has p log-uniform in [1e-250, 1e-3], and upper_tail.json has 1 - p log-uniform in [1e-16, 1e-3]. mu is in [-10, 10] and sigma in [0.1, 10.1].
  • On develop, 472 of the 500 lower-tail points return -Infinity, and 499 assertions fail in test.quantile.js.
  • With this change, the lower tail is within 13 ULP (median 1), and the tests allow 16. The upper tail is within 90 ULP, where develop is at 96. That error is the cancellation in mu + sigma*z when the result is near zero, not the tail evaluation, and the tests allow 96. The existing R fixture test is unchanged and passes.
• The lower-tail fixture stops at 1e-250 because math/base/special/erfcinv itself loses accuracy below about 1e-260 (up to ~500 ULP near 1e-300, against 2048-bit erfcinv). That's a separate issue in erfcinv, which I'll look at on its own.
• Tests pass in test.quantile.js, test.factory.js and test.native.js, and in the lognormal/quantile, normal/ctor, lognormal/ctor and stats/fligner-test dependents.

Checklist

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

• Read, understood, and followed the contributing guidelines.

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

[stdlib-bot]

Coverage Report

Package	Statements	Branches	Functions	Lines
stats/base/dists/normal/quantile	$\\color{green}316/316$ $\\color{green}+100.00\\%$	$\\color{green}27/27$ $\\color{green}+100.00\\%$	$\\color{green}4/4$ $\\color{green}+100.00\\%$	$\\color{green}316/316$ $\\color{green}+100.00\\%$

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