fix: avoid overflow and underflow for scales far from 1 in stats/base/dists/rayleigh/logpdf

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

@@ -27,6 +27,7 @@ var pow = require( '@stdlib/math/base/special/pow' );
 var ln = require( '@stdlib/math/base/special/ln' );
 var PINF = require( '@stdlib/constants/float64/pinf' );
 var NINF = require( '@stdlib/constants/float64/ninf' );
+var FLOAT64_SMALLEST_NORMAL = require( '@stdlib/constants/float64/smallest-normal' );
 
 
 // MAIN //
@@ -46,7 +47,9 @@ var NINF = require( '@stdlib/constants/float64/ninf' );
 * // returns ~-0.936
 */
 function factory( sigma ) {
+	var normal;
 	var s2i;
+	var lns;
 	var s2;
 	if ( isnan( sigma ) || sigma < 0.0 ) {
 		return constantFunction( NaN );
@@ -56,6 +59,8 @@ function factory( sigma ) {
 	}
 	s2 = pow( sigma, 2.0 );
 	s2i = 1.0 / s2;
+	normal = ( s2 >= FLOAT64_SMALLEST_NORMAL && 2.0 * s2 < PINF );
+	lns = ln( sigma );
 	return logpdf;
 
 	/**
@@ -70,13 +75,32 @@ function factory( sigma ) {
 	* // returns <number>
 	*/
 	function logpdf( x ) {
+		var lq;
+		var x2;
+		var r;
+		var q;
 		if ( isnan( x ) ) {
 			return NaN;
 		}
 		if ( x < 0.0 || x === PINF ) {
 			return NINF;
 		}
-		return ln( s2i * x ) - (pow( x, 2.0 ) / ( 2.0 * s2 ));
+		x2 = pow( x, 2.0 );
+		q = s2i * x;
+		if ( normal && x2 < PINF && q >= FLOAT64_SMALLEST_NORMAL ) {
+			return ln( q ) - (x2 / ( 2.0 * s2 ));
+		}
+		// For `sigma` or `x` far from `1`, forming `sigma^2` or `x^2` overflows or underflows, so evaluate via the ratio `x/sigma` instead. `ln( x/sigma^2 )` is taken directly while that ratio is normal, as subtracting two large logarithms cancels when `x` is close to `sigma^2`...
+		r = x / sigma;
+		q = r / sigma;
+		if ( q >= FLOAT64_SMALLEST_NORMAL && q < PINF ) {
+			lq = ln( q );
+		} else if ( r >= FLOAT64_SMALLEST_NORMAL ) {
+			lq = ln( r ) - lns;
+		} else {
+			lq = ( ln( x ) - lns ) - lns;
+		}
+		return lq - ( 0.5 * r * r );
 	}
 }
 

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

@@ -25,6 +25,7 @@ var pow = require( '@stdlib/math/base/special/pow' );
 var ln = require( '@stdlib/math/base/special/ln' );
 var PINF = require( '@stdlib/constants/float64/pinf' );
 var NINF = require( '@stdlib/constants/float64/ninf' );
+var FLOAT64_SMALLEST_NORMAL = require( '@stdlib/constants/float64/smallest-normal' );
 
 
 // MAIN //
@@ -63,7 +64,11 @@ var NINF = require( '@stdlib/constants/float64/ninf' );
 */
 function logpdf( x, sigma ) {
 	var s2i;
+	var lq;
 	var s2;
+	var x2;
+	var r;
+	var q;
 	if (
 		isnan( x ) ||
 		isnan( sigma ) ||
@@ -79,7 +84,27 @@ function logpdf( x, sigma ) {
 	}
 	s2 = pow( sigma, 2.0 );
 	s2i = 1.0 / s2;
-	return ln( s2i * x ) - ( pow( x, 2.0 ) / ( 2.0 * s2 ) );
+	x2 = pow( x, 2.0 );
+	q = s2i * x;
+	if (
+		s2 >= FLOAT64_SMALLEST_NORMAL &&
+		2.0 * s2 < PINF &&
+		x2 < PINF &&
+		q >= FLOAT64_SMALLEST_NORMAL
+	) {
+		return ln( q ) - ( x2 / ( 2.0 * s2 ) );
+	}
+	// For `sigma` or `x` far from `1`, forming `sigma^2` or `x^2` overflows or underflows, so evaluate via the ratio `x/sigma` instead. `ln( x/sigma^2 )` is taken directly while that ratio is normal, as subtracting two large logarithms cancels when `x` is close to `sigma^2`...
+	r = x / sigma;
+	q = r / sigma;
+	if ( q >= FLOAT64_SMALLEST_NORMAL && q < PINF ) {
+		lq = ln( q );
+	} else if ( r >= FLOAT64_SMALLEST_NORMAL ) {
+		lq = ln( r ) - ln( sigma );
+	} else {
+		lq = ( ln( x ) - ln( sigma ) ) - ln( sigma );
+	}
+	return lq - ( 0.5 * r * r );
 }
 
 

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

@@ -43,7 +43,8 @@
         "@stdlib/math/base/special/pow",
         "@stdlib/math/base/special/ln",
         "@stdlib/constants/float64/pinf",
-        "@stdlib/constants/float64/ninf"
+        "@stdlib/constants/float64/ninf",
+        "@stdlib/constants/float64/smallest-normal"
       ]
     },
     {
@@ -63,7 +64,8 @@
         "@stdlib/math/base/special/pow",
         "@stdlib/math/base/special/ln",
         "@stdlib/constants/float64/pinf",
-        "@stdlib/constants/float64/ninf"
+        "@stdlib/constants/float64/ninf",
+        "@stdlib/constants/float64/smallest-normal"
       ]
     },
     {
@@ -82,7 +84,8 @@
         "@stdlib/math/base/special/pow",
         "@stdlib/math/base/special/ln",
         "@stdlib/constants/float64/pinf",
-        "@stdlib/constants/float64/ninf"
+        "@stdlib/constants/float64/ninf",
+        "@stdlib/constants/float64/smallest-normal"
       ]
     }
   ]

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

