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

`(x/lambda)^k` underflows to `0` for small `x`, so the function returned
`-Infinity` although the log-CDF is `k*ln(x/lambda)` to working
precision there. Return `k*ln(x/lambda)` once it is below `-36`. The
fixtures are now generated from the textbook formula in extended
precision.

Author: Abhist17

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

@@ -71,13 +71,19 @@ function factory( k, lambda ) {
 	* // returns <number>
 	*/
 	function logcdf( x ) {
+		var lt;
 		var p;
 		if ( isnan( x ) ) {
 			return NaN;
 		}
 		if ( x < 0.0 ) {
 			return NINF;
 		}
+		lt = k * ln( x / lambda );
+		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

@@ -73,6 +73,7 @@ var NINF = require( '@stdlib/constants/float64/ninf' );
 * // returns NaN
 */
 function logcdf( x, k, lambda ) {
+	var lt;
 	var p;
 	if (
 		isnan( k ) ||
@@ -85,6 +86,11 @@ function logcdf( x, k, lambda ) {
 	if ( x < 0.0 ) {
 		return NINF;
 	}
+	lt = k * ln( x / lambda );
+	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/src/main.c

@@ -39,6 +39,7 @@
 * // returns ~-0.145
 */
 double stdlib_base_dists_weibull_logcdf( const double x, const double k, const double lambda ) {
+	double lt;
 	double p;
 	if (
 		stdlib_base_is_nan( x ) ||
@@ -51,6 +52,11 @@ double stdlib_base_dists_weibull_logcdf( const double x, const double k, const d
 	if ( x < 0.0 ) {
 		return STDLIB_CONSTANT_FLOAT64_NINF;
 	}
+	lt = k * stdlib_base_ln( x / lambda );
+	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/both_large.json

@@ -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( -300.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/fixtures/julia/large_scale.json

@@ -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();
 });