fix: avoid underflow close to mu in stats/base/dists/levy/logcdf

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

@@ -23,6 +23,7 @@
 var constantFunction = require( '@stdlib/utils/constant-function' );
 var isnan = require( '@stdlib/math/base/assert/is-nan' );
 var erfc = require( '@stdlib/math/base/special/erfc' );
+var erfcx = require( '@stdlib/math/base/special/erfcx' );
 var sqrt = require( '@stdlib/math/base/special/sqrt' );
 var ln = require( '@stdlib/math/base/special/ln' );
 var NINF = require( '@stdlib/constants/float64/ninf' );
@@ -68,15 +69,21 @@ function factory( mu, c ) {
 	* // returns <number>
 	*/
 	function logcdf( x ) {
+		var w;
 		var z;
 		if ( isnan( x ) ) {
 			return NaN;
 		}
 		if ( x < mu ) {
 			return NINF;
 		}
-		z = sqrt( c / ( 2.0 * ( x-mu ) ) );
-		return ln( erfc( z ) );
+		w = c / ( 2.0 * ( x-mu ) );
+		z = sqrt( w );
+		if ( z < 1.0 ) {
+			return ln( erfc( z ) );
+		}
+		// `erfc(z)` underflows for `z > ~26.5`, so use `ln(erfc(z)) = ln(erfcx(z)) - z^2`:
+		return ln( erfcx( z ) ) - w;
 	}
 }
 

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

@@ -22,6 +22,7 @@
 
 var isnan = require( '@stdlib/math/base/assert/is-nan' );
 var erfc = require( '@stdlib/math/base/special/erfc' );
+var erfcx = require( '@stdlib/math/base/special/erfcx' );
 var sqrt = require( '@stdlib/math/base/special/sqrt' );
 var NINF = require( '@stdlib/constants/float64/ninf' );
 var ln = require( '@stdlib/math/base/special/ln' );
