fix: avoid cancellation near the edges of the support in stats/base/dists/cosine/logpdf

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

@@ -77,7 +77,9 @@ function factory( mu, s ) {
 			return NINF;
 		}
 		z = ( x - mu ) / s;
-		return ln( 1.0 + cospi( z ) ) - ln( 2.0 * s );
+
+		// `1 + cos(pi*z)` cancels near the edges of the support (`z = -1` and `z = 1`); `1 + cos(pi*z) = 2*cos(pi*z/2)^2` has no subtraction:
+		return ( 2.0 * ln( cospi( z / 2.0 ) ) ) - ln( s );
 	}
 }
 

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

@@ -90,7 +90,9 @@ function logpdf( x, mu, s ) {
 		return NINF;
 	}
 	z = ( x - mu ) / s;
-	return ln( 1.0 + cospi( z ) ) - ln( 2.0 * s );
+
+	// `1 + cos(pi*z)` cancels near the edges of the support (`z = -1` and `z = 1`); `1 + cos(pi*z) = 2*cos(pi*z/2)^2` has no subtraction:
+	return ( 2.0 * ln( cospi( z / 2.0 ) ) ) - ln( s );
 }
 
 

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

@@ -54,5 +54,7 @@ double stdlib_base_dists_cosine_logpdf( const double x, const double mu, const d
 		return STDLIB_CONSTANT_FLOAT64_NINF;
 	}
 	const double p = ( x - mu ) / s;
-	return stdlib_base_ln( 1.0 + stdlib_base_cospi( p ) ) - stdlib_base_ln( 2.0 * s );
+
+	// `1 + cos(pi*z)` cancels near the edges of the support (`z = -1` and `z = 1`); `1 + cos(pi*z) = 2*cos(pi*z/2)^2` has no subtraction:
+	return ( 2.0 * stdlib_base_ln( stdlib_base_cospi( p / 2.0 ) ) ) - stdlib_base_ln( s );
 }

lib/node_modules/@stdlib/stats/base/dists/cosine/logpdf/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/cosine/logpdf/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: logpdf, Cosine
 import JSON
 
 """
 	gen( x, mu, s, name )
 
 Generate fixture data and write to file.
 
+The log-density is evaluated in extended precision from the textbook formula `ln((1 + cos(pi*z)) / (2s))`, with `z = (x-mu)/s`.
+
 # Arguments
 
 * `x`: input value
@@ -41,9 +42,9 @@ julia> gen( x, mu, s, "data.json" );
 ```
 """
 function gen( x, mu, s, name )
-	z = Array{Float64}( undef, length(x) );
-	for i in eachindex(x)
-		z[ i ] = logpdf( Cosine( mu[i], s[i] ), x[i] );
+	z = setprecision( BigFloat, 2048 ) do
+		zz = ( BigFloat.( x ) .- BigFloat.( mu ) ) ./ BigFloat.( s );
+		Float64.( log.( ( 1 .+ cos.( BigFloat( pi ) .* zz ) ) ./ ( 2 .* BigFloat.( s ) ) ) )
 	end
 
 	# Store data to be written to file as a collection:
@@ -70,20 +71,27 @@ file = @__FILE__;
 # Extract the directory in which this file resides:
 dir = dirname( file );
 
-# Negative mean:
-x = rand( 1000 ) .* 10.0 .- 20.0;
+# Negative mean (`x` inside the support):
 mu = rand( 1000 ) .* -10.0;
-s = rand( 1000 ) .* 5.0;
+s = ( rand( 1000 ) .* 5.0 ) .+ 0.1;
+x = mu .+ ( s .* ( ( rand( 1000 ) .* 2.0 ) .- 1.0 ) );
 gen( x, mu, s, "negative_mean.json" );
 
-# Positive mean:
-x = rand( 1000 ) .* 10.0 .- 20.0;
+# Positive mean (`x` inside the support):
 mu = rand( 1000 ) .* 10.0;
-s = rand( 1000 ) .* 5.0;
+s = ( rand( 1000 ) .* 5.0 ) .+ 0.1;
+x = mu .+ ( s .* ( ( rand( 1000 ) .* 2.0 ) .- 1.0 ) );
 gen( x, mu, s, "positive_mean.json" );
 
-# Large variance:
-x = rand( 1000 ) .* 5.0;
+# Large variance (`x` inside the support):
 mu = rand( 1000 );
