fix: use log1p above the mode in stats/base/dists/triangular/logcdf

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

@@ -24,6 +24,7 @@ var constantFunction = require( '@stdlib/utils/constant-function' );
 var isnan = require( '@stdlib/math/base/assert/is-nan' );
 var pow = require( '@stdlib/math/base/special/pow' );
 var ln = require( '@stdlib/math/base/special/ln' );
+var log1p = require( '@stdlib/math/base/special/log1p' );
 var NINF = require( '@stdlib/constants/float64/ninf' );
 
 
@@ -88,7 +89,7 @@ function factory( a, b, c ) {
 		}
 		// Case: x > c
 		if ( x < b ) {
-			return ln( 1.0 - ( pow( b - x, 2.0 ) / denom2 ) );
+			return log1p( -pow( b - x, 2.0 ) / denom2 );
 		}
 		// Case: x >= b
 		return 0.0;

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

@@ -23,6 +23,7 @@
 var isnan = require( '@stdlib/math/base/assert/is-nan' );
 var pow = require( '@stdlib/math/base/special/pow' );
 var ln = require( '@stdlib/math/base/special/ln' );
+var log1p = require( '@stdlib/math/base/special/log1p' );
 var NINF = require( '@stdlib/constants/float64/ninf' );
 
 
@@ -99,7 +100,7 @@ function logcdf( x, a, b, c ) {
 	}
 	// Case: x > c
 	if ( x < b ) {
-		return ln( 1.0 - ( pow( b - x, 2.0 ) / denom2 ) );
+		return log1p( -pow( b - x, 2.0 ) / denom2 );
 	}
 	// Case: x >= b
 	return 0.0;

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

