fix: avoid underflow for small x in stats/base/dists/weibull/logcdf

lib/node_modules/@stdlib/stats/base/dists/weibull/logcdf/README.md

@@ -68,7 +68,7 @@ y = logcdf( -Infinity, 4.0, 2.0 );
 // returns -Infinity
 
 y = logcdf( +Infinity, 4.0, 2.0 );
-// returns 0.0
+// returns -0.0
 ```
 
 If provided `NaN` as any argument, the function returns `NaN`.

lib/node_modules/@stdlib/stats/base/dists/weibull/logcdf/docs/repl.txt

@@ -31,7 +31,7 @@
     > y = {{alias}}( -1.0, 2.0, 2.0 )
     -Infinity
     > y = {{alias}}( {{alias:@stdlib/constants/float64/pinf}}, 4.0, 2.0 )
-    0.0
+    -0.0
     > y = {{alias}}( {{alias:@stdlib/constants/float64/ninf}}, 4.0, 2.0 )
     -Infinity
     > y = {{alias}}( NaN, 0.0, 1.0 )

lib/node_modules/@stdlib/stats/base/dists/weibull/logcdf/docs/types/index.d.ts

@@ -52,7 +52,7 @@ interface LogCDF {
 	*
 	* @example
 	* var y = logcdf( +Infinity, 4.0, 2.0 );
-	* // returns 0.0
+	* // returns -0.0
 	*
 	* @example
 	* var y = logcdf( -Infinity, 4.0, 2.0 );

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

@@ -28,6 +28,7 @@ var exp = require( '@stdlib/math/base/special/exp' );
 var pow = require( '@stdlib/math/base/special/pow' );
 var ln = require( '@stdlib/math/base/special/ln' );
 var LNHALF = require( '@stdlib/constants/float64/ln-half' );
+var FLOAT64_SMALLEST_NORMAL = require( '@stdlib/constants/float64/smallest-normal' );
 var NINF = require( '@stdlib/constants/float64/ninf' );
 
 
@@ -71,13 +72,25 @@ function factory( k, lambda ) {
 	* // returns <number>
 	*/
 	function logcdf( x ) {
+		var lt;
+		var r;
 		var p;
 		if ( isnan( x ) ) {
 			return NaN;
 		}
 		if ( x < 0.0 ) {
 			return NINF;
 		}
+		r = x / lambda;
+		if ( r < FLOAT64_SMALLEST_NORMAL ) {
+			lt = k * ( ln( x ) - ln( lambda ) );
+		} else {
+			lt = k * ln( r );
+		}
+		if ( lt < -36.0 ) {
+			// For `t = (x/lambda)^k < e^-36`, `ln(1 - e^-t) = ln(t) + ln(1 - t/2 + ...)` and the correction is below half an ULP of `ln(t)`. Using `ln(t)` directly avoids `t` underflowing to `0`:
+			return lt;
+		}
 		p = -pow( x / lambda, k );
 		return ( p < LNHALF ) ? log1p( -exp( p ) ) : ln( -expm1( p ) );
 	}

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

@@ -27,6 +27,7 @@ var exp = require( '@stdlib/math/base/special/exp' );
 var pow = require( '@stdlib/math/base/special/pow' );
 var ln = require( '@stdlib/math/base/special/ln' );
 var LNHALF = require( '@stdlib/constants/float64/ln-half' );
+var FLOAT64_SMALLEST_NORMAL = require( '@stdlib/constants/float64/smallest-normal' );
 var NINF = require( '@stdlib/constants/float64/ninf' );
 
 
@@ -50,7 +51,7 @@ var NINF = require( '@stdlib/constants/float64/ninf' );
 *
 * @example
 * var y = logcdf( +Infinity, 4.0, 2.0 );
-* // returns 0.0
+* // returns -0.0
 *
 * @example
 * var y = logcdf( -Infinity, 4.0, 2.0 );
@@ -73,6 +74,8 @@ var NINF = require( '@stdlib/constants/float64/ninf' );
 * // returns NaN
 */
 function logcdf( x, k, lambda ) {
+	var lt;
+	var r;
 	var p;
 	if (
 		isnan( k ) ||
@@ -85,6 +88,16 @@ function logcdf( x, k, lambda ) {
 	if ( x < 0.0 ) {
 		return NINF;
 	}
+	r = x / lambda;
+	if ( r < FLOAT64_SMALLEST_NORMAL ) {
+		lt = k * ( ln( x ) - ln( lambda ) );
+	} else {
+		lt = k * ln( r );
+	}
+	if ( lt < -36.0 ) {
+		// For `t = (x/lambda)^k < e^-36`, `ln(1 - e^-t) = ln(t) + ln(1 - t/2 + ...)` and the correction is below half an ULP of `ln(t)`. Using `ln(t)` directly avoids `t` underflowing to `0`:
+		return lt;
+	}
 	p = -pow( x / lambda, k );
 	return ( p < LNHALF ) ? log1p( -exp( p ) ) : ln( -expm1( p ) );
 }

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

@@ -44,7 +44,7 @@ var addon = require( './../src/addon.node' );
 *
 * @example
 * var y = logcdf( +Infinity, 4.0, 2.0 );
-* // returns 0.0
+* // returns -0.0
 *
 * @example
 * var y = logcdf( -Infinity, 4.0, 2.0 );

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

@@ -47,6 +47,7 @@
         "@stdlib/math/base/special/ln",
         "@stdlib/math/base/special/pow",
         "@stdlib/constants/float64/ln-half",
+        "@stdlib/constants/float64/smallest-normal",
         "@stdlib/constants/float64/ninf"
       ]
     },
@@ -71,6 +72,7 @@
         "@stdlib/math/base/special/ln",
         "@stdlib/math/base/special/pow",
         "@stdlib/constants/float64/ln-half",
+        "@stdlib/constants/float64/smallest-normal",
         "@stdlib/constants/float64/ninf"
       ]
     },
@@ -94,6 +96,7 @@
         "@stdlib/math/base/special/ln",
         "@stdlib/math/base/special/pow",
         "@stdlib/constants/float64/ln-half",
+        "@stdlib/constants/float64/smallest-normal",
         "@stdlib/constants/float64/ninf"
       ]
     }

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

@@ -24,6 +24,7 @@
 #include "stdlib/math/base/special/exp.h"
 #include "stdlib/math/base/special/ln.h"
 #include "stdlib/constants/float64/ln_half.h"
+#include "stdlib/constants/float64/smallest_normal.h"
 #include "stdlib/constants/float64/ninf.h"
 
 /**
@@ -39,6 +40,8 @@
 * // returns ~-0.145
 */
 double stdlib_base_dists_weibull_logcdf( const double x, const double k, const double lambda ) {
+	double lt;
+	double r;
 	double p;
 	if (
 		stdlib_base_is_nan( x ) ||
@@ -51,6 +54,16 @@ double stdlib_base_dists_weibull_logcdf( const double x, const double k, const d
 	if ( x < 0.0 ) {
 		return STDLIB_CONSTANT_FLOAT64_NINF;
 	}
+	r = x / lambda;
+	if ( r < STDLIB_CONSTANT_FLOAT64_SMALLEST_NORMAL ) {
+		lt = k * ( stdlib_base_ln( x ) - stdlib_base_ln( lambda ) );
+	} else {
+		lt = k * stdlib_base_ln( r );
+	}
+	if ( lt < -36.0 ) {
+		// For `t = (x/lambda)^k < e^-36`, `ln(1 - e^-t) = ln(t) + ln(1 - t/2 + ...)` and the correction is below half an ULP of `ln(t)`. Using `ln(t)` directly avoids `t` underflowing to `0`:
+		return lt;
+	}
 	p = -stdlib_base_pow( x / lambda, k );
 	return ( p < STDLIB_CONSTANT_FLOAT64_LN_HALF ) ? stdlib_base_log1p( -stdlib_base_exp( p ) ) : stdlib_base_ln( -stdlib_base_expm1( p ) );
 }

lib/node_modules/@stdlib/stats/base/dists/weibull/logcdf/test/fixtures/julia/REQUIRE

@@ -1,3 +1,2 @@
-Distributions 0.23.8
 julia 1.5
 JSON 0.21

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

@@ -16,14 +16,15 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-import Distributions: logcdf, Weibull
 import JSON
 
 """
 	gen( x, k, lambda, name )
 
 Generate fixture data and write to file.
 
+The logarithm of the CDF is evaluated in extended precision from the textbook formula `ln(1 - exp(-(x/lambda)^k))`, written with `expm1` so that it stays exact when `(x/lambda)^k` is far below the working precision, as `(x/lambda)^k` underflows in double precision for small `x`.
+
 # Arguments
 
 * `x`: input value
@@ -41,9 +42,8 @@ julia> gen( x, k, lambda, "data.json" );
 ```
 """
 function gen( x, k, lambda, name )
-	z = Array{Float64}( undef, length(x) );
-	for i in eachindex(x)
-		z[ i ] = logcdf( Weibull( k[i], lambda[i] ), x[i] );
+	z = setprecision( BigFloat, 2048 ) do
+		Float64.( log.( .-expm1.( .-( ( BigFloat.( x ) ./ BigFloat.( lambda ) ) .^ BigFloat.( k ) ) ) ) )
 	end
 
 	# Store data to be written to file as a collection:
@@ -87,3 +87,9 @@ x = rand( 1000 ) .* 5.0;
 lambda = ( rand( 1000 ) .* 5.0 ) .+ 10.0;
 k = ( rand( 1000 ) .* 5.0 ) .+ 10.0;
 gen( x, k, lambda, "both_large.json" );
+
+# Small `x`:
+x = exp10.( range( -323.0, stop = -2.0, length = 1000 ) );
+lambda = ( rand( 1000 ) .* 5.0 ) .+ 0.5;
+k = ( rand( 1000 ) .* 5.0 ) .+ 0.5;
+gen( x, k, lambda, "small_x.json" );

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

@@ -31,6 +31,7 @@ var factory = require( './../lib/factory.js' );
 // FIXTURES //
 
 var largeScale = require( './fixtures/julia/large_scale.json' );
+var smallX = require( './fixtures/julia/small_x.json' );
 var largeShape = require( './fixtures/julia/large_shape.json' );
 var bothLarge = require( './fixtures/julia/both_large.json' );
 
@@ -192,7 +193,7 @@ tape( 'the created function evaluates the logcdf for `x` given large `lambda` an
 	for ( i = 0; i < x.length; i++ ) {
 		logcdf = factory( k[i], lambda[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();
 });
@@ -213,7 +214,7 @@ tape( 'the created function evaluates the logcdf for `x` given large shape param
 	for ( i = 0; i < x.length; i++ ) {
 		logcdf = factory( k[i], lambda[i] );
 		y = logcdf( x[i] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 0 ), true, 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 3236 ), true, 'returns expected value' );
 	}
 	t.end();
 });
@@ -234,7 +235,28 @@ tape( 'the created function evaluates the logcdf for `x` given large scale param
 	for ( i = 0; i < x.length; i++ ) {
 		logcdf = factory( k[i], lambda[i] );
 		y = logcdf( x[i] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 0 ), true, 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 2 ), true, 'returns expected value' );
+	}
+	t.end();
+});
+
+tape( 'the created function evaluates the logcdf for `x` close to `0`', function test( t ) {
+	var expected;
+	var logcdf;
+	var lambda;
+	var i;
+	var k;
+	var x;
+	var y;
+
+	expected = smallX.expected;
+	x = smallX.x;
+	lambda = smallX.lambda;
+	k = smallX.k;
+	for ( i = 0; i < x.length; i++ ) {
+		logcdf = factory( k[i], lambda[i] );
+		y = logcdf( x[i] );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 1 ), true, 'returns expected value' );
 	}
 	t.end();
 });

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

@@ -31,6 +31,7 @@ var logcdf = require( './../lib' );
 // FIXTURES //
 
 var largeScale = require( './fixtures/julia/large_scale.json' );
+var smallX = require( './fixtures/julia/small_x.json' );
 var largeShape = require( './fixtures/julia/large_shape.json' );
 var bothLarge = require( './fixtures/julia/both_large.json' );
 
@@ -133,7 +134,7 @@ tape( 'the function evaluates the logcdf for `x` given large `lambda` and `k`',
 	k = bothLarge.k;
 	for ( i = 0; i < x.length; i++ ) {
 		y = logcdf( x[i], k[i], lambda[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();
 });
@@ -152,7 +153,7 @@ tape( 'the function evaluates the logcdf for `x` given large shape parameter `la
 	k = largeShape.k;
 	for ( i = 0; i < x.length; i++ ) {
 		y = logcdf( x[i], k[i], lambda[i] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 0 ), true, 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 3236 ), true, 'returns expected value' );
 	}
 	t.end();
 });
@@ -171,7 +172,26 @@ tape( 'the function evaluates the logcdf for `x` given large scale parameter `k`
 	k = largeScale.k;
 	for ( i = 0; i < x.length; i++ ) {
 		y = logcdf( x[i], k[i], lambda[i] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 0 ), true, 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 2 ), true, 'returns expected value' );
+	}
+	t.end();
+});
+
+tape( 'the function evaluates the logcdf for `x` close to `0`', function test( t ) {
+	var expected;
+	var lambda;
+	var i;
+	var k;
+	var x;
+	var y;
+
+	expected = smallX.expected;
+	x = smallX.x;
+	lambda = smallX.lambda;
+	k = smallX.k;
+	for ( i = 0; i < x.length; i++ ) {
+		y = logcdf( x[i], k[i], lambda[i] );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 1 ), true, 'returns expected value' );
 	}
 	t.end();
 });

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

@@ -32,6 +32,7 @@ var NINF = require( '@stdlib/constants/float64/ninf' );
 // FIXTURES //
 
 var largeScale = require( './fixtures/julia/large_scale.json' );
+var smallX = require( './fixtures/julia/small_x.json' );
 var largeShape = require( './fixtures/julia/large_shape.json' );
 var bothLarge = require( './fixtures/julia/both_large.json' );
 
@@ -142,7 +143,7 @@ tape( 'the function evaluates the logcdf for `x` given large `lambda` and `k`',
 	k = bothLarge.k;
 	for ( i = 0; i < x.length; i++ ) {
 		y = logcdf( x[i], k[i], lambda[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();
 });
@@ -161,7 +162,7 @@ tape( 'the function evaluates the logcdf for `x` given large shape parameter `la
 	k = largeShape.k;
 	for ( i = 0; i < x.length; i++ ) {
 		y = logcdf( x[i], k[i], lambda[i] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 0 ), true, 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 3236 ), true, 'returns expected value' );
 	}
 	t.end();
 });
@@ -180,7 +181,26 @@ tape( 'the function evaluates the logcdf for `x` given large scale parameter `k`
 	k = largeScale.k;
 	for ( i = 0; i < x.length; i++ ) {
 		y = logcdf( x[i], k[i], lambda[i] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 0 ), true, 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 2 ), true, 'returns expected value' );
+	}
+	t.end();
+});
+
+tape( 'the function evaluates the logcdf for `x` close to `0`', opts, function test( t ) {
+	var expected;
+	var lambda;
+	var i;
+	var k;
+	var x;
+	var y;
+
+	expected = smallX.expected;
+	x = smallX.x;
+	lambda = smallX.lambda;
+	k = smallX.k;
+	for ( i = 0; i < x.length; i++ ) {
+		y = logcdf( x[i], k[i], lambda[i] );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 1 ), true, 'returns expected value' );
 	}
 	t.end();
 });