@@ -67,6 +68,7 @@ var ln = require( '@stdlib/math/base/special/ln' );
 * // returns NaN
 */
 function logcdf( x, mu, c ) {
+	var w;
 	var z;
 	if (
 		isnan( x ) ||
@@ -79,8 +81,13 @@ function logcdf( x, mu, c ) {
 	if ( x < mu ) {
 		return NINF;
 	}
-	z = sqrt( c / ( 2.0 * ( x-mu ) ) );
-	return ln( erfc( z ) );
+	w = c / ( 2.0 * ( x-mu ) );
+	z = sqrt( w );
+	if ( z < 1.0 ) {
+		return ln( erfc( z ) );
+	}
+	// `erfc(z)` underflows for `z > ~26.5`, so use `ln(erfc(z)) = ln(erfcx(z)) - z^2`:
+	return ln( erfcx( z ) ) - w;
 }
 
 

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

@@ -41,6 +41,7 @@
         "@stdlib/math/base/napi/ternary",
         "@stdlib/math/base/assert/is-nan",
         "@stdlib/math/base/special/erfc",
+        "@stdlib/math/base/special/erfcx",
         "@stdlib/math/base/special/sqrt",
         "@stdlib/math/base/special/ln",
         "@stdlib/constants/float64/ninf"
@@ -60,6 +61,7 @@
       "dependencies": [
         "@stdlib/math/base/assert/is-nan",
         "@stdlib/math/base/special/erfc",
+        "@stdlib/math/base/special/erfcx",
         "@stdlib/math/base/special/sqrt",
         "@stdlib/math/base/special/ln",
         "@stdlib/constants/float64/ninf",
@@ -80,6 +82,7 @@
       "dependencies": [
         "@stdlib/math/base/assert/is-nan",
         "@stdlib/math/base/special/erfc",
+        "@stdlib/math/base/special/erfcx",
         "@stdlib/math/base/special/sqrt",
         "@stdlib/math/base/special/ln",
         "@stdlib/constants/float64/ninf",

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

@@ -19,6 +19,7 @@
 #include "stdlib/stats/base/dists/levy/logcdf.h"
 #include "stdlib/math/base/assert/is_nan.h"
 #include "stdlib/math/base/special/erfc.h"
+#include "stdlib/math/base/special/erfcx.h"
 #include "stdlib/math/base/special/sqrt.h"
 #include "stdlib/constants/float64/ninf.h"
 #include "stdlib/math/base/special/ln.h"
@@ -36,6 +37,7 @@
 * // returns ~-0.735
 */
 double stdlib_base_dists_levy_logcdf( const double x, const double mu, const double c ) {
+	double w;
 	double z;
 	if (
 		stdlib_base_is_nan( x ) ||
@@ -48,6 +50,11 @@ double stdlib_base_dists_levy_logcdf( const double x, const double mu, const dou
 	if ( x < mu ) {
 		return STDLIB_CONSTANT_FLOAT64_NINF;
 	}
-	z = stdlib_base_sqrt( c / ( 2.0 * ( x-mu ) ) );
-	return stdlib_base_ln( stdlib_base_erfc( z ) );
+	w = c / ( 2.0 * ( x-mu ) );
+	z = stdlib_base_sqrt( w );
+	if ( z < 1.0 ) {
+		return stdlib_base_ln( stdlib_base_erfc( z ) );
+	}
+	// `erfc(z)` underflows for `z > ~26.5`, so use `ln(erfc(z)) = ln(erfcx(z)) - z^2`:
+	return stdlib_base_ln( stdlib_base_erfcx( z ) ) - w;
 }

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

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

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

@@ -16,14 +16,16 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-import Distributions: logcdf, Levy
+import SpecialFunctions: erfc
 import JSON
 
 """
 	gen( x, mu, c, name )
 
 Generate fixture data and write to file.
 
+The logarithm of the CDF is evaluated in extended precision from the textbook formula `ln(erfc(sqrt(c / (2*(x-mu)))))`, as `erfc` underflows in double precision close to `mu`.
+
 # Arguments
 
 * `x`: input value
@@ -41,9 +43,8 @@ julia> gen( x, mu, c, "data.json" );
 ```
 """
 function gen( x, mu, c, name )
-	z = Array{Float64}( undef, length(x) );
-	for i in eachindex(x)
-		z[ i ] = logcdf( Levy( mu[i], c[i] ), x[i] );
+	z = setprecision( BigFloat, 2048 ) do
+		Float64.( log.( erfc.( sqrt.( BigFloat.( c ) ./ ( 2 .* ( BigFloat.( x ) .- BigFloat.( mu ) ) ) ) ) ) )
 	end
 
 	# Store data to be written to file as a collection:
@@ -87,3 +88,9 @@ mu = rand( 1000 );
 x = rand( 1000 ) .* 15.0 .+ mu .+ 0.1;
 c = ( rand( 1000 ) .* 15.0 ) .+ 1.0;
 gen( x, mu, c, "large_variance.json" );
+
+# Close to `mu`:
+mu = ( rand( 1000 ) .* 20.0 ) .- 10.0;
+x = mu .+ exp10.( range( -15.0, stop = 0.0, length = 1000 ) );
+c = ( rand( 1000 ) .* 5.0 ) .+ 0.1;
+gen( x, mu, c, "near_mu.json" );

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

@@ -31,6 +31,7 @@ var factory = require( './../lib/factory.js' );
 // FIXTURES //
 
 var positiveMean = require( './fixtures/julia/positive_mean.json' );
+var nearMu = require( './fixtures/julia/near_mu.json' );
 var negativeMean = require( './fixtures/julia/negative_mean.json' );
 var largeVariance = require( './fixtures/julia/large_variance.json' );
 
@@ -168,7 +169,7 @@ tape( 'the created function evaluates the logcdf for `x` given positive `mu`', f
 	for ( i = 0; i < x.length; i++ ) {
 		logcdf = factory( mu[i], c[i] );
 		y = logcdf( x[i] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 0 ), true, 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 3 ), true, 'returns expected value' );
 	}
 	t.end();
 });
@@ -189,7 +190,7 @@ tape( 'the created function evaluates the logcdf for `x` given negative `mu`', f
 	for ( i = 0; i < x.length; i++ ) {
 		logcdf = factory( mu[i], c[i] );
 		y = logcdf( x[i] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 0 ), true, 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 3 ), true, 'returns expected value' );
 	}
 	t.end();
 });