@@ -22,6 +22,7 @@
 #include "stdlib/math/base/special/ln.h"
 #include "stdlib/constants/float64/pinf.h"
 #include "stdlib/constants/float64/ninf.h"
+#include "stdlib/constants/float64/smallest_normal.h"
 
 /**
 * Evaluates the natural logarithm of the probability density function (PDF) for a Rayleigh distribution with scale parameter `sigma` at a value `x`.
@@ -36,7 +37,11 @@
 */
 double stdlib_base_dists_rayleigh_logpdf( const double x, const double sigma ) {
 	double s2i;
+	double lq;
 	double s2;
+	double x2;
+	double r;
+	double q;
 	if (
 		stdlib_base_is_nan( x ) ||
 		stdlib_base_is_nan( sigma ) ||
@@ -52,5 +57,25 @@ double stdlib_base_dists_rayleigh_logpdf( const double x, const double sigma ) {
 	}
 	s2 = stdlib_base_pow( sigma, 2.0 );
 	s2i = 1.0 / s2;
-	return stdlib_base_ln( s2i * x ) - stdlib_base_pow( x, 2.0 ) / ( 2.0 * s2 );
+	x2 = stdlib_base_pow( x, 2.0 );
+	q = s2i * x;
+	if (
+		s2 >= STDLIB_CONSTANT_FLOAT64_SMALLEST_NORMAL &&
+		2.0 * s2 < STDLIB_CONSTANT_FLOAT64_PINF &&
+		x2 < STDLIB_CONSTANT_FLOAT64_PINF &&
+		q >= STDLIB_CONSTANT_FLOAT64_SMALLEST_NORMAL
+	) {
+		return stdlib_base_ln( q ) - x2 / ( 2.0 * s2 );
+	}
+	// For `sigma` or `x` far from `1`, forming `sigma^2` or `x^2` overflows or underflows, so evaluate via the ratio `x/sigma` instead. `ln( x/sigma^2 )` is taken directly while that ratio is normal, as subtracting two large logarithms cancels when `x` is close to `sigma^2`...
+	r = x / sigma;
+	q = r / sigma;
+	if ( q >= STDLIB_CONSTANT_FLOAT64_SMALLEST_NORMAL && q < STDLIB_CONSTANT_FLOAT64_PINF ) {
+		lq = stdlib_base_ln( q );
+	} else if ( r >= STDLIB_CONSTANT_FLOAT64_SMALLEST_NORMAL ) {
+		lq = stdlib_base_ln( r ) - stdlib_base_ln( sigma );
+	} else {
+		lq = ( stdlib_base_ln( x ) - stdlib_base_ln( sigma ) ) - stdlib_base_ln( sigma );
+	}
+	return lq - ( 0.5 * r * r );
 }

lib/node_modules/@stdlib/stats/base/dists/rayleigh/logpdf/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: logpdf, Rayleigh
 import JSON
 
 """
@@ -39,9 +38,11 @@ julia> gen( x, sigma, \"data.json\" );
 ```
 """
 function gen( x, sigma, name )
-	z = Array{Float64}( undef, length(x) );
-	for i in eachindex(x)
-		z[ i ] = logpdf( Rayleigh( sigma[i] ), x[i] );
+	# The log PDF is evaluated in extended precision using the textbook formula `ln(x/sigma^2) - x^2/(2*sigma^2)`, as forming `sigma^2` or `x^2` in double precision overflows or underflows for scale parameters far from `1`...
+	z = setprecision( BigFloat, 2048 ) do
+		X = BigFloat.( x );
+		S = BigFloat.( sigma );
+		Float64.( log.( X ./ ( S .^ 2 ) ) .- ( ( X .^ 2 ) ./ ( 2 .* ( S .^ 2 ) ) ) );
 	end
 
 	# Store data to be written to file as a collection:
@@ -81,3 +82,13 @@ gen( x, sigma, "medium_scale.json" );
 x = rand( 1000 ) .* 20.0;
 sigma = rand( 1000 ) .* 5.0;
 gen( x, sigma, "small_scale.json" );
+
+# Tiny scale parameter (`sigma^2` underflows):
+sigma = exp10.( ( rand( 1000 ) .* 45.0 ) .- 200.0 );
+x = sigma .* exp10.( ( rand( 1000 ) .* 4.0 ) .- 3.0 );
+gen( x, sigma, "tiny_scale.json" );
+
+# Huge scale parameter (`sigma^2` overflows):
+sigma = exp10.( ( rand( 1000 ) .* 45.0 ) .+ 155.0 );
+x = sigma .* exp10.( ( rand( 1000 ) .* 4.0 ) .- 3.0 );
+gen( x, sigma, "huge_scale.json" );

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

@@ -32,6 +32,8 @@ var factory = require( './../lib/factory.js' );
 
 var smallScale = require( './fixtures/julia/small_scale.json' );
 var mediumScale = require( './fixtures/julia/medium_scale.json' );
+var tinyScale = require( './fixtures/julia/tiny_scale.json' );
+var hugeScale = require( './fixtures/julia/huge_scale.json' );
 var largeScale = require( './fixtures/julia/large_scale.json' );
 
 
@@ -187,3 +189,45 @@ tape( 'the created function evaluates the logpdf for `x` given large scale param
 	}
 	t.end();
 });
