fix: evaluate the upper tail without cancellation in stats/base/dists/laplace/logcdf

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

@@ -22,7 +22,7 @@
 
 var constantFunction = require( '@stdlib/utils/constant-function' );
 var isnan = require( '@stdlib/math/base/assert/is-nan' );
-var expm1 = require( '@stdlib/math/base/special/expm1' );
+var exp = require( '@stdlib/math/base/special/exp' );
 var log1p = require( '@stdlib/math/base/special/log1p' );
 var LNHALF = require( '@stdlib/constants/float64/ln-half' );
 
@@ -71,7 +71,8 @@ function factory( mu, b ) {
 		if ( x < mu ) {
 			return LNHALF + z;
 		}
-		return LNHALF + log1p( -expm1( -z ) );
+		// `ln(1 - exp(-z)/2)`: adding `ln(1/2)` to `ln(2 - exp(-z))` cancels for large `z` and rounds to `0`, so take `log1p` of the small term directly...
+		return log1p( -0.5 * exp( -z ) );
 	}
 }
 

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

@@ -21,7 +21,7 @@
 // MODULES //
 
 var isnan = require( '@stdlib/math/base/assert/is-nan' );
-var expm1 = require( '@stdlib/math/base/special/expm1' );
+var exp = require( '@stdlib/math/base/special/exp' );
 var log1p = require( '@stdlib/math/base/special/log1p' );
 var LNHALF = require( '@stdlib/constants/float64/ln-half' );
 
@@ -75,7 +75,8 @@ function logcdf( x, mu, b ) {
 	if ( x < mu ) {
 		return LNHALF + z;
 	}
-	return LNHALF + log1p( -expm1( -z ) );
+	// `ln(1 - exp(-z)/2)`: adding `ln(1/2)` to `ln(2 - exp(-z))` cancels for large `z` and rounds to `0`, so take `log1p` of the small term directly...
+	return log1p( -0.5 * exp( -z ) );
 }
 
 

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

@@ -42,7 +42,7 @@
         "@stdlib/math/base/assert/is-nan",
         "@stdlib/math/base/special/log1p",
         "@stdlib/constants/float64/ln-half",
-        "@stdlib/math/base/special/expm1"
+        "@stdlib/math/base/special/exp"
       ]
     },
     {
@@ -61,7 +61,7 @@
         "@stdlib/math/base/special/log1p",
         "@stdlib/constants/float64/ln-half",
         "@stdlib/constants/float64/eps",
-        "@stdlib/math/base/special/expm1"
+        "@stdlib/math/base/special/exp"
       ]
     },
     {
@@ -79,7 +79,7 @@
         "@stdlib/math/base/assert/is-nan",
         "@stdlib/math/base/special/log1p",
         "@stdlib/constants/float64/ln-half",
-        "@stdlib/math/base/special/expm1"
+        "@stdlib/math/base/special/exp"
       ]
     }
   ]

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

@@ -18,7 +18,7 @@
 
 #include "stdlib/stats/base/dists/laplace/logcdf.h"
 #include "stdlib/math/base/assert/is_nan.h"
-#include "stdlib/math/base/special/expm1.h"
+#include "stdlib/math/base/special/exp.h"
 #include "stdlib/math/base/special/log1p.h"
 #include "stdlib/constants/float64/ln_half.h"
 
@@ -48,5 +48,6 @@ double stdlib_base_dists_laplace_logcdf( const double x, const double mu, const
 	if ( x < mu ) {
 		return STDLIB_CONSTANT_FLOAT64_LN_HALF + z;
 	}
-	return STDLIB_CONSTANT_FLOAT64_LN_HALF + stdlib_base_log1p( -stdlib_base_expm1( -z ) );
+	// `ln(1 - exp(-z)/2)`: adding `ln(1/2)` to `ln(2 - exp(-z))` cancels for large `z` and rounds to `0`, so take `log1p` of the small term directly...
+	return stdlib_base_log1p( -0.5 * stdlib_base_exp( -z ) );
 }

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

@@ -16,7 +16,6 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-import Distributions: logcdf, Laplace
 import JSON
 
 """
@@ -41,9 +40,10 @@ julia> gen( x, mu, b, "data.json" );
 ```
 """
 function gen( x, mu, b, name )
