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

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

@@ -25,7 +25,9 @@ var degenerate = require( '@stdlib/stats/base/dists/degenerate/pdf' ).factory;
 var isnan = require( '@stdlib/math/base/assert/is-nan' );
 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 PINF = require( '@stdlib/constants/float64/pinf' );
+var FLOAT64_SMALLEST_NORMAL = require( '@stdlib/constants/float64/smallest-normal' );
 
 
 // MAIN //
@@ -45,7 +47,9 @@ var PINF = require( '@stdlib/constants/float64/pinf' );
 * // returns ~0.392
 */
 function factory( sigma ) {
+	var normal;
 	var s2i;
+	var lns;
 	var s2;
 	if ( isnan( sigma ) || sigma < 0.0 ) {
 		return constantFunction( NaN );
@@ -55,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 pdf;
 
 	/**
@@ -69,13 +75,31 @@ function factory( sigma ) {
 	* // returns <number>
 	*/
 	function pdf( x ) {
+		var x2;
+		var r;
+		var q;
 		if ( isnan( x ) ) {
 			return NaN;
 		}
 		if ( x < 0.0 || x === PINF ) {
 			return 0.0;
 		}
-		return s2i * x * exp( -pow( x, 2.0 ) / ( 2.0 * s2 ) );
+		x2 = pow( x, 2.0 );
+		q = s2i * x;
+		if ( normal && x2 < PINF && q >= FLOAT64_SMALLEST_NORMAL ) {
+			return q * exp( -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...
+		r = x / sigma;
+		if ( r === PINF ) {
+			return 0.0;
+		}
+		q = r / sigma;
+		if ( q < PINF ) {
+			return q * exp( -0.5 * r * r );
+		}
+		// `x/sigma^2` overflows only for a tiny `sigma`, where the exponential decides the result:
+		return exp( ( ln( r ) - lns ) - ( 0.5 * r * r ) );
 	}
 }
 

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

@@ -23,7 +23,9 @@
 var isnan = require( '@stdlib/math/base/assert/is-nan' );
 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 PINF = require( '@stdlib/constants/float64/pinf' );
+var FLOAT64_SMALLEST_NORMAL = require( '@stdlib/constants/float64/smallest-normal' );
 
 
 // MAIN //
@@ -63,6 +65,9 @@ var PINF = require( '@stdlib/constants/float64/pinf' );
 function pdf( x, sigma ) {
 	var s2i;
 	var s2;
+	var x2;
+	var r;
+	var q;
 	if (
 		isnan( x ) ||
 		isnan( sigma ) ||
@@ -78,7 +83,27 @@ function pdf( x, sigma ) {
 	}
 	s2 = pow( sigma, 2.0 );
 	s2i = 1.0 / s2;
-	return s2i * x * exp( -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 q * exp( -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...
+	r = x / sigma;
+	if ( r === PINF ) {
+		return 0.0;
+	}
+	q = r / sigma;
+	if ( q < PINF ) {
+		return q * exp( -0.5 * r * r );
+	}
+	// `x/sigma^2` overflows only for a tiny `sigma`, where the exponential decides the result:
+	return exp( ( ln( r ) - ln( sigma ) ) - ( 0.5 * r * r ) );
 }
 
 

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

@@ -42,7 +42,9 @@
         "@stdlib/math/base/assert/is-nan",
         "@stdlib/math/base/special/exp",
         "@stdlib/math/base/special/pow",
-        "@stdlib/constants/float64/pinf"
+        "@stdlib/constants/float64/pinf",
+        "@stdlib/constants/float64/smallest-normal",
+        "@stdlib/math/base/special/ln"
       ]
     },
     {
@@ -61,7 +63,9 @@
         "@stdlib/math/base/special/exp",
         "@stdlib/constants/float64/eps",
         "@stdlib/math/base/special/pow",
-        "@stdlib/constants/float64/pinf"
+        "@stdlib/constants/float64/pinf",
+        "@stdlib/constants/float64/smallest-normal",
+        "@stdlib/math/base/special/ln"
       ]
     },
     {
@@ -79,7 +83,9 @@
         "@stdlib/math/base/assert/is-nan",
         "@stdlib/math/base/special/exp",
         "@stdlib/math/base/special/pow",
-        "@stdlib/constants/float64/pinf"
+        "@stdlib/constants/float64/pinf",
+        "@stdlib/constants/float64/smallest-normal",
+        "@stdlib/math/base/special/ln"
       ]
     }
   ]

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

@@ -20,7 +20,9 @@
 #include "stdlib/math/base/assert/is_nan.h"
 #include "stdlib/math/base/special/exp.h"
 #include "stdlib/math/base/special/pow.h"
+#include "stdlib/math/base/special/ln.h"
 #include "stdlib/constants/float64/pinf.h"
+#include "stdlib/constants/float64/smallest_normal.h"
 
 /**
 * Evaluates the probability density function (PDF) for a Rayleigh distribution with scale parameter `sigma` at a value `x`.
@@ -36,6 +38,9 @@
 double stdlib_base_dists_rayleigh_pdf( const double x, const double sigma ) {
 	double s2;
 	double s2i;
+	double x2;
+	double r;
+	double q;
 	if (
 		stdlib_base_is_nan( x ) ||
 		stdlib_base_is_nan( sigma ) ||
@@ -51,5 +56,25 @@ double stdlib_base_dists_rayleigh_pdf( const double x, const double sigma ) {
 	}
 	s2 = stdlib_base_pow( sigma, 2.0 );
 	s2i = 1.0 / s2;
-	return s2i * x * stdlib_base_exp( -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 q * stdlib_base_exp( -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...
+	r = x / sigma;
+	if ( r == STDLIB_CONSTANT_FLOAT64_PINF ) {
+		return 0.0;
+	}
+	q = r / sigma;
+	if ( q < STDLIB_CONSTANT_FLOAT64_PINF ) {
+		return q * stdlib_base_exp( -0.5 * r * r );
+	}
+	// `x/sigma^2` overflows only for a tiny `sigma`, where the exponential decides the result:
+	return stdlib_base_exp( ( stdlib_base_ln( r ) - stdlib_base_ln( sigma ) ) - ( 0.5 * r * r ) );
 }

lib/node_modules/@stdlib/stats/base/dists/rayleigh/pdf/test/fixtures/julia/runner.jl

@@ -61,6 +61,40 @@ function gen( x, sigma, name )
 	close( outfile );
 end
 
+"""
+	gen_extended( x, sigma, name )
+
+Generate fixture data in extended precision and write to file.
+
+The PDF is evaluated using the textbook formula `x/sigma^2 * exp(-x^2/(2*sigma^2))` in 2048-bit `BigFloat`, as Distributions.jl evaluates it in `Float64`, where `sigma^2` and `x^2` overflow or underflow for scale parameters far from `1`.
+
+# Arguments
+
+* `x`: input value
+* `sigma`: scale parameter
+* `name::AbstractString`: output filename
+"""
+function gen_extended( x, sigma, name )
+	z = setprecision( BigFloat, 2048 ) do
+		X = BigFloat.( x );
+		S = BigFloat.( sigma );
+		Float64.( ( X ./ ( S .^ 2 ) ) .* exp.( .-( X .^ 2 ) ./ ( 2 .* ( S .^ 2 ) ) ) );
+	end
+
+	# Store data to be written to file as a collection:
+	data = Dict([
+		("x", x),
+		("sigma", sigma),
+		("expected", z)
+	]);
+
+	# Write the data to the output filepath as JSON:
+	outfile = open( joinpath( dir, name ), "w" );
+	write( outfile, JSON.json(data) );
+	write( outfile, "\n" );
+	close( outfile );
+end
+
 # Get the filename:
 file = @__FILE__;
 
@@ -81,3 +115,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 ) .* 3.0 ) .- 3.0 );
+gen_extended( x, sigma, "tiny_scale.json" );
+
+# Huge scale parameter (`sigma^2` overflows):
+sigma = exp10.( ( rand( 1000 ) .* 45.0 ) .+ 155.0 );
+x = sigma .* exp10.( ( rand( 1000 ) .* 3.0 ) .- 3.0 );
+gen_extended( x, sigma, "huge_scale.json" );

lib/node_modules/@stdlib/stats/base/dists/rayleigh/pdf/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' );
 
 
@@ -86,6 +88,20 @@ tape( 'if provided a valid `sigma`, the function returns a function which return
 	t.end();
 });
 
+tape( 'the created function returns `0` when `x/sigma` overflows', function test( t ) {
+	var pdf;
+	var y;
+
+	pdf = factory( 1.0e-10 );
+	y = pdf( 1.0e300 );
+	t.strictEqual( y, 0.0, 'returns expected value' );
+
+	pdf = factory( 1.0e-300 );
+	y = pdf( 1.0e308 );
+	t.strictEqual( y, 0.0, 'returns expected value' );
+	t.end();
+});
+
 tape( 'if provided a negative `sigma`, the created function always returns `NaN`', function test( t ) {
 	var pdf;
 	var y;
@@ -187,3 +203,45 @@ tape( 'the created function evaluates the pdf for `x` given large scale paramete
 	}
 	t.end();
 });
+
+tape( 'the created function evaluates the pdf for `x` given tiny scale parameter `sigma`', function test( t ) {
+	var expected;
+	var sigma;
+	var pdf;
+	var i;
+	var x;
+	var y;
+
+	expected = tinyScale.expected;
+	x = tinyScale.x;
+	sigma = tinyScale.sigma;
+	for ( i = 0; i < x.length; i++ ) {
+		pdf = factory( sigma[i] );
+		y = pdf( 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 pdf for `x` given huge scale parameter `sigma`', function test( t ) {
+	var expected;
+	var sigma;
+	var pdf;
+	var i;
+	var x;
+	var y;
+
+	expected = hugeScale.expected;
+	x = hugeScale.x;
+	sigma = hugeScale.sigma;
+	for ( i = 0; i < x.length; i++ ) {
+		pdf = factory( sigma[i] );
+		y = pdf( 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/pdf/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' );
 
 
@@ -66,6 +68,17 @@ tape( 'if provided `+infinity` for `x` and a valid `sigma`, the function returns
 	t.end();
 });
 
+tape( 'the function returns `0` when `x/sigma` overflows', opts, function test( t ) {
+	var y;
+
+	y = pdf( 1.0e300, 1.0e-10 );
+	t.strictEqual( y, 0.0, 'returns expected value' );
+
+	y = pdf( 1.0e308, 1.0e-300 );
+	t.strictEqual( y, 0.0, 'returns expected value' );
+	t.end();
+});
+
 tape( 'if provided `-infinity` for `x` and a valid `sigma`, the function returns `0`', opts, function test( t ) {
 	var y = pdf( NINF, 1.0 );
 	t.strictEqual( y, 0.0, 'returns expected value' );
@@ -164,3 +177,41 @@ tape( 'the function evaluates the pdf for `x` given large scale parameter `sigma
 	}
 	t.end();
 });
+
+tape( 'the function evaluates the pdf 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 = pdf( 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 pdf 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 = pdf( x[i], sigma[i] );
+		if ( expected[i] !== null ) {
+			t.strictEqual( isAlmostSameValue( y, expected[ i ], 2 ), true, 'returns expected value' );
+		}
+	}
+	t.end();
+});