+
+tape( 'the created function evaluates the logpdf for `x` given tiny scale parameter `sigma`', function test( t ) {
+	var expected;
+	var logpdf;
+	var sigma;
+	var i;
+	var x;
+	var y;
+
+	expected = tinyScale.expected;
+	x = tinyScale.x;
+	sigma = tinyScale.sigma;
+	for ( i = 0; i < x.length; i++ ) {
+		logpdf = factory( sigma[i] );
+		y = logpdf( x[i] );
+		if ( expected[i] !== null ) {
+			t.strictEqual( isAlmostSameValue( y, expected[i], 2 ), true, 'returns expected value' );
+		}
+	}
+	t.end();
+});
+
+tape( 'the created function evaluates the logpdf for `x` given huge scale parameter `sigma`', function test( t ) {
+	var expected;
+	var logpdf;
+	var sigma;
+	var i;
+	var x;
+	var y;
+
+	expected = hugeScale.expected;
+	x = hugeScale.x;
+	sigma = hugeScale.sigma;
+	for ( i = 0; i < x.length; i++ ) {
+		logpdf = factory( sigma[i] );
+		y = logpdf( x[i] );
+		if ( expected[i] !== null ) {
+			t.strictEqual( isAlmostSameValue( y, expected[i], 2 ), true, 'returns expected value' );
+		}
+	}
+	t.end();
+});

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

@@ -33,6 +33,8 @@ var logpdf = require( './../lib' );
 
 var smallScale = require( './fixtures/julia/small_scale.json' );
 var mediumScale = require( './fixtures/julia/medium_scale.json' );
+var tinyScale = require( './fixtures/julia/tiny_scale.json' );
+var hugeScale = require( './fixtures/julia/huge_scale.json' );
 var largeScale = require( './fixtures/julia/large_scale.json' );
 
 