-	z = Array{Float64}( undef, length(x) );
-	for i in eachindex(x)
-		z[ i ] = logcdf( Laplace( mu[i], b[i] ), x[i] );
+	# The log CDF is evaluated in extended precision, as `ln(1/2) + ln(2 - exp(-z))` cancels in double precision above the mean (Distributions.jl's `logcdf` evaluates it in `Float64` and returns `0` in the upper tail)...
+	z = setprecision( BigFloat, 2048 ) do
+		t = ( BigFloat.( x ) .- BigFloat.( mu ) ) ./ BigFloat.( b );
+		Float64.( [ ti < 0 ? ti + log( BigFloat( 0.5 ) ) : log1p( -0.5 * exp( -ti ) ) for ti in t ] );
 	end
 
 	# Store data to be written to file as a collection:

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

@@ -153,7 +153,7 @@ tape( 'the created function evaluates the logcdf for `x` given positive `mu`', f
 	for ( i = 0; i < x.length; i++ ) {
 		logcdf = factory( mu[i], b[i] );
 		y = logcdf( x[i] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 0 ), true, 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 1 ), true, 'returns expected value' );
 	}
 	t.end();
 });
@@ -174,7 +174,7 @@ tape( 'the created function evaluates the logcdf for `x` given negative `mu`', f
 	for ( i = 0; i < x.length; i++ ) {
 		logcdf = factory( mu[i], b[i] );
 		y = logcdf( x[i] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 0 ), true, 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 1 ), true, 'returns expected value' );
 	}
 	t.end();
 });
@@ -195,7 +195,9 @@ tape( 'the created function evaluates the logcdf for `x` given large variance (
 	for ( i = 0; i < x.length; i++ ) {
 		logcdf = factory( mu[i], b[i] );
 		y = logcdf( x[i] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 0 ), true, 'returns expected value' );
+
+		// Far above the mean the result is about `-exp(-z)/2`, so the rounding of `z = (x-mu)/b` is amplified by `z`; with `z` up to ~120 here that sets the tolerance...
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 160 ), true, 'returns expected value' );
 	}
 	t.end();
 });

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

@@ -103,7 +103,7 @@ tape( 'the function evaluates the logcdf for `x` given positive `mu`', function
 	b = positiveMean.b;
 	for ( i = 0; i < x.length; i++ ) {
 		y = logcdf( x[i], mu[i], b[i] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 0 ), true, 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 1 ), true, 'returns expected value' );
 	}
 	t.end();
 });
@@ -122,7 +122,7 @@ tape( 'the function evaluates the logcdf for `x` given negative `mu`', function
 	b = negativeMean.b;
 	for ( i = 0; i < x.length; i++ ) {
 		y = logcdf( x[i], mu[i], b[i] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 0 ), true, 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 1 ), true, 'returns expected value' );
 	}
 	t.end();
 });
@@ -141,7 +141,9 @@ tape( 'the function evaluates the logcdf for `x` given large variance ( = large
 	b = largeVariance.b;
 	for ( i = 0; i < x.length; i++ ) {
 		y = logcdf( x[i], mu[i], b[i] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 0 ), true, 'returns expected value' );
+
+		// Far above the mean the result is about `-exp(-z)/2`, so the rounding of `z = (x-mu)/b` is amplified by `z`; with `z` up to ~120 here that sets the tolerance...
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 160 ), true, 'returns expected value' );
 	}
 	t.end();
 });

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

@@ -112,7 +112,7 @@ tape( 'the function evaluates the logcdf for `x` given positive `mu`', opts, fun
 	b = positiveMean.b;
 	for ( i = 0; i < x.length; i++ ) {
 		y = logcdf( x[i], mu[i], b[i] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 0 ), true, 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 1 ), true, 'returns expected value' );
 	}
 	t.end();
 });
@@ -131,7 +131,7 @@ tape( 'the function evaluates the logcdf for `x` given negative `mu`', opts, fun
 	b = negativeMean.b;
 	for ( i = 0; i < x.length; i++ ) {
 		y = logcdf( x[i], mu[i], b[i] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 0 ), true, 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 1 ), true, 'returns expected value' );
 	}
 	t.end();
 });
@@ -150,7 +150,9 @@ tape( 'the function evaluates the logcdf for `x` given large variance ( = large
 	b = largeVariance.b;
 	for ( i = 0; i < x.length; i++ ) {
 		y = logcdf( x[i], mu[i], b[i] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 0 ), true, 'returns expected value' );
+
+		// Far above the mean the result is about `-exp(-z)/2`, so the rounding of `z = (x-mu)/b` is amplified by `z`; with `z` up to ~120 here that sets the tolerance...
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 160 ), true, 'returns expected value' );
 	}
 	t.end();
 });