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

lib/node_modules/@stdlib/stats/base/dists/normal/quantile/lib/factory.js

@@ -22,7 +22,7 @@
 
 var constantFunction = require( '@stdlib/utils/constant-function' );
 var degenerate = require( '@stdlib/stats/base/dists/degenerate/quantile' ).factory;
-var erfinv = require( '@stdlib/math/base/special/erfinv' );
+var erfcinv = require( '@stdlib/math/base/special/erfcinv' );
 var isnan = require( '@stdlib/math/base/assert/is-nan' );
 var sqrt = require( '@stdlib/math/base/special/sqrt' );
 
@@ -72,7 +72,7 @@ function factory( mu, sigma ) {
 		if ( isnan( p ) || p < 0.0 || p > 1.0 ) {
 			return NaN;
 		}
-		return A + ( B * erfinv( (2.0*p) - 1.0 ) );
+		return A - ( B * erfcinv( 2.0*p ) );
 	}
 }
 

lib/node_modules/@stdlib/stats/base/dists/normal/quantile/lib/main.js

@@ -20,7 +20,7 @@
 
 // MODULES //
 
-var erfinv = require( '@stdlib/math/base/special/erfinv' );
+var erfcinv = require( '@stdlib/math/base/special/erfcinv' );
 var isnan = require( '@stdlib/math/base/assert/is-nan' );
 var sqrt = require( '@stdlib/math/base/special/sqrt' );
 
