main.js (5265B)
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 isNumber = require( '@stdlib/assert/is-number' ).isPrimitive; 24 var isnan = require( '@stdlib/math/base/assert/is-nan' ); 25 26 27 // MAIN // 28 29 /** 30 * Returns an accumulator function which incrementally computes an unbiased sample covariance. 31 * 32 * ## Method 33 * 34 * - We begin by defining the co-moment \\(C_n\\) 35 * 36 * ```tex 37 * C_n = \sum_{i=1}^{N} ( x_i - \bar{x}_n ) ( y_i - \bar{y}_n ) 38 * ``` 39 * 40 * where \\(\bar{x}_n\\) and \\(\bar{y}_n\\) are the sample means for \\(x\\) and \\(y\\), respectively. 41 * 42 * - Based on Welford's method, we know the update formulas for the sample means are given by 43 * 44 * ```tex 45 * \bar{x}_n = \bar{x}_{n-1} + \frac{x_n - \bar{x}_{n-1}}{n} 46 * ``` 47 * 48 * and 49 * 50 * ```tex 51 * \bar{y}_n = \bar{y}_{n-1} + \frac{y_n - \bar{y}_{n-1}}{n} 52 * ``` 53 * 54 * - Substituting into the equation for \\(C_n\\) and rearranging terms 55 * 56 * ```tex 57 * C_n = C_{n-1} + (x_n - \bar{x}_n) (y_n - \bar{y}_{n-1}) 58 * ``` 59 * 60 * where the apparent asymmetry arises from 61 * 62 * ```tex 63 * x_n - \bar{x}_n = \frac{n-1}{n} (x_n - \bar{x}_{n-1}) 64 * ``` 65 * 66 * and, hence, the update term can be equivalently expressed 67 * 68 * ```tex 69 * \frac{n-1}{n} (x_n - \bar{x}_{n-1}) (y_n - \bar{y}_{n-1}) 70 * ``` 71 * 72 * - The covariance can be defined 73 * 74 * ```tex 75 * \begin{align*} 76 * \operatorname{cov}_n(x,y) &= \frac{C_n}{n} \\ 77 * &= \frac{C_{n-1} + (x_n - \bar{x}_n) (y_n - \bar{y}_{n-1})}{n} \\ 78 * &= \frac{(n-1)\operatorname{cov}_{n-1}(x,y) + (x_n - \bar{x}_n) (y_n - \bar{y}_{n-1})}{n} 79 * \end{align*} 80 * ``` 81 * 82 * - Applying Bessel's correction, we arrive at an update formula for calculating an unbiased sample covariance 83 * 84 * ```tex 85 * \begin{align*} 86 * \operatorname{cov}_n(x,y) &= \frac{n}{n-1}\cdot\frac{(n-1)\operatorname{cov}_{n-1}(x,y) + (x_n - \bar{x}_n) (y_n - \bar{y}_{n-1})}{n} \\ 87 * &= \operatorname{cov}_{n-1}(x,y) + \frac{(x_n - \bar{x}_n) (y_n - \bar{y}_{n-1})}{n-1} \\ 88 * &= \frac{C_{n-1} + (x_n - \bar{x}_n) (y_n - \bar{y}_{n-1})}{n-1} 89 * &= \frac{C_{n-1} + (x_n - \bar{x}_{n-1}) (y_n - \bar{y}_n)}{n-1} 90 * \end{align*} 91 * ``` 92 * 93 * @param {number} [meanx] - mean value 94 * @param {number} [meany] - mean value 95 * @throws {TypeError} first argument must be a number primitive 96 * @throws {TypeError} second argument must be a number primitive 97 * @returns {Function} accumulator function 98 * 99 * @example 100 * var accumulator = incrcovariance(); 101 * 102 * var v = accumulator(); 103 * // returns null 104 * 105 * v = accumulator( 2.0, 1.0 ); 106 * // returns 0.0 107 * 108 * v = accumulator( -5.0, 3.14 ); 109 * // returns ~-7.49 110 * 111 * v = accumulator(); 112 * // returns ~-7.49 113 * 114 * @example 115 * var accumulator = incrcovariance( 2.0, -3.0 ); 116 */ 117 function incrcovariance( meanx, meany ) { 118 var dx; 119 var mx; 120 var my; 121 var C; 122 var N; 123 124 C = 0.0; 125 N = 0; 126 if ( arguments.length ) { 127 if ( !isNumber( meanx ) ) { 128 throw new TypeError( 'invalid argument. First argument must be a number primitive. Value: `' + meanx + '`.' ); 129 } 130 if ( !isNumber( meany ) ) { 131 throw new TypeError( 'invalid argument. Second argument must be a number primitive. Value: `' + meany + '`.' ); 132 } 133 mx = meanx; 134 my = meany; 135 return accumulator2; 136 } 137 mx = 0.0; 138 my = 0.0; 139 return accumulator1; 140 141 /** 142 * If provided input values, the accumulator function returns an updated unbiased sample covariance. If not provided input values, the accumulator function returns the current unbiased sample covariance. 143 * 144 * @private 145 * @param {number} [x] - new value 146 * @param {number} [y] - new value 147 * @returns {(number|null)} unbiased sample covariance or null 148 */ 149 function accumulator1( x, y ) { 150 if ( arguments.length === 0 ) { 151 if ( N === 0 ) { 152 return null; 153 } 154 if ( N === 1 ) { 155 return ( isnan( C ) ) ? NaN : 0.0; 156 } 157 return C / (N-1); 158 } 159 N += 1; 160 dx = x - mx; 161 mx += dx / N; 162 my += ( y-my ) / N; 163 C += dx * ( y-my ); // Note: repeated `y-my` is intentional, as `my` is updated when used here 164 if ( N < 2 ) { 165 return ( isnan( C ) ) ? NaN : 0.0; 166 } 167 return C / (N-1); 168 } 169 170 /** 171 * If provided input values, the accumulator function returns an updated unbiased sample covariance. If not provided input values, the accumulator function returns the current unbiased sample covariance. 172 * 173 * @private 174 * @param {number} [x] - new value 175 * @param {number} [y] - new value 176 * @returns {(number|null)} unbiased sample covariance or null 177 */ 178 function accumulator2( x, y ) { 179 if ( arguments.length === 0 ) { 180 if ( N === 0 ) { 181 return null; 182 } 183 return C / N; 184 } 185 N += 1; 186 C += ( x-mx ) * ( y-my ); 187 return C / N; 188 } 189 } 190 191 192 // EXPORTS // 193 194 module.exports = incrcovariance;