-s = rand( 1000 ) .* 20.0;
+s = ( rand( 1000 ) .* 20.0 ) .+ 1.0;
+x = mu .+ ( s .* ( ( rand( 1000 ) .* 2.0 ) .- 1.0 ) );
 gen( x, mu, s, "large_variance.json" );
+
+# Close to the edges of the support (`mu = 0` and `s` a power of two, so that `z = x/s` is exact and the test measures the error from `1 + cos(pi*z)`):
+s = repeat( [ 0.5, 1.0, 2.0, 4.0 ], 250 );
+d = exp10.( range( -15.0, stop = -1.0, length = 500 ) );
+x = s .* vcat( 1.0 .- d, -1.0 .+ d );
+mu = zeros( 1000 );
+gen( x, mu, s, "edges.json" );

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

@@ -33,6 +33,7 @@ var factory = require( './../lib/factory.js' );
 var positiveMean = require( './fixtures/julia/positive_mean.json' );
 var negativeMean = require( './fixtures/julia/negative_mean.json' );
 var largeVariance = require( './fixtures/julia/large_variance.json' );
+var edges = require( './fixtures/julia/edges.json' );
 
 
 // TESTS //
@@ -183,7 +184,7 @@ tape( 'the created function evaluates the logpdf for `x` given positive `mu`', f
 		logpdf = factory( mu[i], s[i] );
 		y = logpdf( x[i] );
 		if ( expected[i] !== null ) {
-			t.strictEqual( isAlmostSameValue( y, expected[i], 0 ), true, 'returns expected value' );
+			t.strictEqual( isAlmostSameValue( y, expected[i], 1062 ), true, 'returns expected value' );
 		}
 	}
 	t.end();
@@ -206,7 +207,7 @@ tape( 'the created function evaluates the logpdf for `x` given negative `mu`', f
 		logpdf = factory( mu[i], s[i] );
 		y = logpdf( x[i] );
 		if ( expected[i] !== null ) {
-			t.strictEqual( isAlmostSameValue( y, expected[i], 1 ), true, 'returns expected value' );
+			t.strictEqual( isAlmostSameValue( y, expected[i], 1530 ), true, 'returns expected value' );
 		}
 	}
 	t.end();
