Refactor to better handle NaN values

Author: kgryte

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

@@ -11,9 +11,6 @@
     accumulator function returns the current moving corrected sample standard
     deviation.
 
-    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 corrected sample standard deviation values will
     have less statistical support than subsequent corrected sample standard
     deviation values, as `W` values are needed to fill the window buffer. Until

lib/node_modules/@stdlib/stats/incr/mstdev/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' );
 var sqrt = require( '@stdlib/math/base/special/sqrt' );
 
 
@@ -104,6 +105,8 @@ function incrmstdev( W, mean ) {
 	* // returns ~5.29
 	*/
 	function accumulator1( x ) {
+		var k;
+		var v;
 		if ( arguments.length === 0 ) {
 			if ( N === 0 ) {
 				return null;
@@ -119,9 +122,13 @@ function incrmstdev( 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 ) ) {
+			M2 = NaN;
+		}
+		// Case: initial window...
+		else if ( N < W ) {
+			buf[ i ] = x; // update buffer
 			N += 1;
 			delta = x - mu;
 			mu += delta / N;
@@ -131,18 +138,41 @@ function incrmstdev( W, mean ) {
 			}
 			return sqrt( M2 / (N-1) );
 		}
-		// N = W = 1
+		// Case: N = W = 1
 		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...
+		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 ) ) {
+						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 sqrt( M2 / n );
 	}
 
@@ -154,6 +184,7 @@ function incrmstdev( W, mean ) {
 	* @returns {(number|null)} corrected sample standard deviation or null
 	*/
 	function accumulator2( x ) {
+		var k;
 		if ( arguments.length === 0 ) {
 			if ( N === 0 ) {
 				return null;
@@ -166,18 +197,40 @@ function incrmstdev( 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 ) ) {
+			M2 = NaN;
+		}
+		// Case: initial window...
+		else if ( N < W ) {
+			buf[ i ] = x; // update buffer
 			N += 1;
 			delta = x - mu;
 			M2 += delta * delta;
 			return sqrt( M2 / N );
 		}
-		// Update the existing window...
-		tmp = buf[ i ];
+		// Case: outgoing value is NaN, and, thus, we need to compute the accumulated values...
+		if ( isnan( buf[ i ] ) ) {
+			M2 = 0.0;
+			for ( k = 0; k < W; k++ ) {
+				if ( k !== i ) {
+					if ( isnan( buf[ k ] ) ) {
+						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 sqrt( M2 / W );
 	}
 }