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

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

@@ -29,7 +29,10 @@ 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 LN2 = require( '@stdlib/constants/float64/ln-two' );
+var FLOAT64_SMALLEST_NORMAL = require( '@stdlib/constants/float64/smallest-normal' );
 var NINF = require( '@stdlib/constants/float64/ninf' );
+var PINF = require( '@stdlib/constants/float64/pinf' );
 
 
 // MAIN //
@@ -71,14 +74,33 @@ function factory( sigma ) {
 	* // returns <number>
 	*/
 	function logcdf( x ) {
+		var lt;
+		var x2;
+		var r;
 		var p;
 		if ( isnan( x ) ) {
 			return NaN;
 		}
 		if ( x < 0.0 ) {
 			return NINF;
 		}
-		p = -pow( x, 2.0 ) / ( 2.0 * s2 );
+		r = x / sigma;
+		if ( r < FLOAT64_SMALLEST_NORMAL ) {
+			lt = ( 2.0*( ln( x ) - ln( sigma ) ) ) - LN2;
+		} else {
+			lt = ( 2.0*ln( r ) ) - LN2;
+		}
+		if ( lt < -36.0 ) {
+			// For `t = x^2/(2*sigma^2) < 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 `x^2` underflowing to `0`:
+			return lt;
+		}
+		x2 = pow( x, 2.0 );
+		if ( s2 >= FLOAT64_SMALLEST_NORMAL && 2.0 * s2 < PINF && x2 < PINF ) {
+			p = -x2 / ( 2.0 * s2 );
+		} else {
+			// Where `sigma^2` or `x^2` would overflow or underflow, take `t` from the ratio `r` instead; elsewhere the direct form is kept, as `r*r` rounds slightly worse for large `r`...
+			p = -0.5 * r * r;
+		}
 		return ( p < LNHALF ) ? log1p( -exp( p ) ) : ln( -expm1( p ) );
 	}
 }

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

@@ -27,7 +27,10 @@ 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 LN2 = require( '@stdlib/constants/float64/ln-two' );
+var FLOAT64_SMALLEST_NORMAL = require( '@stdlib/constants/float64/smallest-normal' );
 var NINF = require( '@stdlib/constants/float64/ninf' );
+var PINF = require( '@stdlib/constants/float64/pinf' );
 
 
 // MAIN //
@@ -66,6 +69,9 @@ var NINF = require( '@stdlib/constants/float64/ninf' );
 */
 function logcdf( x, sigma ) {
 	var s2;
+	var x2;
+	var lt;
+	var r;
 	var p;
 	if (
 		isnan( x ) ||
@@ -81,7 +87,23 @@ function logcdf( x, sigma ) {
 		return NINF;
 	}
 	s2 = pow( sigma, 2.0 );
-	p = -pow( x, 2.0 ) / ( 2.0 * s2 );
+	r = x / sigma;
+	if ( r < FLOAT64_SMALLEST_NORMAL ) {
+		lt = ( 2.0*( ln( x ) - ln( sigma ) ) ) - LN2;
+	} else {
+		lt = ( 2.0*ln( r ) ) - LN2;
+	}
+	if ( lt < -36.0 ) {
+		// For `t = x^2/(2*sigma^2) < 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 `x^2` underflowing to `0`:
+		return lt;
+	}
+	x2 = pow( x, 2.0 );
+	if ( s2 >= FLOAT64_SMALLEST_NORMAL && 2.0 * s2 < PINF && x2 < PINF ) {
+		p = -x2 / ( 2.0 * s2 );
+	} else {
+		// Where `sigma^2` or `x^2` would overflow or underflow, take `t` from the ratio `r` instead; elsewhere the direct form is kept, as `r*r` rounds slightly worse for large `r`...
+		p = -0.5 * r * r;
+	}
 	return ( p < LNHALF ) ? log1p( -exp( p ) ) : ln( -expm1( p ) );
 }
 

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

@@ -46,7 +46,10 @@
         "@stdlib/math/base/special/pow",
         "@stdlib/math/base/special/ln",
         "@stdlib/constants/float64/ln-half",
-        "@stdlib/constants/float64/ninf"
+        "@stdlib/constants/float64/ln-two",
+        "@stdlib/constants/float64/smallest-normal",
+        "@stdlib/constants/float64/ninf",
+        "@stdlib/constants/float64/pinf"
       ]
     },
     {
@@ -69,7 +72,10 @@
         "@stdlib/math/base/special/pow",
         "@stdlib/math/base/special/ln",
         "@stdlib/constants/float64/ln-half",
-        "@stdlib/constants/float64/ninf"
+        "@stdlib/constants/float64/ln-two",
+        "@stdlib/constants/float64/smallest-normal",
+        "@stdlib/constants/float64/ninf",
+        "@stdlib/constants/float64/pinf"
       ]
     },
     {
@@ -91,7 +97,10 @@
         "@stdlib/math/base/special/pow",
         "@stdlib/math/base/special/ln",
         "@stdlib/constants/float64/ln-half",
-        "@stdlib/constants/float64/ninf"
+        "@stdlib/constants/float64/ln-two",
+        "@stdlib/constants/float64/smallest-normal",
+        "@stdlib/constants/float64/ninf",
+        "@stdlib/constants/float64/pinf"
       ]
     }
   ]

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