@@ -40,6 +40,7 @@
       "dependencies": [
         "@stdlib/math/base/napi/quaternary",
         "@stdlib/math/base/assert/is-nan",
+        "@stdlib/math/base/special/log1p",
         "@stdlib/constants/float64/ninf",
         "@stdlib/math/base/special/ln",
         "@stdlib/math/base/special/pow"
@@ -58,6 +59,7 @@
       "libpath": [],
       "dependencies": [
         "@stdlib/math/base/assert/is-nan",
+        "@stdlib/math/base/special/log1p",
         "@stdlib/constants/float64/eps",
         "@stdlib/constants/float64/ninf",
         "@stdlib/math/base/special/ln",
@@ -78,6 +80,7 @@
       "libpath": [],
       "dependencies": [
         "@stdlib/math/base/assert/is-nan",
+        "@stdlib/math/base/special/log1p",
         "@stdlib/constants/float64/eps",
         "@stdlib/constants/float64/ninf",
         "@stdlib/math/base/special/ln",

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

@@ -21,6 +21,7 @@
 #include "stdlib/constants/float64/ninf.h"
 #include "stdlib/math/base/special/ln.h"
 #include "stdlib/math/base/special/pow.h"
+#include "stdlib/math/base/special/log1p.h"
 
 /**
 * Evaluates the natural logarithm of the cumulative distribution function (CDF) for a triangular distribution with lower limit `a`, upper limit `b`, and mode `c` at a value `x`.
@@ -53,7 +54,7 @@ double stdlib_base_dists_triangular_logcdf( const double x, const double a, cons
 		return ( 2.0 * stdlib_base_ln( x - a ) - stdlib_base_ln( ( b - a ) * ( c - a ) ) );
 	}
 	if ( x < b ) {
-		return stdlib_base_ln( 1.0 - ( stdlib_base_pow( b - x, 2.0 ) / ( ( b - a ) * ( b - c ) ) ) );
+		return stdlib_base_log1p( -stdlib_base_pow( b - x, 2.0 ) / ( ( b - a ) * ( b - c ) ) );
 	}
 	return 0.0;
 }

lib/node_modules/@stdlib/stats/base/dists/triangular/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/triangular/logcdf/test/fixtures/julia/runner.jl

@@ -16,14 +16,15 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-import Distributions: logcdf, TriangularDist
 import JSON
 
 """
 	gen( x, a, b, c, name )
 
 Generate fixture data and write to file.
 
+The logarithm of the CDF is evaluated in extended precision from the textbook formula, `ln((x-a)^2 / ((b-a)(c-a)))` for `x <= c` and `ln(1 - (b-x)^2 / ((b-a)(b-c)))` for `x > c`, as forming `1 - (b-x)^2/...` in double precision loses the digits of a result close to `0`.
+
 # Arguments
 
 * `x`: input value
@@ -44,8 +45,18 @@ julia> gen( x, a, b, c, "data.json" );
 """
 function gen( x, a, b, c, name )
 	z = Array{Float64}( undef, length(x) );
-	for i in eachindex(x)
-		z[ i ] = logcdf( TriangularDist( a[i], b[i], c[i] ), x[i] );
+	setprecision( BigFloat, 2048 ) do
+		for i in eachindex(x)
+			xi = BigFloat( x[i] );
+			ai = BigFloat( a[i] );
+			bi = BigFloat( b[i] );
+			ci = BigFloat( c[i] );
+			if xi <= ci
+				z[ i ] = Float64( log( ( xi-ai )^2 / ( ( bi-ai ) * ( ci-ai ) ) ) );
+			else
+				z[ i ] = Float64( log( 1 - ( ( bi-xi )^2 / ( ( bi-ai ) * ( bi-ci ) ) ) ) );
+			end
+		end
 	end
 
 	# Store data to be written to file as a collection:
@@ -93,3 +104,10 @@ b = ( rand( 1000 ) .* 80.0 ) .+ a;
 c = a .+ ( b .- a ) .* rand();
 x = a .+ ( b .- a ) .* rand();
 gen( x, a, b, c, "large_range.json" );
+
+# `x` close to `b`:
+a = rand( 1000 ) .* 20.0;
+b = ( rand( 1000 ) .* 20.0 ) .+ a .+ 1.0;
+c = a .+ ( ( b .- a ) .* rand( 1000 ) );
+x = b .- ( ( b .- c ) .* exp10.( range( -8.0, stop = -1.0, length = 1000 ) ) );
+gen( x, a, b, c, "near_b.json" );

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

@@ -31,6 +31,7 @@ var factory = require( './../lib/factory.js' );
 // FIXTURES //
 
 var smallRange = require( './fixtures/julia/small_range.json' );
+var nearB = require( './fixtures/julia/near_b.json' );
 var mediumRange = require( './fixtures/julia/medium_range.json' );
 var largeRange = require( './fixtures/julia/large_range.json' );
 
@@ -171,7 +172,7 @@ tape( 'the created function evaluates the logcdf for `x` given a small range `b
 	for ( i = 0; i < x.length; i++ ) {
 		logcdf = factory( a[i], b[i], c[i] );
 		y = logcdf( x[i] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 6 ), true, 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 4 ), true, 'returns expected value' );
 	}
 	t.end();
 });
@@ -194,7 +195,7 @@ tape( 'the created function evaluates the logcdf for `x` given a medium range `b
 	for ( i = 0; i < x.length; i++ ) {
 		logcdf = factory( a[i], b[i], c[i] );
 		y = logcdf( x[i] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 0 ), true, 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 4 ), true, 'returns expected value' );
 	}
 	t.end();
 });
@@ -217,7 +218,30 @@ tape( 'the created function evaluates the logcdf for `x` given a large range `b
 	for ( i = 0; i < x.length; i++ ) {
 		logcdf = factory( a[i], b[i], c[i] );
 		y = logcdf( x[i] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 11 ), true, 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 4 ), true, 'returns expected value' );
+	}
+	t.end();
+});
+
+tape( 'the created function evaluates the logcdf for `x` close to `b`', function test( t ) {
+	var expected;
+	var logcdf;
+	var x;
+	var a;
+	var b;
+	var c;
+	var y;
+	var i;
+
+	expected = nearB.expected;
+	x = nearB.x;
+	a = nearB.a;
+	b = nearB.b;
+	c = nearB.c;
+	for ( i = 0; i < x.length; i++ ) {
+		logcdf = factory( a[i], b[i], c[i] );
+		y = logcdf( x[i] );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 2 ), true, 'returns expected value' );
 	}
 	t.end();
 });

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

@@ -31,6 +31,7 @@ var logcdf = require( './../lib' );
 // FIXTURES //
 
 var smallRange = require( './fixtures/julia/small_range.json' );
+var nearB = require( './fixtures/julia/near_b.json' );
 var mediumRange = require( './fixtures/julia/medium_range.json' );
 var largeRange = require( './fixtures/julia/large_range.json' );
 
@@ -105,7 +106,7 @@ tape( 'the function evaluates the logcdf for `x` given a small range `b - a`', f
 	c = smallRange.c;
 	for ( i = 0; i < x.length; i++ ) {
 		y = logcdf( x[i], a[i], b[i], c[i] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 6 ), true, 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 4 ), true, 'returns expected value' );
 	}
 	t.end();
 });
@@ -126,7 +127,7 @@ tape( 'the function evaluates the logcdf for `x` given a medium range `b - a`',
 	c = mediumRange.c;
 	for ( i = 0; i < x.length; i++ ) {
 		y = logcdf( x[i], a[i], b[i], c[i] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 0 ), true, 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 4 ), true, 'returns expected value' );
 	}
 	t.end();
 });
@@ -147,7 +148,28 @@ tape( 'the function evaluates the logcdf for `x` given a large range `b - a`', f
 	c = largeRange.c;
 	for ( i = 0; i < x.length; i++ ) {
 		y = logcdf( x[i], a[i], b[i], c[i] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 11 ), true, 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 4 ), true, 'returns expected value' );
+	}
+	t.end();
+});
+
+tape( 'the function evaluates the logcdf for `x` close to `b`', function test( t ) {
+	var expected;
+	var x;
+	var a;
+	var b;
+	var c;
+	var y;
+	var i;
+
+	expected = nearB.expected;
+	x = nearB.x;
+	a = nearB.a;
+	b = nearB.b;
+	c = nearB.c;
+	for ( i = 0; i < x.length; i++ ) {
+		y = logcdf( x[i], a[i], b[i], c[i] );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 2 ), true, 'returns expected value' );
 	}
 	t.end();
 });

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

@@ -32,6 +32,7 @@ var tryRequire = require( '@stdlib/utils/try-require' );
 // FIXTURES //
 
 var smallRange = require( './fixtures/julia/small_range.json' );
+var nearB = require( './fixtures/julia/near_b.json' );
 var mediumRange = require( './fixtures/julia/medium_range.json' );
 var largeRange = require( './fixtures/julia/large_range.json' );
 
@@ -114,7 +115,7 @@ tape( 'the function evaluates the logcdf for `x` given a small range `b - a`', o
 	c = smallRange.c;
 	for ( i = 0; i < x.length; i++ ) {
 		y = logcdf( x[i], a[i], b[i], c[i] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 6 ), true, 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 4 ), true, 'returns expected value' );
 	}
 	t.end();
 });
@@ -135,7 +136,7 @@ tape( 'the function evaluates the logcdf for `x` given a medium range `b - a`',
 	c = mediumRange.c;
 	for ( i = 0; i < x.length; i++ ) {
 		y = logcdf( x[i], a[i], b[i], c[i] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 0 ), true, 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 4 ), true, 'returns expected value' );
 	}
 	t.end();
 });
@@ -156,7 +157,28 @@ tape( 'the function evaluates the logcdf for `x` given a large range `b - a`', o
 	c = largeRange.c;
 	for ( i = 0; i < x.length; i++ ) {
 		y = logcdf( x[i], a[i], b[i], c[i] );
-		t.strictEqual( isAlmostSameValue( y, expected[ i ], 11 ), true, 'returns expected value' );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 4 ), true, 'returns expected value' );
+	}
+	t.end();
+});
+
+tape( 'the function evaluates the logcdf for `x` close to `b`', opts, function test( t ) {
+	var expected;
+	var x;
+	var a;
+	var b;
+	var c;
+	var y;
+	var i;
+
+	expected = nearB.expected;
+	x = nearB.x;
+	a = nearB.a;
+	b = nearB.b;
+	c = nearB.c;
+	for ( i = 0; i < x.length; i++ ) {
+		y = logcdf( x[i], a[i], b[i], c[i] );
+		t.strictEqual( isAlmostSameValue( y, expected[ i ], 2 ), true, 'returns expected value' );
 	}
 	t.end();
 });