@@ -210,7 +211,49 @@ tape( 'the created function evaluates the logcdf for `x` given large variance (
 	for ( i = 0; i < x.length; i++ ) {
 		logcdf = factory( mu[i], c[i] );
 		y = logcdf( x[i] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 0 ), true, 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 3 ), true, 'returns expected value' );
 	}
 	t.end();
 });
+
+tape( 'the created function evaluates the logcdf for `x` close to `mu`', function test( t ) {
+	var expected;
+	var logcdf;
+	var mu;
+	var c;
+	var i;
+	var x;
+	var y;
+
+	expected = nearMu.expected;
+	x = nearMu.x;
+	mu = nearMu.mu;
+	c = nearMu.c;
+	for ( i = 0; i < x.length; i++ ) {
+		logcdf = factory( mu[i], c[i] );
+		y = logcdf( x[i] );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 2 ), true, 'returns expected value' );
+	}
+	t.end();
+});
+
+tape( 'the created function evaluates the logcdf for `x` very close to `mu`', function test( t ) {
+	var logcdf;
+	var y;
+
+	logcdf = factory( 0.0, 1.0 );
+	y = logcdf( 1.0e-20 );
+	t.strictEqual( isAlmostSameValue( y, -5.0e19, 0 ), true, 'returns expected value' );
+
+	y = logcdf( 1.0e-200 );
+	t.strictEqual( isAlmostSameValue( y, -5.0e199, 0 ), true, 'returns expected value' );
+
+	t.end();
+});
+
+tape( 'if provided `x` equal to `mu`, the created function returns `-infinity`', function test( t ) {
+	var logcdf = factory( 2.0, 1.0 );
+	var y = logcdf( 2.0 );
+	t.strictEqual( y, NINF, 'returns expected value' );
+	t.end();
+});

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

@@ -31,6 +31,7 @@ var logcdf = require( './../lib' );
 // FIXTURES //
 
 var positiveMean = require( './fixtures/julia/positive_mean.json' );
+var nearMu = require( './fixtures/julia/near_mu.json' );
 var negativeMean = require( './fixtures/julia/negative_mean.json' );
 var largeVariance = require( './fixtures/julia/large_variance.json' );
 
@@ -116,7 +117,7 @@ tape( 'the function evaluates the logcdf for `x` given positive `mu`', function
 	c = positiveMean.c;
 	for ( i = 0; i < x.length; i++ ) {
 		y = logcdf( x[i], mu[i], c[i] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 0 ), true, 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 3 ), true, 'returns expected value' );
 	}
 	t.end();
 });
@@ -135,7 +136,7 @@ tape( 'the function evaluates the logcdf for `x` given negative `mu`', function
 	c = negativeMean.c;
 	for ( i = 0; i < x.length; i++ ) {
 		y = logcdf( x[i], mu[i], c[i] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 0 ), true, 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 3 ), true, 'returns expected value' );
 	}
 	t.end();
 });
@@ -154,7 +155,44 @@ tape( 'the function evaluates the logcdf for `x` given large variance ( = large
 	c = largeVariance.c;
 	for ( i = 0; i < x.length; i++ ) {
 		y = logcdf( x[i], mu[i], c[i] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 0 ), true, 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 3 ), true, 'returns expected value' );
 	}
 	t.end();
 });