@@ -180,3 +182,57 @@ tape( 'the function evaluates the logpdf for `x` given large scale parameter `si
 	}
 	t.end();
 });
+
+tape( 'the function evaluates the logpdf for `x` given tiny scale parameter `sigma`', function test( t ) {
+	var expected;
+	var delta;
+	var sigma;
+	var tol;
+	var i;
+	var x;
+	var y;
+
+	expected = tinyScale.expected;
+	x = tinyScale.x;
+	sigma = tinyScale.sigma;
+	for ( i = 0; i < x.length; i++ ) {
+		y = logpdf( x[i], sigma[i] );
+		if ( expected[i] !== null ) {
+			if ( y === expected[i] ) {
+				t.strictEqual( y, expected[i], 'x: '+x[i]+', sigma: '+sigma[i]+', y: '+y+', expected: '+expected[i] );
+			} else {
+				delta = abs( y - expected[ i ] );
+				tol = 3.0 * EPS * abs( expected[ i ] );
+				t.ok( delta <= tol, 'within tolerance. x: '+x[ i ]+'. sigma: '+sigma[i]+'. y: '+y+'. E: '+expected[ i ]+'. Δ: '+delta+'. tol: '+tol+'.' );
+			}
+		}
+	}
+	t.end();
+});
+
+tape( 'the function evaluates the logpdf for `x` given huge scale parameter `sigma`', function test( t ) {
+	var expected;
+	var delta;
+	var sigma;
+	var tol;
+	var i;
+	var x;
+	var y;
+
+	expected = hugeScale.expected;
+	x = hugeScale.x;
+	sigma = hugeScale.sigma;
+	for ( i = 0; i < x.length; i++ ) {
+		y = logpdf( x[i], sigma[i] );
+		if ( expected[i] !== null ) {
+			if ( y === expected[i] ) {
+				t.strictEqual( y, expected[i], 'x: '+x[i]+', sigma: '+sigma[i]+', y: '+y+', expected: '+expected[i] );
+			} else {
+				delta = abs( y - expected[ i ] );
+				tol = 3.0 * EPS * abs( expected[ i ] );
+				t.ok( delta <= tol, 'within tolerance. x: '+x[ i ]+'. sigma: '+sigma[i]+'. y: '+y+'. E: '+expected[ i ]+'. Δ: '+delta+'. tol: '+tol+'.' );
+			}
+		}
+	}
+	t.end();
+});

lib/node_modules/@stdlib/stats/base/dists/rayleigh/logpdf/test/test.native.js

@@ -33,6 +33,8 @@ var NINF = require( '@stdlib/constants/float64/ninf' );
 
 var smallScale = require( './fixtures/julia/small_scale.json' );
 var mediumScale = require( './fixtures/julia/medium_scale.json' );
+var tinyScale = require( './fixtures/julia/tiny_scale.json' );
+var hugeScale = require( './fixtures/julia/huge_scale.json' );
 var largeScale = require( './fixtures/julia/large_scale.json' );
 
 
@@ -164,3 +166,41 @@ tape( 'the function evaluates the logpdf for `x` given large scale parameter `si
 	}
 	t.end();
 });
+
+tape( 'the function evaluates the logpdf for `x` given tiny scale parameter `sigma`', opts, function test( t ) {
+	var expected;
+	var sigma;
+	var i;
+	var x;
+	var y;
+
+	expected = tinyScale.expected;
+	x = tinyScale.x;
+	sigma = tinyScale.sigma;
+	for ( i = 0; i < x.length; i++ ) {
+		y = logpdf( x[i], sigma[i] );
+		if ( expected[i] !== null ) {
+			t.strictEqual( isAlmostSameValue( y, expected[i], 2 ), true, 'returns expected value' );
+		}
+	}
+	t.end();
+});
+
+tape( 'the function evaluates the logpdf for `x` given huge scale parameter `sigma`', opts, function test( t ) {
+	var expected;
+	var sigma;
+	var i;
+	var x;
+	var y;
+
+	expected = hugeScale.expected;
+	x = hugeScale.x;
+	sigma = hugeScale.sigma;
+	for ( i = 0; i < x.length; i++ ) {
+		y = logpdf( x[i], sigma[i] );
+		if ( expected[i] !== null ) {
+			t.strictEqual( isAlmostSameValue( y, expected[i], 2 ), true, 'returns expected value' );
+		}
+	}
+	t.end();
+});