Refactor to better handle NaN values

Author: kgryte

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

@@ -10,9 +10,6 @@
     unbiased sample variance. If not provided a value, the accumulator function
     returns the current moving unbiased sample variance.
 
-    If provided `NaN` or a value which, when used in computations, results in
-    `NaN`, the accumulated value is `NaN` for all future invocations.
-
     The first `W-1` returned unbiased sample variance values will have less
     statistical support than subsequent unbiased sample variance values, as `W`
     values are needed to fill the window buffer. Until the window is full, the

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

@@ -22,6 +22,7 @@
 
 var isPositiveInteger = require( '@stdlib/assert/is-positive-integer' ).isPrimitive;
 var isNumber = require( '@stdlib/assert/is-number' ).isPrimitive;
+var isnan = require( '@stdlib/math/base/assert/is-nan' );
 
 
 // MAIN //
@@ -164,6 +165,8 @@ function incrmvariance( W, mean ) {
 	* // returns 28.0
 	*/
 	function accumulator1( x ) {
+		var k;
+		var v;
 		if ( arguments.length === 0 ) {
 			if ( N === 0 ) {
 				return null;
@@ -179,9 +182,14 @@ function incrmvariance( W, mean ) {
 		// Update the index for managing the circular buffer:
 		i = (i+1) % W;
 
-		// Determine if we should update the initial window...
-		if ( N < W ) {
-			buf[ i ] = x;
+		// Case: incoming value is NaN, the sliding second moment is automatically NaN...
+		if ( isnan( x ) ) {
+			N = W; // explicitly set to avoid `N < W` branch
+			M2 = NaN;
+		}
+		// Case: initial window...
+		else if ( N < W ) {
+			buf[ i ] = x; // update buffer
 			N += 1;
 			delta = x - mu;
 			mu += delta / N;
@@ -191,18 +199,42 @@ function incrmvariance( W, mean ) {
 			}
 			return M2 / (N-1);
 		}
-		// N = W = 1
-		if ( N === 1 ) {
+		// Case: N = W = 1
+		else if ( N === 1 ) {
 			return 0.0;
 		}
-		// Update the existing window...
-		tmp = buf[ i ];
+		// Case: outgoing value is NaN, and, thus, we need to compute the accumulated values...
+		else if ( isnan( buf[ i ] ) ) {
+			N = 1;
+			mu = x;
+			M2 = 0.0;
+			for ( k = 0; k < W; k++ ) {
+				if ( k !== i ) {
+					v = buf[ k ];
+					if ( isnan( v ) ) {
+						N = W; // explicitly set to avoid `N < W` branch
+						M2 = NaN;
+						break; // second moment is automatically NaN, so no need to continue
+					}
+					N += 1;
+					delta = v - mu;
+					mu += delta / N;
+					M2 += delta * (v - mu);
+				}
+			}
+		}
+		// Case: neither the current second moment nor the incoming value are NaN, so we need to update the accumulated values...
+		else if ( isnan( M2 ) === false ) {
+			tmp = buf[ i ];
+			delta = x - tmp;
+			d1 = tmp - mu;
+			mu += delta / W;
+			d2 = x - mu;
+			M2 += delta * (d1 + d2);
+		}
+		// Case: the current second moment is NaN, so nothing to do until the buffer no longer contains NaN values...
+
 		buf[ i ] = x;
-		delta = x - tmp;
-		d1 = tmp - mu;
-		mu += delta / W;
-		d2 = x - mu;
-		M2 += delta * (d1 + d2);
 		return M2 / n;
 	}
 
@@ -214,6 +246,7 @@ function incrmvariance( W, mean ) {
 	* @returns {(number|null)} unbiased sample variance or null
 	*/
 	function accumulator2( x ) {
+		var k;
 		if ( arguments.length === 0 ) {
 			if ( N === 0 ) {
 				return null;
@@ -226,18 +259,42 @@ function incrmvariance( W, mean ) {
 		// Update the index for managing the circular buffer:
 		i = (i+1) % W;
 
-		// Determine if we should update the initial window...
-		if ( N < W ) {
-			buf[ i ] = x;
+		// Case: incoming value is NaN, the sliding second moment is automatically NaN...
+		if ( isnan( x ) ) {
+			N = W; // explicitly set to avoid `N < W` branch
+			M2 = NaN;
+		}
+		// Case: initial window...
+		else if ( N < W ) {
+			buf[ i ] = x; // update buffer
 			N += 1;
 			delta = x - mu;
 			M2 += delta * delta;
 			return M2 / N;
 		}
-		// Update the existing window...
-		tmp = buf[ i ];
+		// Case: outgoing value is NaN, and, thus, we need to compute the accumulated values...
+		else if ( isnan( buf[ i ] ) ) {
+			M2 = 0.0;
+			for ( k = 0; k < W; k++ ) {
+				if ( k !== i ) {
+					if ( isnan( buf[ k ] ) ) {
+						N = W; // explicitly set to avoid `N < W` branch
+						M2 = NaN;
+						break; // second moment is automatically NaN, so no need to continue
+					}
+					delta = buf[ k ] - mu;
+					M2 += delta * delta;
+				}
+			}
+		}
+		// Case: neither the current second moment nor the incoming value are NaN, so we need to update the accumulated values...
+		else if ( isnan( M2 ) === false ) {
+			tmp = buf[ i ];
+			M2 += ( x-tmp ) * ( x-mu + tmp-mu );
+		}
+		// Case: the current second moment is NaN, so nothing to do until the buffer no longer contains NaN values...
+
 		buf[ i ] = x;
-		M2 += ( x-tmp ) * ( x-mu + tmp-mu );
 		return M2 / W;
 	}
 }