+
+tape( 'the function evaluates the logcdf for `x` close to `mu`', function test( t ) {
+	var expected;
+	var mu;
+	var c;
+	var i;
+	var x;
+	var y;
+
+	expected = nearMu.expected;
+	x = nearMu.x;
+	mu = nearMu.mu;
+	c = nearMu.c;
+	for ( i = 0; i < x.length; i++ ) {
+		y = logcdf( x[i], mu[i], c[i] );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 2 ), true, 'returns expected value' );
+	}
+	t.end();
+});
+
+tape( 'the function evaluates the logcdf for `x` very close to `mu`', function test( t ) {
+	var y;
+
+	y = logcdf( 1.0e-20, 0.0, 1.0 );
+	t.strictEqual( isAlmostSameValue( y, -5.0e19, 0 ), true, 'returns expected value' );
+
+	y = logcdf( 1.0e-200, 0.0, 1.0 );
+	t.strictEqual( isAlmostSameValue( y, -5.0e199, 0 ), true, 'returns expected value' );
+
+	t.end();
+});
+
+tape( 'if provided `x` equal to `mu`, the function returns `-infinity`', function test( t ) {
+	var y = logcdf( 2.0, 2.0, 1.0 );
+	t.strictEqual( y, NINF, 'returns expected value' );
+	t.end();
+});

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

@@ -40,6 +40,7 @@ var opts = {
 // FIXTURES //
 
 var positiveMean = require( './fixtures/julia/positive_mean.json' );
+var nearMu = require( './fixtures/julia/near_mu.json' );
 var negativeMean = require( './fixtures/julia/negative_mean.json' );
 var largeVariance = require( './fixtures/julia/large_variance.json' );
 
@@ -125,7 +126,7 @@ tape( 'the function evaluates the logcdf for `x` given positive `mu`', opts, fun
 	c = positiveMean.c;
 	for ( i = 0; i < x.length; i++ ) {
 		y = logcdf( x[i], mu[i], c[i] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 0 ), true, 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 3 ), true, 'returns expected value' );
 	}
 	t.end();
 });
@@ -144,7 +145,7 @@ tape( 'the function evaluates the logcdf for `x` given negative `mu`', opts, fun
 	c = negativeMean.c;
 	for ( i = 0; i < x.length; i++ ) {
 		y = logcdf( x[i], mu[i], c[i] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 0 ), true, 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 3 ), true, 'returns expected value' );
 	}
 	t.end();
 });
@@ -163,7 +164,44 @@ tape( 'the function evaluates the logcdf for `x` given large variance ( = large
 	c = largeVariance.c;
 	for ( i = 0; i < x.length; i++ ) {
 		y = logcdf( x[i], mu[i], c[i] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 0 ), true, 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 3 ), true, 'returns expected value' );
 	}
 	t.end();
 });
+
+tape( 'the function evaluates the logcdf for `x` close to `mu`', opts, function test( t ) {
+	var expected;
+	var mu;
+	var c;
+	var i;
+	var x;
+	var y;
+
+	expected = nearMu.expected;
+	x = nearMu.x;
+	mu = nearMu.mu;
+	c = nearMu.c;
+	for ( i = 0; i < x.length; i++ ) {
+		y = logcdf( x[i], mu[i], c[i] );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 2 ), true, 'returns expected value' );
+	}
+	t.end();
+});
+
+tape( 'the function evaluates the logcdf for `x` very close to `mu`', opts, function test( t ) {
+	var y;
+
+	y = logcdf( 1.0e-20, 0.0, 1.0 );
+	t.strictEqual( isAlmostSameValue( y, -5.0e19, 0 ), true, 'returns expected value' );
+
+	y = logcdf( 1.0e-200, 0.0, 1.0 );
+	t.strictEqual( isAlmostSameValue( y, -5.0e199, 0 ), true, 'returns expected value' );
+
+	t.end();
+});
+
+tape( 'if provided `x` equal to `mu`, the function returns `-infinity`', opts, function test( t ) {
+	var y = logcdf( 2.0, 2.0, 1.0 );
+	t.strictEqual( y, NINF, 'returns expected value' );
+	t.end();
+});