main.js (8548B)
1 /** 2 * @license Apache-2.0 3 * 4 * Copyright (c) 2018 The Stdlib Authors. 5 * 6 * Licensed under the Apache License, Version 2.0 (the "License"); 7 * you may not use this file except in compliance with the License. 8 * You may obtain a copy of the License at 9 * 10 * http://www.apache.org/licenses/LICENSE-2.0 11 * 12 * Unless required by applicable law or agreed to in writing, software 13 * distributed under the License is distributed on an "AS IS" BASIS, 14 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 15 * See the License for the specific language governing permissions and 16 * limitations under the License. 17 */ 18 19 'use strict'; 20 21 // MODULES // 22 23 var isPositiveInteger = require( '@stdlib/assert/is-positive-integer' ).isPrimitive; 24 var isNumber = require( '@stdlib/assert/is-number' ).isPrimitive; 25 var isnan = require( '@stdlib/math/base/assert/is-nan' ); 26 27 28 // MAIN // 29 30 /** 31 * Returns an accumulator function which incrementally computes a moving unbiased sample variance. 32 * 33 * ## Method 34 * 35 * - Let \\(W\\) be a window of \\(N\\) elements over which we want to compute an unbiased sample variance. 36 * 37 * - 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 38 * 39 * ```tex 40 * \Delta s^2 = s_{i+1}^2 - s_{i}^2 41 * ``` 42 * 43 * - If we multiply both sides by \\(N-1\\), 44 * 45 * ```tex 46 * (N-1)(\Delta s^2) = (N-1)s_{i+1}^2 - (N-1)s_{i}^2 47 * ``` 48 * 49 * - If we substitute the definition of the unbiased sample variance having the form 50 * 51 * ```tex 52 * \begin{align*} 53 * s^2 &= \frac{1}{N-1} \biggl( \sum_{i=1}^{N} (x_i - \bar{x})^2 \biggr) \\ 54 * &= \frac{1}{N-1} \biggl( \sum_{i=1}^{N} (x_i^2 - 2\bar{x}x_i + \bar{x}^2) \biggr) \\ 55 * &= \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) \\ 56 * &= \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) \\ 57 * &= \frac{1}{N-1} \biggl( \sum_{i=1}^{N} x_i^2 - 2N\bar{x}^2 + N\bar{x}^2 \biggr) \\ 58 * &= \frac{1}{N-1} \biggl( \sum_{i=1}^{N} x_i^2 - N\bar{x}^2 \biggr) 59 * \end{align*} 60 * ``` 61 * 62 * we return 63 * 64 * ```tex 65 * (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) 66 * ``` 67 * 68 * - This can be further simplified by recognizing that subtracting the sums reduces to \\(x_N^2 - x_0^2\\); in which case, 69 * 70 * ```tex 71 * \begin{align*} 72 * (N-1)(\Delta s^2) &= x_N^2 - x_0^2 - N\bar{x}_{i+1}^2 + N\bar{x}_{i}^2 \\ 73 * &= x_N^2 - x_0^2 - N(\bar{x}_{i+1}^2 - \bar{x}_{i}^2) \\ 74 * &= x_N^2 - x_0^2 - N(\bar{x}_{i+1} - \bar{x}_{i})(\bar{x}_{i+1} + \bar{x}_{i}) 75 * \end{align*} 76 * ``` 77 * 78 * - Recognizing that the difference of means can be expressed 79 * 80 * ```tex 81 * \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} 82 * ``` 83 * 84 * and substituting into the equation above 85 * 86 * ```tex 87 * (N-1)(\Delta s^2) = x_N^2 - x_0^2 - (x_N - x_0)(\bar{x}_{i+1} + \bar{x}_{i}) 88 * ``` 89 * 90 * - Rearranging terms gives us the update equation 91 * 92 * ```tex 93 * \begin{align*} 94 * (N-1)(\Delta s^2) &= (x_N - x_0)(x_N + x_0) - (x_N - x_0)(\bar{x}_{i+1} + \bar{x}_{i}) 95 * &= (x_N - x_0)(x_N + x_0 - \bar{x}_{i+1} - \bar{x}_{i}) \\ 96 * &= (x_N - x_0)(x_N - \bar{x}_{i+1} + x_0 - \bar{x}_{i}) 97 * \end{align*} 98 * ``` 99 * 100 * @param {PositiveInteger} W - window size 101 * @param {number} [mean] - mean value 102 * @throws {TypeError} first argument must be a positive integer 103 * @throws {TypeError} second argument must be a number primitive 104 * @returns {Function} accumulator function 105 * 106 * @example 107 * var accumulator = incrmvariance( 3 ); 108 * 109 * var s2 = accumulator(); 110 * // returns null 111 * 112 * s2 = accumulator( 2.0 ); 113 * // returns 0.0 114 * 115 * s2 = accumulator( -5.0 ); 116 * // returns 24.5 117 * 118 * s2 = accumulator( 3.0 ); 119 * // returns 19.0 120 * 121 * s2 = accumulator( 5.0 ); 122 * // returns 28.0 123 * 124 * s2 = accumulator(); 125 * // returns 28.0 126 * 127 * @example 128 * var accumulator = incrmvariance( 3, -2.0 ); 129 */ 130 function incrmvariance( W, mean ) { 131 var delta; 132 var buf; 133 var tmp; 134 var M2; 135 var mu; 136 var d1; 137 var d2; 138 var N; 139 var n; 140 var i; 141 if ( !isPositiveInteger( W ) ) { 142 throw new TypeError( 'invalid argument. Must provide a positive integer. Value: `' + W + '`.' ); 143 } 144 buf = new Array( W ); 145 n = W - 1; 146 M2 = 0.0; 147 i = -1; 148 N = 0; 149 if ( arguments.length > 1 ) { 150 if ( !isNumber( mean ) ) { 151 throw new TypeError( 'invalid argument. Must provide a number primitive. Value: `' + mean + '`.' ); 152 } 153 mu = mean; 154 return accumulator2; 155 } 156 mu = 0.0; 157 return accumulator1; 158 159 /** 160 * If provided a value, the accumulator function returns an updated unbiased sample variance. If not provided a value, the accumulator function returns the current unbiased sample variance. 161 * 162 * @private 163 * @param {number} [x] - input value 164 * @returns {(number|null)} unbiased sample variance or null 165 */ 166 function accumulator1( x ) { 167 var k; 168 var v; 169 if ( arguments.length === 0 ) { 170 if ( N === 0 ) { 171 return null; 172 } 173 if ( N === 1 ) { 174 return ( isnan( M2 ) ) ? NaN : 0.0; 175 } 176 if ( N < W ) { 177 return M2 / (N-1); 178 } 179 return M2 / n; 180 } 181 // Update the index for managing the circular buffer: 182 i = (i+1) % W; 183 184 // Case: incoming value is NaN, the sliding second moment is automatically NaN... 185 if ( isnan( x ) ) { 186 N = W; // explicitly set to avoid `N < W` branch 187 M2 = NaN; 188 } 189 // Case: initial window... 190 else if ( N < W ) { 191 buf[ i ] = x; // update buffer 192 N += 1; 193 delta = x - mu; 194 mu += delta / N; 195 M2 += delta * (x - mu); 196 if ( N === 1 ) { 197 return 0.0; 198 } 199 return M2 / (N-1); 200 } 201 // Case: N = W = 1 202 else if ( N === 1 ) { 203 M2 = 0.0; 204 return M2; 205 } 206 // Case: outgoing value is NaN, and, thus, we need to compute the accumulated values... 207 else if ( isnan( buf[ i ] ) ) { 208 N = 1; 209 mu = x; 210 M2 = 0.0; 211 for ( k = 0; k < W; k++ ) { 212 if ( k !== i ) { 213 v = buf[ k ]; 214 if ( isnan( v ) ) { 215 N = W; // explicitly set to avoid `N < W` branch 216 M2 = NaN; 217 break; // second moment is automatically NaN, so no need to continue 218 } 219 N += 1; 220 delta = v - mu; 221 mu += delta / N; 222 M2 += delta * (v - mu); 223 } 224 } 225 } 226 // Case: neither the current second moment nor the incoming value are NaN, so we need to update the accumulated values... 227 else if ( isnan( M2 ) === false ) { 228 tmp = buf[ i ]; 229 delta = x - tmp; 230 d1 = tmp - mu; 231 mu += delta / W; 232 d2 = x - mu; 233 M2 += delta * (d1 + d2); 234 } 235 // Case: the current second moment is NaN, so nothing to do until the buffer no longer contains NaN values... 236 237 buf[ i ] = x; 238 return M2 / n; 239 } 240 241 /** 242 * If provided a value, the accumulator function returns an updated unbiased sample variance. If not provided a value, the accumulator function returns the current unbiased sample variance. 243 * 244 * @private 245 * @param {number} [x] - input value 246 * @returns {(number|null)} unbiased sample variance or null 247 */ 248 function accumulator2( x ) { 249 var k; 250 if ( arguments.length === 0 ) { 251 if ( N === 0 ) { 252 return null; 253 } 254 if ( N < W ) { 255 return M2 / N; 256 } 257 return M2 / W; 258 } 259 // Update the index for managing the circular buffer: 260 i = (i+1) % W; 261 262 // Case: incoming value is NaN, the sliding second moment is automatically NaN... 263 if ( isnan( x ) ) { 264 N = W; // explicitly set to avoid `N < W` branch 265 M2 = NaN; 266 } 267 // Case: initial window... 268 else if ( N < W ) { 269 buf[ i ] = x; // update buffer 270 N += 1; 271 delta = x - mu; 272 M2 += delta * delta; 273 return M2 / N; 274 } 275 // Case: outgoing value is NaN, and, thus, we need to compute the accumulated values... 276 else if ( isnan( buf[ i ] ) ) { 277 M2 = 0.0; 278 for ( k = 0; k < W; k++ ) { 279 if ( k !== i ) { 280 if ( isnan( buf[ k ] ) ) { 281 N = W; // explicitly set to avoid `N < W` branch 282 M2 = NaN; 283 break; // second moment is automatically NaN, so no need to continue 284 } 285 delta = buf[ k ] - mu; 286 M2 += delta * delta; 287 } 288 } 289 } 290 // Case: neither the current second moment nor the incoming value are NaN, so we need to update the accumulated values... 291 else if ( isnan( M2 ) === false ) { 292 tmp = buf[ i ]; 293 M2 += ( x-tmp ) * ( x-mu + tmp-mu ); 294 } 295 // Case: the current second moment is NaN, so nothing to do until the buffer no longer contains NaN values... 296 297 buf[ i ] = x; 298 return M2 / W; 299 } 300 } 301 302 303 // EXPORTS // 304 305 module.exports = incrmvariance;