function incrmmeanvar( out, window ) {
	var meanvar;
	var delta;
	var buf;
	var tmp;
	var M2;
	var mu;
	var d1;
	var d2;
	var W;
	var N;
	var n;
	var i;
	if ( arguments.length === 1 ) {
		meanvar = [ 0.0, 0.0 ];
		W = out;
	} else {
		if ( !isArrayLike( out ) ) {
			throw new TypeError( format( 'invalid argument. Output argument must be an array-like object. Value: `%s`.', out ) );
		}
		meanvar = out;
		W = window;
	}
	if ( !isPositiveInteger( W ) ) {
		throw new TypeError( format( 'invalid argument. Window size must be a positive integer. Value: `%s`.', W ) );
	}
	buf = new Float64Array( W );
	n = W - 1;
	M2 = 0.0;
	mu = 0.0;
	i = -1;
	N = 0;

	return accumulator;

	/**
	* If provided a value, the accumulator function returns updated accumulated values. If not provided a value, the accumulator function returns the current accumulated values.
	*
	* @private
	* @param {number} [x] - input value
	* @returns {(ArrayLikeObject|null)} output array or null
	*/
	function accumulator( x ) {
		var k;
		var v;
		if ( arguments.length === 0 ) {
			if ( N === 0 ) {
				return null;
			}
			meanvar[ 0 ] = mu;
			if ( N === 1 ) {
				if ( isnan( M2 ) ) {
					meanvar[ 1 ] = NaN;
				} else {
					meanvar[ 1 ] = 0.0;
				}
			} else if ( N < W ) {
				meanvar[ 1 ] = M2 / (N-1);
			} else {
				meanvar[ 1 ] = M2 / n;
			}
			return meanvar;
		}
		// Update the index for managing the circular buffer:
		i = (i+1) % W;

		// Case: incoming value is NaN, the sliding second moment is automatically NaN...
		if ( isnan( x ) ) {
			N = W; // explicitly set to avoid `N < W` branch
			mu = NaN;
			M2 = NaN;
		}
		// Case: initial window...
		else if ( N < W ) {
			buf[ i ] = x; // update buffer
			N += 1;
			delta = x - mu;
			mu += delta / N;
			M2 += delta * (x - mu);

			meanvar[ 0 ] = mu;
			if ( N === 1 ) {
				meanvar[ 1 ] = 0.0;
			} else {
				meanvar[ 1 ] = M2 / (N-1);
			}
			return meanvar;
		}
		// Case: N = W = 1
		else if ( N === 1 ) {
			mu = x;
			M2 = 0.0;
			meanvar[ 0 ] = x;
			meanvar[ 1 ] = 0.0;
			return meanvar;
		}
		// 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
						mu = NaN;
						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;

		meanvar[ 0 ] = mu;
		meanvar[ 1 ] = M2 / n;
		return meanvar;
	}
}
Returns an accumulator function which incrementally computes a moving arithmetic mean and unbiased sample variance.

## Method

-   Let \\(W\\) be a window of \\(N\\) elements over which we want to compute an unbiased sample variance.

-   The difference between the unbiased sample variance in a window \\(W_i\\) and the unbiased sample variance in a window \\(W_{i+1})\\) is given by

    ```tex
    \Delta s^2 = s_{i+1}^2 - s_{i}^2
    ```

-   If we multiply both sides by \\(N-1\\),

    ```tex
    (N-1)(\Delta s^2) = (N-1)s_{i+1}^2 - (N-1)s_{i}^2
    ```

-   If we substitute the definition of the unbiased sample variance having the form

    ```tex
    \begin{align*}
    s^2 &= \frac{1}{N-1} \biggl( \sum_{i=1}^{N} (x_i - \bar{x})^2 \biggr) \\
        &= \frac{1}{N-1} \biggl( \sum_{i=1}^{N} (x_i^2 - 2\bar{x}x_i + \bar{x}^2) \biggr) \\
        &= \frac{1}{N-1} \biggl( \sum_{i=1}^{N} x_i^2 - 2\bar{x} \sum_{i=1}^{N} x_i + \sum_{i=1}^{N} \bar{x}^2) \biggr) \\
        &= \frac{1}{N-1} \biggl( \sum_{i=1}^{N} x_i^2 - \frac{2N\bar{x}\sum_{i=1}^{N} x_i}{N} + N\bar{x}^2 \biggr) \\
        &= \frac{1}{N-1} \biggl( \sum_{i=1}^{N} x_i^2 - 2N\bar{x}^2 + N\bar{x}^2 \biggr) \\
        &= \frac{1}{N-1} \biggl( \sum_{i=1}^{N} x_i^2 - N\bar{x}^2 \biggr)
    \end{align*}
    ```

    we return

    ```tex
    (N-1)(\Delta s^2) = \biggl(\sum_{k=1}^N x_k^2 - N\bar{x}_{i+1}^2 \biggr) - \biggl(\sum_{k=0}^{N-1} x_k^2 - N\bar{x}_{i}^2 \biggr)
    ```

-   This can be further simplified by recognizing that subtracting the sums reduces to \\(x_N^2 - x_0^2\\); in which case,

    ```tex
    \begin{align*}
    (N-1)(\Delta s^2) &= x_N^2 - x_0^2 - N\bar{x}_{i+1}^2 + N\bar{x}_{i}^2 \\
    &= x_N^2 - x_0^2 - N(\bar{x}_{i+1}^2 - \bar{x}_{i}^2) \\
    &= x_N^2 - x_0^2 - N(\bar{x}_{i+1} - \bar{x}_{i})(\bar{x}_{i+1} + \bar{x}_{i})
    \end{align*}
    ```

-   Recognizing that the difference of means can be expressed

    ```tex
    \bar{x}_{i+1} - \bar{x}_i = \frac{1}{N} \biggl( \sum_{k=1}^N x_k - \sum_{k=0}^{N-1} x_k \biggr) = \frac{x_N - x_0}{N}
    ```

    and substituting into the equation above

    ```tex
    (N-1)(\Delta s^2) = x_N^2 - x_0^2 - (x_N - x_0)(\bar{x}_{i+1} + \bar{x}_{i})
    ```

-   Rearranging terms gives us the update equation

    ```tex
    \begin{align*}
    (N-1)(\Delta s^2) &= (x_N - x_0)(x_N + x_0) - (x_N - x_0)(\bar{x}_{i+1} + \bar{x}_{i})
    &= (x_N - x_0)(x_N + x_0 - \bar{x}_{i+1} - \bar{x}_{i}) \\
    &= (x_N - x_0)(x_N - \bar{x}_{i+1} + x_0 - \bar{x}_{i})
    \end{align*}
    ```

@param {Collection} [out] - output array
@param {PositiveInteger} window - window size
@throws {TypeError} output argument must be array-like
@throws {TypeError} window size must be a positive integer
@returns {Function} accumulator function

@example
var accumulator = incrmmeanvar( 3 );

var v = accumulator();
// returns null

v = accumulator( 2.0 );
// returns [ 2.0, 0.0 ]

v = accumulator( -5.0 );
// returns [ -1.5, 24.5 ]

v = accumulator( 3.0 );
// returns [ 0.0, 19.0 ]

v = accumulator( 5.0 );
// returns [ 1.0, 28.0 ]

v = accumulator();
// returns [ 1.0, 28.0 ]