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