fix: clamp moving Pearson correlation distance to [0,2]

lib/node_modules/@stdlib/stats/incr/mpcorrdist/README.md

@@ -108,7 +108,7 @@ r = accumulator();
 
 -   Input values are **not** type checked. If provided `NaN` or a value which, when used in computations, results in `NaN`, the accumulated value is `NaN` for **at least** `W-1` future invocations. If non-numeric inputs are possible, you are advised to type check and handle accordingly **before** passing the value to the accumulator function.
 -   As `W` (x,y) pairs are needed to fill the window buffer, the first `W-1` returned values are calculated from smaller sample sizes. Until the window is full, each returned value is calculated from all provided values.
--   Due to limitations inherent in representing numeric values using floating-point format (i.e., the inability to represent numeric values with infinite precision), the [sample correlation distance][pearson-correlation] between perfectly correlated random variables may **not** be `0` or `2`. In fact, the [sample correlation distance][pearson-correlation] is **not** guaranteed to be strictly on the interval `[0,2]`. Any computed distance should, however, be within floating-point roundoff error.
+-   The computed [sample correlation distance][pearson-correlation] is clamped to the interval `[0,2]`. Due to floating-point rounding in the underlying Welford accumulation, the raw distance may deviate from this interval by a ULP or two; the clamp ensures the returned value always satisfies the mathematical constraint.
 
 </section>
 

lib/node_modules/@stdlib/stats/incr/mpcorrdist/docs/repl.txt

@@ -6,12 +6,10 @@
     The correlation distance is defined as one minus the Pearson product-moment
     correlation coefficient and, thus, resides on the interval [0,2].
 
-    However, due to limitations inherent in representing numeric values using
-    floating-point format (i.e., the inability to represent numeric values with
-    infinite precision), the correlation distance between perfectly correlated
-    random variables may *not* be `0` or `2`. In fact, the correlation distance
-    is *not* guaranteed to be strictly on the interval [0,2]. Any computed
-    distance should, however, be within floating-point roundoff error.
+    The computed correlation distance is clamped to the interval [0,2]. Due to
+    floating-point rounding in the underlying Welford accumulation, the raw
+    distance may deviate from this interval by a ULP or two; the clamp ensures
+    the returned value always satisfies the mathematical constraint.
 
     The `W` parameter defines the number of values over which to compute the
     moving sample correlation distance.

lib/node_modules/@stdlib/stats/incr/mpcorrdist/lib/main.js

@@ -100,15 +100,30 @@ function incrmpcorrdist( W, meanx, meany ) {
 	* @returns {(number|null)} sample correlation distance or null
 	*/
 	function accumulator( x, y ) {
+		var d;
 		var r;
 		if ( arguments.length === 0 ) {
 			r = pcorr();
 			if ( r === null ) {
 				return r;
 			}
-			return 1.0 - r;
+			d = 1.0 - r;
+			if ( d < 0.0 ) {
+				return 0.0;
+			}
+			if ( d > 2.0 ) {
+				return 2.0;
+			}
+			return d;
+		}
+		d = 1.0 - pcorr( x, y );
+		if ( d < 0.0 ) {
+			return 0.0;
+		}
+		if ( d > 2.0 ) {
+			return 2.0;
 		}
-		return 1.0 - pcorr( x, y );
+		return d;
 	}
 }
 

lib/node_modules/@stdlib/stats/incr/mpcorrdist/test/test.js

@@ -376,14 +376,14 @@ tape( 'the accumulator function computes a moving sample Pearson product-moment
 			if ( actual === expected ) {
 				t.strictEqual( actual, expected, 'returns expected value. dataset: '+i+'. window: '+j+'.' );
 			} else {
-				if ( actual < 0.0 ) {
-					actual = 0.0; // NOTE: this addresses occasional negative values due to accumulated floating-point error. Based on observation, typically `|actual| ≅ |expected|`, but `actual < 0` and `expected > 0`, suggesting that a sign got "flipped" along the way due to, e.g., operations which theoretically should compute to the same value, but do not due to floating-point error.
-				}
 				delta = abs( actual - expected );
 				if ( expected === 0.0 || actual === 0.0 ) {
-					tol = 10.0 * EPS;
+					tol = 1.0e2 * EPS;
 				} else {
 					tol = 1.0e6 * EPS * abs( expected );
+					if ( tol < 1.0e2 * EPS ) {
+						tol = 1.0e2 * EPS;
+					}
 				}
 				t.strictEqual( delta <= tol, true, 'dataset: '+i+'. window: '+j+'. expected: '+expected+'. actual: '+actual+'. tol: '+tol+'. delta: '+delta+'.' );
 			}
@@ -435,6 +435,9 @@ tape( 'the accumulator function computes a moving sample Pearson product-moment
 			} else {
 				delta = abs( actual - expected );
 				tol = 1.0e6 * EPS * abs( expected );
+				if ( tol < 1.0e2 * EPS ) {
+					tol = 1.0e2 * EPS;
+				}
 				t.strictEqual( delta <= tol, true, 'dataset: '+i+'. window: '+j+'. expected: '+expected+'. actual: '+actual+'. tol: '+tol+'. delta: '+delta+'.' );
 			}
 		}
@@ -817,3 +820,30 @@ tape( 'if provided `NaN`, the accumulated value is `NaN` for at least `W` invoca
 	}
 	t.end();
 });
+
+tape( 'the accumulator function clamps the correlation distance to [0,2] when floating-point rounding causes the coefficient to fall outside [-1,1]', function test( t ) {
+	var acc;
+	var v;
+	var i;
+
+	// Welford accumulation with nearly anti-correlated data at scale 10 yields r < -1 at the 8th datum, making `1 - r` exceed 2.0 before clamping.
+	acc = incrmpcorrdist( 2 );
+	for ( i = 0; i < 8; i++ ) {
+		v = acc( i * 10, ( -i * 10 ) + 1.0e-8 );
+		if ( v !== null ) {
+			t.strictEqual( v <= 2.0, true, 'returns value at most 2 at i='+i );
+		}
+	}
+	t.strictEqual( acc(), 2.0, 'no-argument query returns 2.0 after clamping' );
+
+	// Welford accumulation with nearly positively correlated data at scale 10 yields r > 1, making `1 - r` fall below 0.0 before clamping.
+	acc = incrmpcorrdist( 2 );
+	for ( i = 0; i < 14; i++ ) {
+		v = acc( i * 10, ( i * 10 ) + 1.0e-8 );
+		if ( v !== null ) {
+			t.strictEqual( v >= 0.0, true, 'returns value at least 0 at i='+i );
+		}
+	}
+	t.strictEqual( acc(), 0.0, 'no-argument query returns 0.0 after clamping' );
+	t.end();
+});