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

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

@@ -64,6 +64,7 @@ function factory( mu, s ) {
 	* // returns <number>
 	*/
 	function pdf( x ) {
+		var c;
 		var z;
 		if ( isnan( x ) ) {
 			return NaN;
@@ -75,7 +76,10 @@ function factory( mu, s ) {
 			return 0.0;
 		}
 		z = ( x - mu ) / s;
-		return ( 1.0 + cospi( z ) ) / ( 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:
+		c = cospi( z / 2.0 );
+		return ( c * c ) / s;
 	}
 }
 

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

@@ -69,6 +69,7 @@ var PINF = require( '@stdlib/constants/float64/pinf' );
 * // returns Infinity
 */
 function pdf( x, mu, s ) {
+	var c;
 	var z;
 	if (
 		isnan( x ) ||
@@ -88,7 +89,10 @@ function pdf( x, mu, s ) {
 		return 0.0;
 	}
 	z = ( x - mu ) / s;
-	return ( 1.0 + cospi( z ) ) / ( 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:
+	c = cospi( z / 2.0 );
+	return ( c * c ) / s;
 }
 
 

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

@@ -52,5 +52,8 @@ double stdlib_base_dists_cosine_pdf( const double x, const double mu, const doub
 		return 0.0;
 	}
 	const double p = ( x - mu ) / s;
-	return ( 1.0 + stdlib_base_cospi( p ) ) / ( 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:
+	const double c = stdlib_base_cospi( p / 2.0 );
+	return ( c * c ) / s;
 }

lib/node_modules/@stdlib/stats/base/dists/cosine/pdf/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/pdf/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: pdf, Cosine
 import JSON
 
 """
 	gen( x, mu, s, name )
 
 Generate fixture data and write to file.
 
+The density is evaluated in extended precision from the textbook formula `(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 ] = pdf( Cosine( mu[i], s[i] ), x[i] );
+	z = setprecision( BigFloat, 2048 ) do
+		zz = ( BigFloat.( x ) .- BigFloat.( mu ) ) ./ BigFloat.( s );
+		Float64.( ( 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/pdf/test/test.factory.js

@@ -34,6 +34,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 //
@@ -184,7 +185,7 @@ tape( 'the created function evaluates the pdf for `x` given positive `mu`', func
 		pdf = factory( mu[i], s[i] );
 		y = pdf( x[i] );
 		if ( expected[i] !== null && !( expected[i] === 0.0 && y < EPS ) ) {
-			t.strictEqual( isAlmostSameValue( y, expected[ i ], 0 ), true, 'returns expected value' );
+			t.strictEqual( isAlmostSameValue( y, expected[ i ], 255 ), true, 'returns expected value' );
 		}
 	}
 	t.end();
@@ -207,7 +208,7 @@ tape( 'the created function evaluates the pdf for `x` given negative `mu`', func
 		pdf = factory( mu[i], s[i] );
 		y = pdf( x[i] );
 		if ( expected[i] !== null && !( expected[i] === 0.0 && y < EPS ) ) {
-			t.strictEqual( isAlmostSameValue( y, expected[ i ], 4 ), true, 'returns expected value' );
+			t.strictEqual( isAlmostSameValue( y, expected[ i ], 1742 ), true, 'returns expected value' );
 		}
 	}
 	t.end();
@@ -230,7 +231,30 @@ tape( 'the created function evaluates the pdf for `x` given large variance ( = l
 		pdf = factory( mu[i], s[i] );
 		y = pdf( x[i] );
 		if ( expected[i] !== null ) {
-			t.strictEqual( isAlmostSameValue( y, expected[ i ], 2 ), true, 'returns expected value' );
+			t.strictEqual( isAlmostSameValue( y, expected[ i ], 370 ), true, 'returns expected value' );
+		}
+	}
+	t.end();
+});
+
+tape( 'the created function evaluates the pdf for `x` close to the edges of the support', function test( t ) {
+	var expected;
+	var pdf;
+	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++ ) {
+		pdf = factory( mu[i], s[i] );
+		y = pdf( x[i] );
+		if ( expected[i] !== null && !( expected[i] === 0.0 && y < EPS ) ) {
+			t.strictEqual( isAlmostSameValue( y, expected[ i ], 3 ), true, 'returns expected value' );
 		}
 	}
 	t.end();

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

@@ -35,6 +35,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 //
@@ -149,7 +150,7 @@ tape( 'the function evaluates the pdf for `x` given positive `mu`', opts, functi
 	for ( i = 0; i < x.length; i++ ) {
 		y = pdf( x[i], mu[i], s[i] );
 		if ( expected[i] !== null && !( expected[i] === 0.0 && y < EPS ) ) {
-			t.strictEqual( isAlmostSameValue( y, expected[ i ], 0 ), true, 'returns expected value' );
+			t.strictEqual( isAlmostSameValue( y, expected[ i ], 255 ), true, 'returns expected value' );
 		}
 	}
 	t.end();
@@ -170,7 +171,7 @@ tape( 'the function evaluates the pdf for `x` given negative `mu`', opts, functi
 	for ( i = 0; i < x.length; i++ ) {
 		y = pdf( x[i], mu[i], s[i] );
 		if ( expected[i] !== null && !( expected[i] === 0.0 && y < EPS ) ) {
-			t.strictEqual( isAlmostSameValue( y, expected[ i ], 4 ), true, 'returns expected value' );
+			t.strictEqual( isAlmostSameValue( y, expected[ i ], 1742 ), true, 'returns expected value' );
 		}
 	}
 	t.end();
@@ -191,7 +192,28 @@ tape( 'the function evaluates the pdf for `x` given large variance ( = large `s`
 	for ( i = 0; i < x.length; i++ ) {
 		y = pdf( x[i], mu[i], s[i] );
 		if ( expected[i] !== null ) {
-			t.strictEqual( isAlmostSameValue( y, expected[ i ], 2 ), true, 'returns expected value' );
+			t.strictEqual( isAlmostSameValue( y, expected[ i ], 370 ), true, 'returns expected value' );
+		}
+	}
+	t.end();
+});
+
+tape( 'the function evaluates the pdf 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 = pdf( x[i], mu[i], s[i] );
+		if ( expected[i] !== null && !( expected[i] === 0.0 && y < EPS ) ) {
+			t.strictEqual( isAlmostSameValue( y, expected[ i ], 3 ), true, 'returns expected value' );
 		}
 	}
 	t.end();

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

@@ -34,6 +34,7 @@ var pdf = 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 //
@@ -140,7 +141,7 @@ tape( 'the function evaluates the pdf for `x` given positive `mu`', function tes
 	for ( i = 0; i < x.length; i++ ) {
 		y = pdf( x[i], mu[i], s[i] );
 		if ( expected[i] !== null && !( expected[i] === 0.0 && y < EPS ) ) {
-			t.strictEqual( isAlmostSameValue( y, expected[ i ], 0 ), true, 'returns expected value' );
+			t.strictEqual( isAlmostSameValue( y, expected[ i ], 255 ), true, 'returns expected value' );
 		}
 	}
 	t.end();
@@ -161,7 +162,7 @@ tape( 'the function evaluates the pdf for `x` given negative `mu`', function tes
 	for ( i = 0; i < x.length; i++ ) {
 		y = pdf( x[i], mu[i], s[i] );
 		if ( expected[i] !== null && !( expected[i] === 0.0 && y < EPS ) ) {
-			t.strictEqual( isAlmostSameValue( y, expected[ i ], 4 ), true, 'returns expected value' );
+			t.strictEqual( isAlmostSameValue( y, expected[ i ], 1742 ), true, 'returns expected value' );
 		}
 	}
 	t.end();
@@ -182,7 +183,28 @@ tape( 'the function evaluates the pdf for `x` given large variance ( = large `s`
 	for ( i = 0; i < x.length; i++ ) {
 		y = pdf( x[i], mu[i], s[i] );
 		if ( expected[i] !== null ) {
-			t.strictEqual( isAlmostSameValue( y, expected[ i ], 2 ), true, 'returns expected value' );
+			t.strictEqual( isAlmostSameValue( y, expected[ i ], 370 ), true, 'returns expected value' );
+		}
+	}
+	t.end();
+});
+
+tape( 'the function evaluates the pdf 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 = pdf( x[i], mu[i], s[i] );
+		if ( expected[i] !== null && !( expected[i] === 0.0 && y < EPS ) ) {
+			t.strictEqual( isAlmostSameValue( y, expected[ i ], 3 ), true, 'returns expected value' );
 		}
 	}
 	t.end();