fix: evaluate the lower tail from p in stats/base/dists/laplace/quantile

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

@@ -21,9 +21,7 @@
 // MODULES //
 
 var constantFunction = require( '@stdlib/utils/constant-function' );
-var signum = require( '@stdlib/math/base/special/signum' );
 var isnan = require( '@stdlib/math/base/assert/is-nan' );
-var abs = require( '@stdlib/math/base/special/abs' );
 var ln = require( '@stdlib/math/base/special/ln' );
 
 
@@ -69,7 +67,11 @@ function factory( mu, b ) {
 		if ( isnan( p ) || p < 0.0 || p > 1.0 ) {
 			return NaN;
 		}
-		return mu - ( b * signum( p-0.5 ) * ln( 1.0 - (2.0 * abs( p-0.5 )) ) );
+		// Work with the tail probability (`p` or `1-p`, which is exact for `p >= 0.5`) rather than `1 - 2|p-0.5|`, which rounds away the low-order bits of a small `p`:
+		if ( p < 0.5 ) {
+			return mu + ( b * ln( 2.0*p ) );
+		}
+		return mu - ( b * ln( 2.0*(1.0-p) ) );
 	}
 }
 

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

@@ -20,9 +20,7 @@
 
 // MODULES //
 
-var signum = require( '@stdlib/math/base/special/signum' );
 var isnan = require( '@stdlib/math/base/assert/is-nan' );
-var abs = require( '@stdlib/math/base/special/abs' );
 var ln = require( '@stdlib/math/base/special/ln' );
 
 
@@ -80,7 +78,11 @@ function quantile( p, mu, b ) {
 	) {
 		return NaN;
 	}
-	return mu - ( b * signum( p-0.5 ) * ln( 1.0 - ( 2.0 * abs( p-0.5 ) ) ) );
+	// Work with the tail probability (`p` or `1-p`, which is exact for `p >= 0.5`) rather than `1 - 2|p-0.5|`, which rounds away the low-order bits of a small `p`:
+	if ( p < 0.5 ) {
+		return mu + ( b * ln( 2.0*p ) );
+	}
+	return mu - ( b * ln( 2.0*(1.0-p) ) );
 }
 
 

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

@@ -39,9 +39,7 @@
       "libpath": [],
       "dependencies": [
         "@stdlib/math/base/napi/ternary",
-        "@stdlib/math/base/special/signum",
         "@stdlib/math/base/assert/is-nan",
-        "@stdlib/math/base/special/abs",
         "@stdlib/math/base/special/ln"
       ]
     },
@@ -57,9 +55,7 @@
       "libraries": [],
       "libpath": [],
       "dependencies": [
-        "@stdlib/math/base/special/signum",
         "@stdlib/math/base/assert/is-nan",
-        "@stdlib/math/base/special/abs",
         "@stdlib/constants/float64/eps",
         "@stdlib/math/base/special/ln"
       ]
@@ -76,9 +72,7 @@
       "libraries": [],
       "libpath": [],
       "dependencies": [
-        "@stdlib/math/base/special/signum",
         "@stdlib/math/base/assert/is-nan",
-        "@stdlib/math/base/special/abs",
         "@stdlib/math/base/special/ln"
       ]
     }

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

@@ -17,9 +17,7 @@
 */
 
 #include "stdlib/stats/base/dists/laplace/quantile.h"
-#include "stdlib/math/base/special/signum.h"
 #include "stdlib/math/base/assert/is_nan.h"
