main.js (17310B)
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 sample Pearson product-moment correlation coefficient. 34 * 35 * ## Method 36 * 37 * - Let \\(W\\) be a window of \\(N\\) elements over which we want to compute a sample Pearson product-moment correlation coefficient. 38 * 39 * - We begin by defining the covariance \\( \operatorname{cov}_n(x,y) \\) for a window \\(n\\) as follows 40 * 41 * ```tex 42 * \operatorname{cov}_n(x,y) &= \frac{C_n}{n} 43 * ``` 44 * 45 * where \\(C_n\\) is the co-moment, which is defined as 46 * 47 * ```tex 48 * C_n = \sum_{i=1}^{N} ( x_i - \bar{x}_n ) ( y_i - \bar{y}_n ) 49 * ``` 50 * 51 * and where \\(\bar{x}_n\\) and \\(\bar{y}_n\\) are the sample means for \\(x\\) and \\(y\\), respectively, and \\(i=1\\) specifies the first element in a window. 52 * 53 * - The sample mean is computed using the canonical formula 54 * 55 * ```tex 56 * \bar{x}_n = \frac{1}{N} \sum_{i=1}^{N} x_i 57 * ``` 58 * 59 * which, taking into account a previous window, can be expressed 60 * 61 * ```tex 62 * \begin{align*} 63 * \bar{x}_n &= \frac{1}{N} \biggl( \sum_{i=0}^{N-1} x_i - x_0 + x_N \biggr) \\ 64 * &= \bar{x}_{n-1} + \frac{x_N - x_0}{N} 65 * \end{align*} 66 * ``` 67 * 68 * where \\(x_0\\) is the first value in the previous window. 69 * 70 * - We can substitute into the co-moment equation 71 * 72 * ```tex 73 * \begin{align*} 74 * C_n &= \sum_{i=1}^{N} ( x_i - \bar{x}_n ) ( y_i - \bar{y}_n ) \\ 75 * &= \sum_{i=1}^{N} \biggl( x_i - \bar{x}_{n-1} - \frac{x_N - x_0}{N} \biggr) \biggl( y_i - \bar{y}_{n-1} - \frac{y_N - y_0}{N} \biggr) \\ 76 * &= \sum_{i=1}^{N} \biggl( \Delta x_{i,n-1} - \frac{x_N - x_0}{N} \biggr) \biggl( \Delta y_{i,n-1} - \frac{y_N - y_0}{N} \biggr) 77 * \end{align*} 78 * ``` 79 * 80 * where 81 * 82 * ```tex 83 * \Delta x_{i,k} = x_i - \bar{x}_{k} 84 * ``` 85 * 86 * - We can subsequently expand terms and apply a summation identity 87 * 88 * ```tex 89 * \begin{align*} 90 * C_n &= \sum_{i=1}^{N} \Delta x_{i,n-1} \Delta y_{i,n-1} - \sum_{i=1}^{N} \Delta x_{i,n-1} \biggl( \frac{y_N - y_0}{N} \biggr) - \sum_{i=1}^{N} \Delta y_{i,n-1} \biggl( \frac{x_N - x_0}{N} \biggr) + \sum_{i=1}^{N} \biggl( \frac{x_N - x_0}{N} \biggr) \biggl( \frac{y_N - y_0}{N} \biggr) \\ 91 * &= \sum_{i=1}^{N} \Delta x_{i,n-1} \Delta y_{i,n-1} - \biggl( \frac{y_N - y_0}{N} \biggr) \sum_{i=1}^{N} \Delta x_{i,n-1} - \biggl( \frac{x_N - x_0}{N} \biggr) \sum_{i=1}^{N} \Delta y_{i,n-1} + \frac{(x_N - x_0)(y_N - y_0)}{N} 92 * \end{align*} 93 * ``` 94 * 95 * - Let us first consider the second term which we can reorganize as follows 96 * 97 * ```tex 98 * \begin{align*} 99 * \biggl( \frac{y_N - y_0}{N} \biggr) \sum_{i=1}^{N} \Delta x_{i,n-1} &= \biggl( \frac{y_N - y_0}{N} \biggr) \sum_{i=1}{N} ( x_i - \bar{x}_{n-1}) \\ 100 * &= \biggl( \frac{y_N - y_0}{N} \biggr) \sum_{i=1}^{N} x_i - \biggl( \frac{y_N - y_0}{N} \biggr) \sum_{i=1}^{N} \bar{x}_{n-1} \\ 101 * &= (y_N - y_0) \bar{x}_{n} - (y_N - y_0)\bar{x}_{n-1} \\ 102 * &= (y_N - y_0) (\bar{x}_{n} - \bar{x}_{n-1}) \\ 103 * &= \frac{(x_N - x_0)(y_N - y_0)}{N} 104 * \end{align*} 105 * ``` 106 * 107 * - The third term can be reorganized in a manner similar to the second term such that 108 * 109 * ```tex 110 * \biggl( \frac{x_N - x_0}{N} \biggr) \sum_{i=1}^{N} \Delta y_{i,n-1} = \frac{(x_N - x_0)(y_N - y_0)}{N} 111 * ``` 112 * 113 * - Substituting back into the equation for the co-moment 114 * 115 * ```tex 116 * \begin{align*} 117 * C_n &= \sum_{i=1}^{N} \Delta x_{i,n-1} \Delta y_{i,n-1} - \frac{(x_N - x_0)(y_N - y_0)}{N} - \frac{(x_N - x_0)(y_N - y_0)}{N} + \frac{(x_N - x_0)(y_N - y_0)}{N} \\ 118 * &= \sum_{i=1}^{N} \Delta x_{i,n-1} \Delta y_{i,n-1} - \frac{(x_N - x_0)(y_N - y_0)}{N} 119 * \end{align*} 120 * ``` 121 * 122 * - Let us now consider the first term which we can modify as follows 123 * 124 * ```tex 125 * \begin{align*} 126 * \sum_{i=1}^{N} \Delta x_{i,n-1} \Delta y_{i,n-1} &= \sum_{i=1}^{N} (x_i - \bar{x}_{n-1})(y_i - \bar{y}_{n-1}) \\ 127 * &= \sum_{i=1}^{N-1} (x_i - \bar{x}_{n-1})(y_i - \bar{y}_{n-1}) + (x_N - \bar{x}_{n-1})(y_N - \bar{y}_{n-1}) \\ 128 * &= \sum_{i=1}^{N-1} (x_i - \bar{x}_{n-1})(y_i - \bar{y}_{n-1}) + (x_N - \bar{x}_{n-1})(y_N - \bar{y}_{n-1}) + (x_0 - \bar{x}_{n-1})(y_0 - \bar{y}_{n-1}) - (x_0 - \bar{x}_{n-1})(y_0 - \bar{y}_{n-1}) \\ 129 * &= \sum_{i=0}^{N-1} (x_i - \bar{x}_{n-1})(y_i - \bar{y}_{n-1}) + (x_N - \bar{x}_{n-1})(y_N - \bar{y}_{n-1}) - (x_0 - \bar{x}_{n-1})(y_0 - \bar{y}_{n-1}) 130 * \end{align*} 131 * ``` 132 * 133 * where we recognize that the first term equals the co-moment for the previous window 134 * 135 * ```tex 136 * C_{n-1} = \sum_{i=0}^{N-1} (x_i - \bar{x}_{n-1})(y_i - \bar{y}_{n-1}) 137 * ``` 138 * 139 * In which case, 140 * 141 * ```tex 142 * \begin{align*} 143 * \sum_{i=1}^{N} \Delta x_{i,n-1} \Delta y_{i,n-1} &= C_{n-1} + (x_N - \bar{x}_{n-1})(y_N - \bar{y}_{n-1}) - (x_0 - \bar{x}_{n-1})(y_0 - \bar{y}_{n-1}) \\ 144 * &= C_{n-1} + \Delta x_{N,n-1} \Delta y_{N,n-1} - \Delta x_{0,n-1} \Delta y_{0,n-1} 145 * \end{align*} 146 * ``` 147 * 148 * - Substituting into the equation for the co-moment 149 * 150 * ```tex 151 * C_n = C_{n-1} + \Delta x_{N,n-1} \Delta y_{N,n-1} - \Delta x_{0,n-1} \Delta y_{0,n-1} - \frac{(x_N - x_0)(y_N - y_0)}{N} 152 * ``` 153 * 154 * - We can make one further modification to the last term 155 * 156 * ```tex 157 * \begin{align*} 158 * \frac{(x_N - x_0)(y_N - y_0)}{N} &= \frac{(x_N - \bar{x}_{n-1} - x_0 + \bar{x}_{n-1})(y_N - \bar{y}_{n-1} - y_0 + \bar{y}_{n-1})}{N} \\ 159 * &= \frac{(\Delta x_{N,n-1} - \Delta x_{0,n-1})(\Delta y_{N,n-1} - \Delta y_{0,n-1})}{N} 160 * \end{align*} 161 * ``` 162 * 163 * which, upon substitution into the equation for the co-moment, yields 164 * 165 * ```tex 166 * C_n = C_{n-1} + \Delta x_{N,n-1} \Delta y_{N,n-1} - \Delta x_{0,n-1} \Delta y_{0,n-1} - \frac{(\Delta x_{N,n-1} - \Delta x_{0,n-1})(\Delta y_{N,n-1} - \Delta y_{0,n-1})}{N} 167 * ``` 168 * 169 * - To calculate corrected sample standard deviations, 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. 170 * 171 * - The difference between the unbiased sample variance in a window \\(W_{n-1}\\) and the unbiased sample variance in a window \\(W_{n})\\) is given by 172 * 173 * ```tex 174 * \Delta s^2 = s_n^2 - s_{n-1}^2 175 * ``` 176 * 177 * - If we multiply both sides by \\(N-1\\), 178 * 179 * ```tex 180 * (N-1)(\Delta s^2) = (N-1)s_{n}^2 - (N-1)s_{n-1}^2 181 * ``` 182 * 183 * - If we substitute the definition of the unbiased sample variance having the form 184 * 185 * ```tex 186 * \begin{align*} 187 * s^2 &= \frac{1}{N-1} \biggl( \sum_{i=1}^{N} (x_i - \bar{x})^2 \biggr) \\ 188 * &= \frac{1}{N-1} \biggl( \sum_{i=1}^{N} (x_i^2 - 2\bar{x}x_i + \bar{x}^2) \biggr) \\ 189 * &= \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) \\ 190 * &= \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) \\ 191 * &= \frac{1}{N-1} \biggl( \sum_{i=1}^{N} x_i^2 - 2N\bar{x}^2 + N\bar{x}^2 \biggr) \\ 192 * &= \frac{1}{N-1} \biggl( \sum_{i=1}^{N} x_i^2 - N\bar{x}^2 \biggr) 193 * \end{align*} 194 * ``` 195 * 196 * we return 197 * 198 * ```tex 199 * (N-1)(\Delta s^2) = \biggl(\sum_{i=1}^N x_i^2 - N\bar{x}_{n}^2 \biggr) - \biggl(\sum_{i=0}^{N-1} x_i^2 - N\bar{x}_{n-1}^2 \biggr) 200 * ``` 201 * 202 * - This can be further simplified by recognizing that subtracting the sums reduces to \\(x_N^2 - x_0^2\\); in which case, 203 * 204 * ```tex 205 * \begin{align*} 206 * (N-1)(\Delta s^2) &= x_N^2 - x_0^2 - N\bar{x}_{n}^2 + N\bar{x}_{n-1}^2 \\ 207 * &= x_N^2 - x_0^2 - N(\bar{x}_{n}^2 - \bar{x}_{n-1}^2) \\ 208 * &= x_N^2 - x_0^2 - N(\bar{x}_{n} - \bar{x}_{n-1})(\bar{x}_{n} + \bar{x}_{n-1}) 209 * \end{align*} 210 * ``` 211 * 212 * - Recognizing that the difference of means can be expressed 213 * 214 * ```tex 215 * \bar{x}_{n} - \bar{x}_{n-1} = \frac{1}{N} \biggl( \sum_{i=1}^N x_i - \sum_{i=0}^{N-1} x_i \biggr) = \frac{x_N - x_0}{N} 216 * ``` 217 * 218 * and substituting into the equation above 219 * 220 * ```tex 221 * (N-1)(\Delta s^2) = x_N^2 - x_0^2 - (x_N - x_0)(\bar{x}_{n} + \bar{x}_{n-1}) 222 * ``` 223 * 224 * - Rearranging terms gives us the update equation 225 * 226 * ```tex 227 * \begin{align*} 228 * (N-1)(\Delta s^2) &= (x_N - x_0)(x_N + x_0) - (x_N - x_0)(\bar{x}_{n} + \bar{x}_{n-1}) 229 * &= (x_N - x_0)(x_N + x_0 - \bar{x}_{n} - \bar{x}_{n-1}) \\ 230 * &= (x_N - x_0)(x_N - \bar{x}_{n} + x_0 - \bar{x}_{n-1}) 231 * \end{align*} 232 * ``` 233 * 234 * - To compute the corrected sample standard deviation, we apply Bessel's correction and take the square root. 235 * 236 * - The sample Pearson product-moment correlation coefficient can thus be calculated as 237 * 238 * ```tex 239 * r_n(x,y) = \frac{\operatorname{cov}_n(x,y)}{\sigma_{x,n} \sigma_{y,n}} 240 * ``` 241 * 242 * where \\(\sigma_{x,n}\\) and \\(\sigma_{y,n}\\) are the corrected sample standard deviations for \\(x\\) and \\(y\\), respectively, for a window \\(n\\). 243 * 244 * 245 * @param {PositiveInteger} W - window size 246 * @param {number} [meanx] - mean value 247 * @param {number} [meany] - mean value 248 * @throws {TypeError} first argument must be a positive integer 249 * @throws {TypeError} second argument must be a number primitive 250 * @throws {TypeError} third argument must be a number primitive 251 * @returns {Function} accumulator function 252 * 253 * @example 254 * var accumulator = incrmpcorr( 3 ); 255 * 256 * var r = accumulator(); 257 * // returns null 258 * 259 * r = accumulator( 2.0, 1.0 ); 260 * // returns 0.0 261 * 262 * r = accumulator( -5.0, 3.14 ); 263 * // returns ~-1.0 264 * 265 * r = accumulator( 3.0, -1.0 ); 266 * // returns ~-0.925 267 * 268 * r = accumulator( 5.0, -9.5 ); 269 * // returns ~-0.863 270 * 271 * r = accumulator(); 272 * // returns ~-0.863 273 * 274 * @example 275 * var accumulator = incrmpcorr( 3, -2.0, 10.0 ); 276 */ 277 function incrmpcorr( W, meanx, meany ) { 278 var buf; 279 var dx0; 280 var dxN; 281 var dy0; 282 var dyN; 283 var M2x; 284 var M2y; 285 var mx; 286 var my; 287 var sx; 288 var sy; 289 var dx; 290 var dy; 291 var wi; 292 var C; 293 var N; 294 var n; 295 var i; 296 if ( !isPositiveInteger( W ) ) { 297 throw new TypeError( 'invalid argument. First argument must be a positive integer. Value: `' + W + '`.' ); 298 } 299 buf = new Float64Array( 2*W ); // strided array 300 n = W - 1; 301 M2x = 0.0; 302 M2y = 0.0; 303 C = 0.0; 304 wi = -1; 305 N = 0; 306 if ( arguments.length > 1 ) { 307 if ( !isNumber( meanx ) ) { 308 throw new TypeError( 'invalid argument. Second argument must be a number primitive. Value: `' + meanx + '`.' ); 309 } 310 if ( !isNumber( meany ) ) { 311 throw new TypeError( 'invalid argument. Third argument must be a number primitive. Value: `' + meany + '`.' ); 312 } 313 mx = meanx; 314 my = meany; 315 return accumulator2; 316 } 317 mx = 0.0; 318 my = 0.0; 319 return accumulator1; 320 321 /** 322 * If provided a value, the accumulator function returns an updated sample correlation coefficient. If not provided a value, the accumulator function returns the current sample correlation coefficient. 323 * 324 * @private 325 * @param {number} [x] - input value 326 * @param {number} [y] - input value 327 * @returns {(number|null)} sample correlation coefficient or null 328 */ 329 function accumulator1( x, y ) { 330 var v1; 331 var v2; 332 var n1; 333 var k; 334 var j; 335 if ( arguments.length === 0 ) { 336 if ( N === 0 ) { 337 return null; 338 } 339 if ( N === 1 ) { 340 return 0.0; 341 } 342 if ( N < W ) { 343 return ( C/(N-1) ) / ( sx*sy ); 344 } 345 return ( C/n ) / ( sx*sy ); 346 } 347 // Update the window and strided array indices for managing the circular buffer: 348 wi = (wi+1) % W; 349 i = 2 * wi; 350 351 // Case: an incoming value is NaN, the sliding co-moment is automatically NaN... 352 if ( isnan( x ) || isnan( y ) ) { 353 N = W; // explicitly set to avoid `N < W` branch 354 C = NaN; 355 } 356 // Case: initial window... 357 else if ( N < W ) { 358 buf[ i ] = x; // update buffer 359 buf[ i+1 ] = y; 360 361 N += 1; 362 363 dx = x - mx; 364 mx += dx / N; 365 M2x += dx * ( x-mx ); 366 367 dy = y - my; 368 my += dy / N; 369 dyN = y - my; 370 M2y += dy * dyN; 371 372 C += dx * dyN; 373 if ( N === 1 ) { 374 return 0.0; 375 } 376 n1 = N - 1; 377 sx = sqrt( M2x/n1 ); 378 sy = sqrt( M2y/n1 ); 379 return ( C/n1 ) / ( sx*sy ); // Note: why all the dividing by `N`? To avoid overflow. 380 } 381 // Case: N = W = 1 382 else if ( N === 1 ) { 383 return 0.0; 384 } 385 // Case: an outgoing value is NaN, and, thus, we need to compute the accumulated values... 386 else if ( isnan( buf[ i ] ) || isnan( buf[ i+1 ] ) ) { 387 N = 1; 388 mx = x; 389 my = y; 390 M2x = 0.0; 391 M2y = 0.0; 392 C = 0.0; 393 for ( k = 0; k < W; k++ ) { 394 j = 2 * k; // convert to a strided array index 395 if ( j !== i ) { 396 v1 = buf[ j ]; 397 v2 = buf[ j+1 ]; 398 if ( isnan( v1 ) || isnan( v2 ) ) { 399 N = W; // explicitly set to avoid `N < W` branch 400 C = NaN; 401 break; // co-moment is automatically NaN, so no need to continue 402 } 403 N += 1; 404 405 dx = v1 - mx; 406 mx += dx / N; 407 M2x += dx * ( v1-mx ); 408 409 dy = v2 - my; 410 my += dy / N; 411 dyN = v2 - my; 412 M2y += dy * dyN; 413 414 C += dx * dyN; 415 } 416 } 417 } 418 // Case: neither the current co-moment nor the incoming values are NaN, so we need to update the accumulated values... 419 else if ( isnan( C ) === false ) { 420 dx0 = buf[ i ] - mx; 421 dy0 = buf[ i+1 ] - my; 422 dxN = x - mx; 423 dyN = y - my; 424 dx = dxN - dx0; 425 dy = dyN - dy0; 426 mx += dx / W; 427 my += dy / W; 428 M2x += dx * ( dx0+(x-mx) ); 429 M2y += dy * ( dy0+(y-my) ); 430 C += (dxN*dyN) - (dx0*dy0) - ( dx*dy/W ); 431 } 432 // Case: the current co-moment is NaN, so nothing to do until the buffer no longer contains NaN values... 433 buf[ i ] = x; 434 buf[ i+1 ] = y; 435 436 sx = sqrt( M2x/n ); 437 sy = sqrt( M2y/n ); 438 return ( C/n ) / ( sx*sy ); // Note: why all the dividing by `n`? To avoid overflow. 439 } 440 441 /** 442 * If provided a value, the accumulator function returns an updated sample correlation coefficient. If not provided a value, the accumulator function returns the current sample correlation coefficient. 443 * 444 * @private 445 * @param {number} [x] - input value 446 * @param {number} [y] - input value 447 * @returns {(number|null)} sample correlation coefficient or null 448 */ 449 function accumulator2( x, y ) { 450 var k; 451 var j; 452 if ( arguments.length === 0 ) { 453 if ( N === 0 ) { 454 return null; 455 } 456 if ( N < W ) { 457 return ( C/N ) / ( sx*sy ); 458 } 459 return ( C/W ) / ( sx*sy ); 460 } 461 // Update the window and strided array indices for managing the circular buffer: 462 wi = (wi+1) % W; 463 i = 2 * wi; 464 465 // Case: an incoming value is NaN, the sliding co-moment is automatically NaN... 466 if ( isnan( x ) || isnan( y ) ) { 467 N = W; // explicitly set to avoid `N < W` branch 468 C = NaN; 469 } 470 // Case: initial window... 471 else if ( N < W ) { 472 buf[ i ] = x; // update buffer 473 buf[ i+1 ] = y; 474 475 N += 1; 476 dx = x - mx; 477 M2x += dx * dx; 478 dy = y - my; 479 M2y += dy * dy; 480 481 C += dx * dy; 482 sx = sqrt( M2x/N ); 483 sy = sqrt( M2y/N ); 484 return ( C/N ) / ( sx*sy ); // Note: why all the dividing by `N`? To avoid overflow. 485 } 486 // Case: an outgoing value is NaN, and, thus, we need to compute the accumulated values... 487 else if ( isnan( buf[ i ] ) || isnan( buf[ i+1 ] ) ) { 488 M2x = 0.0; 489 M2y = 0.0; 490 C = 0.0; 491 for ( k = 0; k < W; k++ ) { 492 j = 2 * k; // convert to a strided array index 493 if ( j !== i ) { 494 if ( isnan( buf[ j ] ) || isnan( buf[ j+1 ] ) ) { 495 N = W; // explicitly set to avoid `N < W` branch 496 C = NaN; 497 break; // co-moment is automatically NaN, so no need to continue 498 } 499 dx = buf[j] - mx; 500 M2x += dx * dx; 501 dy = buf[j+1] - my; 502 M2y += dy * dy; 503 C += dx * dy; 504 } 505 } 506 } 507 // Case: neither the current co-moment nor the incoming values are NaN, so we need to update the accumulated values... 508 else if ( isnan( C ) === false ) { 509 // Use textbook formulas along with simplification arising from difference of sums: 510 dx0 = buf[ i ] - mx; 511 dxN = x - mx; 512 dy0 = buf[ i+1 ] - my; 513 dyN = y - my; 514 M2x += ( dxN-dx0 ) * ( dxN+dx0 ); 515 M2y += ( dyN-dy0 ) * ( dyN+dy0 ); 516 C += ( dxN*dyN ) - ( dx0*dy0 ); 517 } 518 // Case: the current co-moment is NaN, so nothing to do until the buffer no longer contains NaN values... 519 buf[ i ] = x; 520 buf[ i+1 ] = y; 521 522 sx = sqrt( M2x/W ); 523 sy = sqrt( M2y/W ); 524 return ( C/W ) / ( sx*sy ); // Note: why all the dividing by `W`? To avoid overflow. 525 } 526 } 527 528 529 // EXPORTS // 530 531 module.exports = incrmpcorr;