@@ -24,7 +24,10 @@
 #include "stdlib/math/base/special/pow.h"
 #include "stdlib/math/base/special/ln.h"
 #include "stdlib/constants/float64/ln_half.h"
+#include "stdlib/constants/float64/ln_two.h"
+#include "stdlib/constants/float64/smallest_normal.h"
 #include "stdlib/constants/float64/ninf.h"
+#include "stdlib/constants/float64/pinf.h"
 
 /**
 * Evaluates the natural logarithm of the cumulative distribution function (CDF) for a Rayleigh distribution with scale parameter `sigma` at a value `x`.
@@ -39,6 +42,9 @@
 */
 double stdlib_base_dists_rayleigh_logcdf( const double x, const double sigma ) {
 	double s2;
+	double x2;
+	double lt;
+	double r;
 	double p;
 	if (
 		stdlib_base_is_nan( x ) ||
@@ -54,6 +60,22 @@ double stdlib_base_dists_rayleigh_logcdf( const double x, const double sigma ) {
 		return STDLIB_CONSTANT_FLOAT64_NINF;
 	}
 	s2 = stdlib_base_pow( sigma, 2.0 );
-	p = -stdlib_base_pow( x, 2.0 ) / ( 2.0 * s2 );
+	r = x / sigma;
+	if ( r < STDLIB_CONSTANT_FLOAT64_SMALLEST_NORMAL ) {
+		lt = ( 2.0 * ( stdlib_base_ln( x ) - stdlib_base_ln( sigma ) ) ) - STDLIB_CONSTANT_FLOAT64_LN2;
+	} else {
+		lt = ( 2.0 * stdlib_base_ln( r ) ) - STDLIB_CONSTANT_FLOAT64_LN2;
+	}
+	if ( lt < -36.0 ) {
+		// For `t = x^2/(2*sigma^2) < 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 `x^2` underflowing to `0`:
+		return lt;
+	}
+	x2 = stdlib_base_pow( x, 2.0 );
+	if ( s2 >= STDLIB_CONSTANT_FLOAT64_SMALLEST_NORMAL && 2.0 * s2 < STDLIB_CONSTANT_FLOAT64_PINF && x2 < STDLIB_CONSTANT_FLOAT64_PINF ) {
+		p = -x2 / ( 2.0 * s2 );
+	} else {
+		// Where `sigma^2` or `x^2` would overflow or underflow, take `t` from the ratio `r` instead; elsewhere the direct form is kept, as `r*r` rounds slightly worse for large `r`...
+		p = -0.5 * r * r;
+	}
 	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/rayleigh/logcdf/test/fixtures/julia/runner.jl

@@ -61,6 +61,39 @@ function gen( x, sigma, name )
 	close( outfile );
 end
 
+"""
+	gen_small( x, sigma, name )
+
+Generate fixture data for small `x` in extended precision and write to file. In double precision, `x^2` underflows for small `x`, so the reference value is computed from the textbook formula `ln(1 - exp(-x^2/(2*sigma^2)))`, written with `expm1` so that it stays exact at 2048 bits.
+
+# Arguments
+
+* `x`: input value
+* `sigma`: scale parameter
+* `name::AbstractString`: output filename
+"""
+function gen_small( x, sigma, name )
+	z = setprecision( BigFloat, 2048 ) do
+		Float64.( log.( .-expm1.( .-( BigFloat.( x ) .^ 2 ) ./ ( 2 .* BigFloat.( sigma ) .^ 2 ) ) ) )
+	end
+
+	# Store data to be written to file as a collection:
+	data = Dict([
+		("x", x),
+		("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__;
 
@@ -81,3 +114,8 @@ gen( x, sigma, "medium_scale.json" );
 x = rand( 1000 ) .* 20.0;
 sigma = rand( 1000 ) .* 5.0;
 gen( x, sigma, "small_scale.json" );
+
+# Small `x` (relative to `sigma`, over a wide range of `sigma`):
+sigma = exp10.( ( rand( 1000 ) .* 200.0 ) .- 100.0 );
+x = sigma .* exp10.( range( -150.0, stop = -4.0, length = 1000 ) );
+gen_small( x, sigma, "small_x.json" );

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

@@ -31,6 +31,7 @@ var factory = require( './../lib/factory.js' );
 // FIXTURES //
 
 var smallScale = require( './fixtures/julia/small_scale.json' );
+var smallX = require( './fixtures/julia/small_x.json' );
 var mediumScale = require( './fixtures/julia/medium_scale.json' );
 var largeScale = require( './fixtures/julia/large_scale.json' );
 
@@ -123,6 +124,30 @@ tape( 'if `sigma` equals `0`, the created function evaluates a degenerate distri
 	t.end();
 });
 
+tape( 'the created function evaluates the logcdf when `sigma^2` or `x^2` would overflow or underflow', function test( t ) {
+	var cases;
+	var c;
+	var y;
+	var i;
+
+	// `sigma^2` or `x^2` underflows or overflows for each of these:
+	cases = [
+		[ 1.0e-170, 1.0e-170, -0.9327521295671886 ],
+		[ 1.0e200, 1.0e200, -0.9327521295671886 ],
+		[
+			1.0847105013290125e-163,
+			1.6127895670863117e-161,
+			-10.696803256634931
+		]
+	];
+	for ( i = 0; i < cases.length; i++ ) {
+		c = cases[ i ];
+		y = factory( c[ 1 ] )( c[ 0 ] );
+		t.strictEqual( isAlmostSameValue( y, c[ 2 ], 1 ), true, 'returns expected value' );
+	}
+	t.end();
+});
+
 tape( 'the created function evaluates the logcdf for `x` given small scale parameter `sigma`', function test( t ) {
 	var expected;
 	var logcdf;
@@ -179,3 +204,22 @@ tape( 'the created function evaluates the logcdf for `x` given large scale param
 	}
 	t.end();
 });
+
+tape( 'the created function evaluates the logcdf for `x` close to `0`', function test( t ) {
+	var expected;
+	var logcdf;
+	var sigma;
+	var i;
+	var x;
+	var y;
+
+	expected = smallX.expected;
+	x = smallX.x;
+	sigma = smallX.sigma;
+	for ( i = 0; i < x.length; i++ ) {
+		logcdf = factory( sigma[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/rayleigh/logcdf/test/test.logcdf.js

@@ -31,6 +31,7 @@ var logcdf = require( './../lib' );
 // FIXTURES //
 
 var smallScale = require( './fixtures/julia/small_scale.json' );
+var smallX = require( './fixtures/julia/small_x.json' );
 var mediumScale = require( './fixtures/julia/medium_scale.json' );
 var largeScale = require( './fixtures/julia/large_scale.json' );
 
@@ -90,6 +91,30 @@ tape( 'if provided `sigma` equals `0`, the function evaluates a degenerate distr
 	t.end();
 });
 
+tape( 'the function evaluates the logcdf when `sigma^2` or `x^2` would overflow or underflow', function test( t ) {
+	var cases;
+	var c;
+	var y;
+	var i;
+
+	// `sigma^2` or `x^2` underflows or overflows for each of these:
+	cases = [
+		[ 1.0e-170, 1.0e-170, -0.9327521295671886 ],
+		[ 1.0e200, 1.0e200, -0.9327521295671886 ],
+		[
+			1.0847105013290125e-163,
+			1.6127895670863117e-161,
+			-10.696803256634931
+		]
+	];
+	for ( i = 0; i < cases.length; i++ ) {
+		c = cases[ i ];
+		y = logcdf( c[ 0 ], c[ 1 ] );
+		t.strictEqual( isAlmostSameValue( y, c[ 2 ], 1 ), true, 'returns expected value' );
+	}
+	t.end();
+});
+
 tape( 'the function evaluates the logcdf for `x` given small scale parameter `sigma`', function test( t ) {
 	var expected;
 	var sigma;
@@ -140,3 +165,20 @@ tape( 'the function evaluates the logcdf for `x` given large scale parameter `si
 	}
 	t.end();
 });
+
+tape( 'the function evaluates the logcdf for `x` close to `0`', function test( t ) {
+	var expected;
+	var sigma;
+	var i;
+	var x;
+	var y;
+
+	expected = smallX.expected;
+	x = smallX.x;
+	sigma = smallX.sigma;
+	for ( i = 0; i < x.length; i++ ) {
+		y = logcdf( x[i], sigma[i] );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 1 ), true, 'returns expected value' );
+	}
+	t.end();
+});