@@ -30,6 +30,16 @@ var sqrt = require( '@stdlib/math/base/special/sqrt' );
 /**
 * Evaluates the quantile function for a normal distribution with mean `mu` and standard deviation `sigma` at a probability `p`.
 *
+* ## Method
+*
+* The quantile function is
+*
+* ```tex
+* Q(p) = \mu + \sigma \sqrt{2} \operatorname{erf}^{-1}(2p-1) = \mu - \sigma \sqrt{2} \operatorname{erfc}^{-1}(2p)
+* ```
+*
+* The function evaluates the second form. In the first, \\(2p-1\\) loses the low-order digits of small \\(p\\) and rounds to \\(-1\\) for \\(p < 2^{-54}\\), which returns \\(-\infty\\).
+*
 * @param {Probability} p - input value
 * @param {number} mu - mean
 * @param {NonNegativeNumber} sigma - standard deviation
@@ -87,7 +97,7 @@ function quantile( p, mu, sigma ) {
 	}
 	A = mu;
 	B = sigma * sqrt( 2.0 );
-	return A + (B * erfinv( (2.0*p) - 1.0 ));
+	return A - ( B * erfcinv( 2.0*p ) );
 }
 
 

lib/node_modules/@stdlib/stats/base/dists/normal/quantile/manifest.json

@@ -41,7 +41,7 @@
         "@stdlib/math/base/napi/ternary",
         "@stdlib/math/base/assert/is-nan",
         "@stdlib/math/base/special/sqrt",
-        "@stdlib/math/base/special/erfinv",
+        "@stdlib/math/base/special/erfcinv",
         "@stdlib/constants/float64/sqrt-two"
       ]
     },
@@ -59,7 +59,7 @@
       "dependencies": [
         "@stdlib/math/base/assert/is-nan",
         "@stdlib/math/base/special/sqrt",
-        "@stdlib/math/base/special/erfinv",
+        "@stdlib/math/base/special/erfcinv",
         "@stdlib/constants/float64/eps",
         "@stdlib/constants/float64/sqrt-two"
       ]
@@ -78,7 +78,7 @@
       "dependencies": [
         "@stdlib/math/base/assert/is-nan",
         "@stdlib/math/base/special/sqrt",
-        "@stdlib/math/base/special/erfinv",
+        "@stdlib/math/base/special/erfcinv",
         "@stdlib/constants/float64/eps",
         "@stdlib/constants/float64/sqrt-two"
       ]

lib/node_modules/@stdlib/stats/base/dists/normal/quantile/src/main.c

@@ -19,12 +19,22 @@
 #include "stdlib/stats/base/dists/normal/quantile.h"
 #include "stdlib/math/base/assert/is_nan.h"
 #include "stdlib/math/base/special/sqrt.h"
-#include "stdlib/math/base/special/erfinv.h"
+#include "stdlib/math/base/special/erfcinv.h"
 #include "stdlib/constants/float64/sqrt_two.h"
 
 /**
 * Evaluates the quantile function for a normal distribution with mean `mu` and standard deviation `sigma` at a probability `p`.
 *
+* ## Method
+*
+* The quantile function is
+*
+* ```tex
+* Q(p) = \mu + \sigma \sqrt{2} \operatorname{erf}^{-1}(2p-1) = \mu - \sigma \sqrt{2} \operatorname{erfc}^{-1}(2p)
+* ```
+*
+* The function evaluates the second form. In the first, `2p-1` loses the low-order digits of small `p` and rounds to `-1` for `p < 2^{-54}`, which returns `-\infty`.
+*
 * @param p      probability
 * @param mu     mean
 * @param sigma  standard deviation
@@ -48,5 +58,5 @@ double stdlib_base_dists_normal_quantile( const double p, const double mu, const
 	if ( sigma == 0.0 ) {
 		return mu;
 	}
-	return mu + ( ( sigma * STDLIB_CONSTANT_FLOAT64_SQRT2 ) * stdlib_base_erfinv( ( 2.0 * p ) - 1.0 ) );
+	return mu - ( ( sigma * STDLIB_CONSTANT_FLOAT64_SQRT2 ) * stdlib_base_erfcinv( 2.0 * p ) );
 }

lib/node_modules/@stdlib/stats/base/dists/normal/quantile/test/fixtures/julia/REQUIRE

@@ -0,0 +1,3 @@
+Distributions 0.25
+julia 1.5
+JSON 0.21

lib/node_modules/@stdlib/stats/base/dists/normal/quantile/test/fixtures/julia/runner.jl

@@ -0,0 +1,87 @@
+#!/usr/bin/env julia
+#
+# @license Apache-2.0
+#
+# Copyright (c) 2026 The Stdlib Authors.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#    http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+import Distributions: quantile, Normal
+import JSON
+
+"""
+	gen( p, mu, sigma, name )
+
+Generate fixture data and write to file.
+
+# Arguments
+
+* `p`: probabilities
+* `mu`: mean
+* `sigma`: standard deviation
+* `name::AbstractString`: output filename
+
+# Examples
+
+``` julia
+julia> p = rand( 1000 );
+julia> mu = rand( 1000 );
+julia> sigma = rand( 1000 ) .+ 0.1;
+julia> gen( p, mu, sigma, "data.json" );
+```
+"""
+function gen( p, mu, sigma, name )
+	z = Array{Float64}( undef, length( p ) );
+
+	# Evaluate the quantile function in extended precision, and round the result to double precision:
+	setprecision( BigFloat, 2048 ) do
+		for i in eachindex( p )
+			z[ i ] = Float64( quantile( Normal( BigFloat( mu[ i ] ), BigFloat( sigma[ i ] ) ), BigFloat( p[ i ] ) ) );
+		end
+	end
+
+	# Store data to be written to file as a collection:
+	data = Dict([
+		("p", p),
+		("mu", mu),
+		("sigma", sigma),
+		("expected", z)
+	]);
+
+	# Based on the script directory, create an output filepath:
+	filepath = joinpath( dir, name );
+
+	# Write the data to the output filepath as JSON:
+	outfile = open( filepath, "w" );
+	write( outfile, JSON.json(data) );
+	write( outfile, "\n" );
+	close( outfile );
+end
+
+# Get the filename:
+file = @__FILE__;
+
+# Extract the directory in which this file resides:
+dir = dirname( file );
+
+# Probabilities in the lower tail, down to `1e-250`:
+p = exp10.( rand( 500 ) .* -247.0 .- 3.0 );
+mu = ( rand( 500 ) .* 20.0 ) .- 10.0;
+sigma = ( rand( 500 ) .* 10.0 ) .+ 0.1;
+gen( p, mu, sigma, "lower_tail.json" );
+
+# Probabilities in the upper tail, up to the largest double less than `1`:
+p = 1.0 .- exp10.( rand( 500 ) .* -13.0 .- 3.0 );
+mu = ( rand( 500 ) .* 20.0 ) .- 10.0;
+sigma = ( rand( 500 ) .* 10.0 ) .+ 0.1;
+gen( p, mu, sigma, "upper_tail.json" );

lib/node_modules/@stdlib/stats/base/dists/normal/quantile/test/test.factory.js

@@ -31,6 +31,8 @@ var factory = require( './../lib/factory.js' );
 // FIXTURES //
 
 var data = require( './fixtures/r/data.json' );
+var lowerTail = require( './fixtures/julia/lower_tail.json' );
+var upperTail = require( './fixtures/julia/upper_tail.json' );
 var i;
 for ( i = 0; i < data.expected.length; i++ ) {
 	if ( data.expected[ i ] === 'Inf' ) {
@@ -169,3 +171,45 @@ tape( 'the created function evaluates the quantile for `x` given parameters `mu`
 	}
 	t.end();
 });
+
+tape( 'the created function evaluates the quantile for small probabilities `p` (lower tail)', function test( t ) {
+	var expected;
+	var quantile;
+	var sigma;
+	var mu;
+	var p;
+	var y;
+	var i;
+
+	expected = lowerTail.expected;
+	p = lowerTail.p;
+	mu = lowerTail.mu;
+	sigma = lowerTail.sigma;
+	for ( i = 0; i < p.length; i++ ) {
+		quantile = factory( mu[i], sigma[i] );
+		y = quantile( p[i] );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 16 ), true, 'p: '+p[i]+'. y: '+y+'. E: '+expected[i] );
+	}
+	t.end();
+});
+
+tape( 'the created function evaluates the quantile for probabilities `p` close to one (upper tail)', function test( t ) {
+	var expected;
+	var quantile;
+	var sigma;
+	var mu;
+	var p;
+	var y;
+	var i;
+
+	expected = upperTail.expected;
+	p = upperTail.p;
+	mu = upperTail.mu;
+	sigma = upperTail.sigma;
+	for ( i = 0; i < p.length; i++ ) {
+		quantile = factory( mu[i], sigma[i] );
+		y = quantile( p[i] );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 96 ), true, 'p: '+p[i]+'. y: '+y+'. E: '+expected[i] );
+	}
+	t.end();
+});

lib/node_modules/@stdlib/stats/base/dists/normal/quantile/test/test.native.js

@@ -40,6 +40,8 @@ var opts = {
 // FIXTURES //
 
 var data = require( './fixtures/r/data.json' );
+var lowerTail = require( './fixtures/julia/lower_tail.json' );
+var upperTail = require( './fixtures/julia/upper_tail.json' );
 var i;
 for ( i = 0; i < data.expected.length; i++ ) {
 	if ( data.expected[ i ] === 'Inf' ) {
@@ -137,3 +139,41 @@ tape( 'the function evaluates the quantile for `x` given parameters `mu` and `si
 	}
 	t.end();
 });
+
+tape( 'the function evaluates the quantile for small probabilities `p` (lower tail)', opts, function test( t ) {
+	var expected;
+	var sigma;
+	var mu;
+	var p;
+	var y;
+	var i;
+
+	expected = lowerTail.expected;
+	p = lowerTail.p;
+	mu = lowerTail.mu;
+	sigma = lowerTail.sigma;
+	for ( i = 0; i < p.length; i++ ) {
+		y = quantile( p[i], mu[i], sigma[i] );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 16 ), true, 'p: '+p[i]+'. y: '+y+'. E: '+expected[i] );
+	}
+	t.end();
+});
+
+tape( 'the function evaluates the quantile for probabilities `p` close to one (upper tail)', opts, function test( t ) {
+	var expected;
+	var sigma;
+	var mu;
+	var p;
+	var y;
+	var i;
+
+	expected = upperTail.expected;
+	p = upperTail.p;
+	mu = upperTail.mu;
+	sigma = upperTail.sigma;
+	for ( i = 0; i < p.length; i++ ) {
+		y = quantile( p[i], mu[i], sigma[i] );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 96 ), true, 'p: '+p[i]+'. y: '+y+'. E: '+expected[i] );
+	}
+	t.end();
+});

lib/node_modules/@stdlib/stats/base/dists/normal/quantile/test/test.quantile.js

@@ -31,6 +31,8 @@ var quantile = require( './../lib' );
 // FIXTURES //
 
 var data = require( './fixtures/r/data.json' );
+var lowerTail = require( './fixtures/julia/lower_tail.json' );
+var upperTail = require( './fixtures/julia/upper_tail.json' );
 var i;
 for ( i = 0; i < data.expected.length; i++ ) {
 	if ( data.expected[ i ] === 'Inf' ) {
@@ -128,3 +130,41 @@ tape( 'the function evaluates the quantile for `x` given parameters `mu` and `si
 	}
 	t.end();
 });
+
+tape( 'the function evaluates the quantile for small probabilities `p` (lower tail)', function test( t ) {
+	var expected;
+	var sigma;
+	var mu;
+	var p;
+	var y;
+	var i;
+
+	expected = lowerTail.expected;
+	p = lowerTail.p;
+	mu = lowerTail.mu;
+	sigma = lowerTail.sigma;
+	for ( i = 0; i < p.length; i++ ) {
+		y = quantile( p[i], mu[i], sigma[i] );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 16 ), true, 'p: '+p[i]+'. y: '+y+'. E: '+expected[i] );
+	}
+	t.end();
+});
+
+tape( 'the function evaluates the quantile for probabilities `p` close to one (upper tail)', function test( t ) {
+	var expected;
+	var sigma;
+	var mu;
+	var p;
+	var y;
+	var i;
+
+	expected = upperTail.expected;
+	p = upperTail.p;
+	mu = upperTail.mu;
+	sigma = upperTail.sigma;
+	for ( i = 0; i < p.length; i++ ) {
+		y = quantile( p[i], mu[i], sigma[i] );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 96 ), true, 'p: '+p[i]+'. y: '+y+'. E: '+expected[i] );
+	}
+	t.end();
+});