-#include "stdlib/math/base/special/abs.h"
 #include "stdlib/math/base/special/ln.h"
 
 /**
@@ -45,5 +43,9 @@ double stdlib_base_dists_laplace_quantile( const double p, const double mu, cons
 	) {
 		return 0.0/0.0; // NaN
 	}
-	return mu - ( b * stdlib_base_signum( p-0.5 ) * stdlib_base_ln( 1.0 - ( 2.0 * stdlib_base_abs( p-0.5 ) ) ) );
+	// Work with the tail probability (`p` or `1-p`, which is exact for `p >= 0.5`) rather than `1 - 2|p-0.5|`, which rounds away the low-order bits of a small `p`:
+	if ( p < 0.5 ) {
+		return mu + ( b * stdlib_base_ln( 2.0*p ) );
+	}
+	return mu - ( b * stdlib_base_ln( 2.0*(1.0-p) ) );
 }

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

@@ -64,6 +64,42 @@ function gen( p, mu, b, name )
 	close( outfile );
 end
 
+"""
+	gen_tails( p, mu, b, name )
+
+Generate fixture data for the tails in extended precision and write to file. In double precision, `1 - 2|p-0.5|` rounds away the low-order bits of a small `p` (all of them for `p < 2^-54`), so the reference value is computed from the textbook formula `mu - b*sign(p-0.5)*ln(1 - 2|p-0.5|)` at 2048 bits.
+
+# Arguments
+
+* `p`: input value
+* `mu`: location parameter
+* `b`: scale parameter
+* `name::AbstractString`: output filename
+"""
+function gen_tails( p, mu, b, name )
+	z = setprecision( BigFloat, 2048 ) do
+		q = BigFloat.( p );
+		Float64.( BigFloat.( mu ) .- ( BigFloat.( b ) .* sign.( q .- 0.5 ) .* log.( 1 .- ( 2 .* abs.( q .- 0.5 ) ) ) ) )
+	end
+
+	# Store data to be written to file as a collection:
+	data = Dict([
+		("p", p),
+		("mu", mu),
+		("b", b),
+		("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__;
 
@@ -87,3 +123,9 @@ p = rand( 1000 );
 mu = rand( 1000 );
 b = rand( 1000 ) .* 20.0;
 gen( p, mu, b, "large_variance.json" );
+
+# Tails (`p` and `1-p` log-spaced; `mu = 0`, so that the result does not also measure the rounding of `mu + ...`):
+p = vcat( exp10.( range( -300.0, stop = -1.0, length = 500 ) ), 1.0 .- exp10.( range( -16.0, stop = -1.0, length = 500 ) ) );
+mu = zeros( 1000 );
+b = ( rand( 1000 ) .* 5.0 ) .+ 0.1;
+gen_tails( p, mu, b, "tails.json" );

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

@@ -33,6 +33,7 @@ var factory = require( './../lib/factory.js' );
 var positiveMean = require( './fixtures/julia/positive_mean.json' );
 var negativeMean = require( './fixtures/julia/negative_mean.json' );
 var largeVariance = require( './fixtures/julia/large_variance.json' );
+var tails = require( './fixtures/julia/tails.json' );
 
 
 // TESTS //
@@ -197,3 +198,26 @@ tape( 'the created function evaluates the quantile function at `p` given large v
 	}
 	t.end();
 });
+
+tape( 'the created function evaluates the quantile function at `p` in the tails (`p` or `1-p` close to `0`)', function test( t ) {
+	var expected;
+	var quantile;
+	var mu;
+	var b;
+	var i;
+	var p;
+	var y;
+
+	expected = tails.expected;
+	p = tails.p;
+	mu = tails.mu;
+	b = tails.b;
+	for ( i = 0; i < p.length; i++ ) {
+		quantile = factory( mu[i], b[i] );
+		y = quantile( p[i] );
+		if ( expected[i] !== null ) {
+			t.strictEqual( isAlmostSameValue( y, expected[i], 1 ), true, 'returns expected value' );
+		}
+	}
+	t.end();
+});

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

@@ -42,6 +42,7 @@ var opts = {
 var positiveMean = require( './fixtures/julia/positive_mean.json' );
 var negativeMean = require( './fixtures/julia/negative_mean.json' );
 var largeVariance = require( './fixtures/julia/large_variance.json' );
+var tails = require( './fixtures/julia/tails.json' );
 
 
 // TESTS //
@@ -162,3 +163,24 @@ tape( 'the function evaluates the quantile function at `p` given large variance
 	}
 	t.end();
 });
+
+tape( 'the function evaluates the quantile function at `p` in the tails (`p` or `1-p` close to `0`)', opts, function test( t ) {
+	var expected;
+	var mu;
+	var b;
+	var i;
+	var p;
+	var y;
+
+	expected = tails.expected;
+	p = tails.p;
+	mu = tails.mu;
+	b = tails.b;
+	for ( i = 0; i < p.length; i++ ) {
+		y = quantile( p[i], mu[i], b[i] );
+		if ( expected[i] !== null ) {
+			t.strictEqual( isAlmostSameValue( y, expected[i], 1 ), true, 'returns expected value' );
+		}
+	}
+	t.end();
+});

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

@@ -33,6 +33,7 @@ var quantile = require( './../lib' );
 var positiveMean = require( './fixtures/julia/positive_mean.json' );
 var negativeMean = require( './fixtures/julia/negative_mean.json' );
 var largeVariance = require( './fixtures/julia/large_variance.json' );
+var tails = require( './fixtures/julia/tails.json' );
 
 
 // TESTS //
@@ -153,3 +154,24 @@ tape( 'the function evaluates the quantile function at `p` given large variance
 	}
 	t.end();
 });
+
+tape( 'the function evaluates the quantile function at `p` in the tails (`p` or `1-p` close to `0`)', function test( t ) {
+	var expected;
+	var mu;
+	var b;
+	var i;
+	var p;
+	var y;
+
+	expected = tails.expected;
+	p = tails.p;
+	mu = tails.mu;
+	b = tails.b;
+	for ( i = 0; i < p.length; i++ ) {
+		y = quantile( p[i], mu[i], b[i] );
+		if ( expected[i] !== null ) {
+			t.strictEqual( isAlmostSameValue( y, expected[i], 1 ), true, 'returns expected value' );
+		}
+	}
+	t.end();
+});