@@ -229,7 +230,30 @@ tape( 'the created function evaluates the logpdf for `x` given large variance (
 		logpdf = factory( mu[i], s[i] );
 		y = logpdf( x[i] );
 		if ( expected[i] !== null ) {
-			t.strictEqual( isAlmostSameValue( y, expected[i], 8 ), true, 'returns expected value' );
+			t.strictEqual( isAlmostSameValue( y, expected[i], 261 ), true, 'returns expected value' );
+		}
+	}
+	t.end();
+});
+
+tape( 'the created function evaluates the logpdf for `x` close to the edges of the support', function test( t ) {
+	var expected;
+	var logpdf;
+	var mu;
+	var s;
+	var x;
+	var y;
+	var i;
+
+	expected = edges.expected;
+	x = edges.x;
+	mu = edges.mu;
+	s = edges.s;
+	for ( i = 0; i < x.length; i++ ) {
+		logpdf = factory( mu[i], s[i] );
+		y = logpdf( x[i] );
+		if ( expected[i] !== null ) {
+			t.strictEqual( isAlmostSameValue( y, expected[i], 1 ), true, 'returns expected value' );
 		}
 	}
 	t.end();

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

@@ -33,6 +33,7 @@ var logpdf = require( './../lib' );
 var positiveMean = require( './fixtures/julia/positive_mean.json' );
 var negativeMean = require( './fixtures/julia/negative_mean.json' );
 var largeVariance = require( './fixtures/julia/large_variance.json' );
+var edges = require( './fixtures/julia/edges.json' );
 
 
 // TESTS //
@@ -139,7 +140,7 @@ tape( 'the function evaluates the logpdf for `x` given positive `mu`', function
 	for ( i = 0; i < x.length; i++ ) {
 		y = logpdf( x[i], mu[i], s[i] );
 		if ( expected[i] !== null ) {
-			t.strictEqual( isAlmostSameValue( y, expected[i], 0 ), true, 'returns expected value' );
+			t.strictEqual( isAlmostSameValue( y, expected[i], 1062 ), true, 'returns expected value' );
 		}
 	}
 	t.end();
@@ -160,7 +161,7 @@ tape( 'the function evaluates the logpdf for `x` given negative `mu`', function
 	for ( i = 0; i < x.length; i++ ) {
 		y = logpdf( x[i], mu[i], s[i] );
 		if ( expected[i] !== null ) {
-			t.strictEqual( isAlmostSameValue( y, expected[i], 1 ), true, 'returns expected value' );
+			t.strictEqual( isAlmostSameValue( y, expected[i], 1530 ), true, 'returns expected value' );
 		}
 	}
 	t.end();
@@ -181,7 +182,28 @@ tape( 'the function evaluates the logpdf for `x` given large variance ( = large
 	for ( i = 0; i < x.length; i++ ) {
 		y = logpdf( x[i], mu[i], s[i] );
 		if ( expected[i] !== null ) {
-			t.strictEqual( isAlmostSameValue( y, expected[i], 8 ), true, 'returns expected value' );
+			t.strictEqual( isAlmostSameValue( y, expected[i], 261 ), true, 'returns expected value' );
+		}
+	}
+	t.end();
+});
+
+tape( 'the function evaluates the logpdf for `x` close to the edges of the support', function test( t ) {
+	var expected;
+	var mu;
+	var x;
+	var s;
+	var y;
+	var i;
+
+	expected = edges.expected;
+	x = edges.x;
+	mu = edges.mu;
+	s = edges.s;
+	for ( i = 0; i < x.length; i++ ) {
+		y = logpdf( x[i], mu[i], s[i] );
+		if ( expected[i] !== null ) {
+			t.strictEqual( isAlmostSameValue( y, expected[i], 1 ), true, 'returns expected value' );
 		}
 	}
 	t.end();

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

@@ -34,6 +34,7 @@ var tryRequire = require( '@stdlib/utils/try-require' );
 var positiveMean = require( './fixtures/julia/positive_mean.json' );
 var negativeMean = require( './fixtures/julia/negative_mean.json' );
 var largeVariance = require( './fixtures/julia/large_variance.json' );
+var edges = require( './fixtures/julia/edges.json' );
 
 
 // VARIABLES //
@@ -148,7 +149,7 @@ tape( 'the function evaluates the logpdf for `x` given positive `mu`', opts, fun
 	for ( i = 0; i < x.length; i++ ) {
 		y = logpdf( x[i], mu[i], s[i] );
 		if ( expected[i] !== null ) {
-			t.strictEqual( isAlmostSameValue( y, expected[i], 0 ), true, 'returns expected value' );
+			t.strictEqual( isAlmostSameValue( y, expected[i], 1062 ), true, 'returns expected value' );
 		}
 	}
 	t.end();
@@ -169,7 +170,7 @@ tape( 'the function evaluates the logpdf for `x` given negative `mu`', opts, fun
 	for ( i = 0; i < x.length; i++ ) {
 		y = logpdf( x[i], mu[i], s[i] );
 		if ( expected[i] !== null ) {
-			t.strictEqual( isAlmostSameValue( y, expected[i], 1 ), true, 'returns expected value' );
+			t.strictEqual( isAlmostSameValue( y, expected[i], 1530 ), true, 'returns expected value' );
 		}
 	}
 	t.end();
@@ -190,7 +191,28 @@ tape( 'the function evaluates the logpdf for `x` given large variance ( = large
 	for ( i = 0; i < x.length; i++ ) {
 		y = logpdf( x[i], mu[i], s[i] );
 		if ( expected[i] !== null ) {
-			t.strictEqual( isAlmostSameValue( y, expected[i], 8 ), true, 'returns expected value' );
+			t.strictEqual( isAlmostSameValue( y, expected[i], 261 ), true, 'returns expected value' );
+		}
+	}
+	t.end();
+});
+
+tape( 'the function evaluates the logpdf for `x` close to the edges of the support', opts, function test( t ) {
+	var expected;
+	var mu;
+	var x;
+	var s;
+	var y;
+	var i;
+
+	expected = edges.expected;
+	x = edges.x;
+	mu = edges.mu;
+	s = edges.s;
+	for ( i = 0; i < x.length; i++ ) {
+		y = logpdf( x[i], mu[i], s[i] );
+		if ( expected[i] !== null ) {
+			t.strictEqual( isAlmostSameValue( y, expected[i], 1 ), true, 'returns expected value' );
 		}
 	}
 	t.end();