main.js (9170B)
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 coefficient of variation (CV). 34 * 35 * ## Method 36 * 37 * - Let \\(W\\) be a window of \\(N\\) elements over which we want to compute the coefficient of variation (CV), which is defined as the sample standard deviation divided by the sample mean. 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 for the unbiased sample variance 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 = incrmcv( 3 ); 112 * 113 * var cv = accumulator(); 114 * // returns null 115 * 116 * cv = accumulator( 2.0 ); 117 * // returns 0.0 118 * 119 * cv = accumulator( 1.0 ); 120 * // returns ~0.47 121 * 122 * cv = accumulator( 3.0 ); 123 * // returns 0.5 124 * 125 * cv = accumulator( 7.0 ); 126 * // returns ~0.83 127 * 128 * cv = accumulator(); 129 * // returns ~0.83 130 * 131 * @example 132 * var accumulator = incrmcv( 3, 2.0 ); 133 */ 134 function incrmcv( 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 accumulated value. If not provided a value, the accumulator function returns the current accumulated value. 165 * 166 * @private 167 * @param {number} [x] - input value 168 * @returns {(number|null)} accumulated value 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/mu; 179 } 180 if ( N < W ) { 181 return sqrt( M2/(N-1) ) / mu; 182 } 183 return sqrt( M2/n ) / mu; 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 mu = NaN; 192 M2 = NaN; 193 } 194 // Case: initial window... 195 else if ( N < W ) { 196 buf[ i ] = x; // update buffer 197 N += 1; 198 delta = x - mu; 199 mu += delta / N; 200 M2 += delta * (x - mu); 201 if ( N === 1 ) { 202 return 0.0 / mu; 203 } 204 return sqrt( M2/(N-1) ) / mu; 205 } 206 // Case: N = W = 1 207 else if ( N === 1 ) { 208 mu = x; 209 M2 = 0.0; 210 return M2 / mu; 211 } 212 // Case: outgoing value is NaN, and, thus, we need to compute the accumulated values... 213 else if ( isnan( buf[ i ] ) ) { 214 N = 1; 215 mu = x; 216 M2 = 0.0; 217 for ( k = 0; k < W; k++ ) { 218 if ( k !== i ) { 219 v = buf[ k ]; 220 if ( isnan( v ) ) { 221 N = W; // explicitly set to avoid `N < W` branch 222 mu = NaN; 223 M2 = NaN; 224 break; // second moment is automatically NaN, so no need to continue 225 } 226 N += 1; 227 delta = v - mu; 228 mu += delta / N; 229 M2 += delta * (v - mu); 230 } 231 } 232 } 233 // Case: neither the current second moment nor the incoming value are NaN, so we need to update the accumulated values... 234 else if ( isnan( M2 ) === false ) { 235 tmp = buf[ i ]; 236 delta = x - tmp; 237 d1 = tmp - mu; 238 mu += delta / W; 239 d2 = x - mu; 240 M2 += delta * (d1 + d2); 241 } 242 // Case: the current second moment is NaN, so nothing to do until the buffer no longer contains NaN values... 243 buf[ i ] = x; 244 245 return sqrt( M2/n ) / mu; 246 } 247 248 /** 249 * If provided a value, the accumulator function returns an updated accumulated value. If not provided a value, the accumulator function returns the current accumulated value. 250 * 251 * @private 252 * @param {number} [x] - input value 253 * @returns {(number|null)} accumulated value or null 254 */ 255 function accumulator2( x ) { 256 var k; 257 if ( arguments.length === 0 ) { 258 if ( N === 0 ) { 259 return null; 260 } 261 if ( N < W ) { 262 return sqrt( M2/N ) / mu; 263 } 264 return sqrt( M2/W ) / mu; 265 } 266 // Update the index for managing the circular buffer: 267 i = (i+1) % W; 268 269 // Case: incoming value is NaN, the sliding second moment is automatically NaN... 270 if ( isnan( x ) ) { 271 N = W; // explicitly set to avoid `N < W` branch 272 M2 = NaN; 273 } 274 // Case: initial window... 275 else if ( N < W ) { 276 buf[ i ] = x; // update buffer 277 N += 1; 278 delta = x - mu; 279 M2 += delta * delta; 280 return sqrt( M2/N ) / mu; 281 } 282 // Case: outgoing value is NaN, and, thus, we need to compute the accumulated values... 283 else if ( isnan( buf[ i ] ) ) { 284 M2 = 0.0; 285 for ( k = 0; k < W; k++ ) { 286 if ( k !== i ) { 287 if ( isnan( buf[ k ] ) ) { 288 N = W; // explicitly set to avoid `N < W` branch 289 M2 = NaN; 290 break; // second moment is automatically NaN, so no need to continue 291 } 292 delta = buf[ k ] - mu; 293 M2 += delta * delta; 294 } 295 } 296 } 297 // Case: neither the current second moment nor the incoming value are NaN, so we need to update the accumulated values... 298 else if ( isnan( M2 ) === false ) { 299 tmp = buf[ i ]; 300 M2 += ( x-tmp ) * ( x-mu + tmp-mu ); 301 } 302 // Case: the current second moment is NaN, so nothing to do until the buffer no longer contains NaN values... 303 buf[ i ] = x; 304 305 return sqrt( M2/W ) / mu; 306 } 307 } 308 309 310 // EXPORTS // 311 312 module.exports = incrmcv;