|
| 1 | +/** |
| 2 | +* @license Apache-2.0 |
| 3 | +* |
| 4 | +* Copyright (c) 2025 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 sqrt = require( '@stdlib/math/base/special/sqrt' ); |
| 25 | +var format = require( '@stdlib/string/format' ); |
| 26 | + |
| 27 | + |
| 28 | +// MAIN // |
| 29 | + |
| 30 | +/** |
| 31 | +* Returns an accumulator function which incrementally computes a weighted standard deviation. |
| 32 | +* |
| 33 | +* ## Method |
| 34 | +* |
| 35 | +* - The weighted standard deviation is defined as |
| 36 | +* |
| 37 | +* ```tex |
| 38 | +* \sigma = \sqrt{\frac{\sum_{i=0}^{n-1} w_i (x_i - \mu)^2}{\sum_{i=0}^{n-1} w_i}} |
| 39 | +* ``` |
| 40 | +* |
| 41 | +* where \( w_i \) are the weights and \( \mu \) is the weighted mean. |
| 42 | +* |
| 43 | +* - To derive an incremental formula for computing a weighted standard deviation, let: |
| 44 | +* |
| 45 | +* ```tex |
| 46 | +* W_n = \sum_{i=1}^{n} w_i |
| 47 | +* ``` |
| 48 | +* |
| 49 | +* - Accordingly: |
| 50 | +* |
| 51 | +* ```tex |
| 52 | +* \begin{align*} |
| 53 | +* \sigma^2 &= \frac{1}{W_n} \sum_{i=1}^{n} w_i (x_i - \mu)^2 \\ |
| 54 | +* &= \frac{1}{W_n} \sum_{i=1}^{n} w_i x_i^2 - \mu^2 |
| 55 | +* \end{align*} |
| 56 | +* ``` |
| 57 | +* |
| 58 | +* - Let: |
| 59 | +* |
| 60 | +* ```tex |
| 61 | +* M2_n = W_n \sigma^2_n = \sum_{i=1}^{n} w_i x_i^2 - W_n \mu_n^2 |
| 62 | +* ``` |
| 63 | +* |
| 64 | +* - To compute the incremental update for \(M2_n\): |
| 65 | +* |
| 66 | +* ```tex |
| 67 | +* \begin{align*} |
| 68 | +* M2_n - M2_{n-1} &= \sum_{i=1}^{n} w_i x_i^2 - W_n \mu_n^2 - \sum_{i=1}^{n-1} w_i x_i^2 + W_{n-1} \mu_{n-1}^2 \\ |
| 69 | +* &= w_n x_n^2 - W_n \mu_n^2 + W_{n-1} \mu_{n-1}^2 \\ |
| 70 | +* &= w_n x_n^2 - W_n \mu_n^2 + (W_n - w_n) \mu_{n-1}^2 \\ |
| 71 | +* &= w_n (x_n^2 - \mu_{n-1}^2) + W_n (\mu_{n-1}^2 - \mu_n^2) |
| 72 | +* \end{align*} |
| 73 | +* ``` |
| 74 | +* |
| 75 | +* - Simplify further: |
| 76 | +* |
| 77 | +* ```tex |
| 78 | +* \begin{align*} |
| 79 | +* M2_n - M2_{n-1} &= w_n (x_n - \mu_{n-1}) (x_n - \mu_n) \\ |
| 80 | +* M2_n &= M2_{n-1} + w_n (x_n - \mu_{n-1}) (x_n - \mu_n) |
| 81 | +* \end{align*} |
| 82 | +* ``` |
| 83 | +* |
| 84 | +* - Finally, compute the weighted standard deviation: |
| 85 | +* |
| 86 | +* ```tex |
| 87 | +* \sigma_n = \sqrt{\frac{M2_n}{W_n}} |
| 88 | +* ``` |
| 89 | +* |
| 90 | +* - Incrementally updating the mean: |
| 91 | +* |
| 92 | +* ```tex |
| 93 | +* \begin{align*} |
| 94 | +* \mu_n &= \mu_{n-1} + \frac{w_n}{W_n} (x_n - \mu_{n-1}) |
| 95 | +* \end{align*} |
| 96 | +* ``` |
| 97 | +* |
| 98 | +* Notation: |
| 99 | +* |
| 100 | +* - \(W_n = \sum_{i=1}^{n} w_i\) is the cumulative sum of weights up to \(n\). |
| 101 | +* - \(\mu_n = \frac{1}{W_n} \sum_{i=1}^{n} w_i x_i\) is the weighted mean at step \(n\). |
| 102 | +* - \(M2_n\) is the cumulative sum of weighted squared deviations. |
| 103 | +* |
| 104 | +* @param {number} [mean] - mean value |
| 105 | +* @throws {TypeError} must provide a number primitive |
| 106 | +* @returns {Function} accumulator function |
| 107 | +* |
| 108 | +* @example |
| 109 | +* var accumulator = incrwstdev(); |
| 110 | +* |
| 111 | +* var stdev = accumulator(); |
| 112 | +* // returns null |
| 113 | +* |
| 114 | +* stdev = accumulator( 2.0, 1.0 ); |
| 115 | +* // returns 0.0 |
| 116 | +* |
| 117 | +* stdev = accumulator( 3.0, 2.0 ); |
| 118 | +* // returns ~0.471 |
| 119 | +* |
| 120 | +* stdev = accumulator( 2.0, 0.1 ); |
| 121 | +* // returns ~0.478 |
| 122 | +* |
| 123 | +* stdev = accumulator(); |
| 124 | +* // returns ~0.478 |
| 125 | +* |
| 126 | +* @example |
| 127 | +* var accumulator = incrwstdev( 3.0 ); |
| 128 | +*/ |
| 129 | +function incrwstdev( mean ) { |
| 130 | + var delta; |
| 131 | + var wsum; |
| 132 | + var FLG; |
| 133 | + var M2; |
| 134 | + var mu; |
| 135 | + |
| 136 | + wsum = 0.0; |
| 137 | + M2 = 0.0; |
| 138 | + if ( arguments.length ) { |
| 139 | + if ( !isNumber( mean ) ) { |
| 140 | + throw new TypeError( format( 'invalid argument. Must provide a number. Value: `%s`.', mean ) ); |
| 141 | + } |
| 142 | + mu = mean; |
| 143 | + return accumulator2; |
| 144 | + } |
| 145 | + mu = 0.0; |
| 146 | + return accumulator1; |
| 147 | + |
| 148 | + /** |
| 149 | + * 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. |
| 150 | + * |
| 151 | + * @private |
| 152 | + * @param {number} [x] - value |
| 153 | + * @param {number} [w] - weight |
| 154 | + * @returns {(number|null)} weighted standard deviation or null |
| 155 | + */ |
| 156 | + function accumulator1( x, w ) { |
| 157 | + if ( arguments.length === 0 ) { |
| 158 | + if ( FLG === void 0 ) { |
| 159 | + return null; |
| 160 | + } |
| 161 | + return sqrt( M2 / wsum ); |
| 162 | + } |
| 163 | + FLG = true; |
| 164 | + delta = x - mu; |
| 165 | + wsum += w; |
| 166 | + mu += ( w / wsum ) * delta; |
| 167 | + M2 += w * delta * ( x - mu ); |
| 168 | + return sqrt( M2 / wsum ); |
| 169 | + } |
| 170 | + |
| 171 | + /** |
| 172 | + * 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. |
| 173 | + * |
| 174 | + * @private |
| 175 | + * @param {number} [x] - value |
| 176 | + * @param {number} [w] - weight |
| 177 | + * @returns {(number|null)} weighted standard deviation or null |
| 178 | + */ |
| 179 | + function accumulator2( x, w ) { |
| 180 | + if ( arguments.length === 0 ) { |
| 181 | + if ( FLG === void 0 ) { |
| 182 | + return null; |
| 183 | + } |
| 184 | + return sqrt( M2 / wsum ); |
| 185 | + } |
| 186 | + FLG = true; |
| 187 | + delta = x - mu; |
| 188 | + wsum += w; |
| 189 | + M2 += w * delta * delta; |
| 190 | + return sqrt( M2 / wsum ); |
| 191 | + } |
| 192 | +} |
| 193 | + |
| 194 | + |
| 195 | +// EXPORTS // |
| 196 | + |
| 197 | +module.exports = incrwstdev; |
0 commit comments