main.js (8637B)
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 hasOwnProp = require( '@stdlib/assert/has-own-property' ); 24 var isObject = require( '@stdlib/assert/is-plain-object' ); 25 var isPositiveInteger = require( '@stdlib/assert/is-positive-integer' ).isPrimitive; 26 var isBoolean = require( '@stdlib/assert/is-boolean' ).isPrimitive; 27 var copy = require( '@stdlib/utils/copy' ); 28 var setReadOnly = require( '@stdlib/utils/define-read-only-property' ); 29 var setReadOnlyAccessor = require( '@stdlib/utils/define-read-only-accessor' ); 30 var max = require( '@stdlib/math/base/special/max' ); 31 var sqrt = require( '@stdlib/math/base/special/sqrt' ); 32 var roundn = require( '@stdlib/math/base/special/roundn' ); 33 var tQuantile = require( './../../../base/dists/t/quantile' ); 34 var Float64Array = require( '@stdlib/array/float64' ); 35 var validate = require( './validate.js' ); 36 var defaults = require( './defaults.json' ); 37 var incrmminmax = require( './minmax.js' ); 38 var incrmmeanstdev = require( './meanstdev.js' ); 39 40 41 // MAIN // 42 43 /** 44 * Returns an accumulator function which incrementally performs a moving Grubbs' test for detecting outliers. 45 * 46 * @param {PositiveInteger} W - window size 47 * @param {Options} [options] - function options 48 * @param {number} [options.alpha=0.05] - significance level 49 * @param {string} [options.alternative='two-sided'] - alternative hypothesis ('two-sided', 'min', 'max') 50 * @throws {TypeError} first argument must be a positive integer 51 * @throws {RangeError} first argument must be greater than or equal to 3 52 * @throws {TypeError} options argument must be an object 53 * @throws {TypeError} must provide valid options 54 * @throws {RangeError} `alpha` option must be on the interval `[0,1]` 55 * @returns {Function} accumulator function 56 * 57 * @example 58 * var rnorm = require( '@stdlib/random/base/normal' ); 59 * 60 * var accumulator; 61 * var opts; 62 * var i; 63 * 64 * accumulator = incrmgrubbs( 20, opts ); 65 * 66 * for ( i = 0; i < 200; i++ ) { 67 * res = accumulator( rnorm( 10.0, 5.0 ) ); 68 * } 69 */ 70 function incrmgrubbs( W ) { 71 var meanstdev; 72 var results; 73 var minmax; 74 var opts; 75 var err; 76 var buf; 77 var sig; 78 var mm; 79 var ms; 80 var tc; 81 var gc; 82 var df; 83 var N; 84 var G; 85 var i; 86 87 if ( !isPositiveInteger( W ) ) { 88 throw new TypeError( 'invalid argument. Window size must be a positive integer. Value: `' + W + '`.' ); 89 } 90 if ( W < 3 ) { 91 throw new RangeError( 'invalid argument. Window size must be greater than or equal to 3. Value: `' + W + '`.' ); 92 } 93 opts = copy( defaults ); 94 if ( arguments.length > 1 ) { 95 err = validate( opts, arguments[ 1 ] ); 96 if ( err ) { 97 throw err; 98 } 99 } 100 buf = new Float64Array( W ); 101 df = W - 2; 102 gc = 0.0; 103 G = 0.0; 104 N = 0; 105 i = -1; 106 107 // Compute the critical values: 108 if ( opts.alternative === 'min' ) { 109 sig = opts.alpha / W; 110 } else if ( opts.alternative === 'max' ) { 111 sig = opts.alpha / W; 112 } else { // two-sided 113 sig = opts.alpha / (2*W); 114 } 115 tc = tQuantile( 1.0-sig, df ); 116 gc = (W-1)*tc / sqrt( W*(df+(tc*tc)) ); 117 118 // Initialize statistics accumulators: 119 mm = [ 0.0, 0.0 ]; 120 minmax = incrmminmax( mm, W, buf ); 121 122 ms = [ 0.0, 0.0 ]; 123 meanstdev = incrmmeanstdev( ms, W, buf ); 124 125 // Initialize the results object: 126 results = {}; 127 setReadOnlyAccessor( results, 'rejected', getRejected ); 128 setReadOnly( results, 'alpha', opts.alpha ); 129 setReadOnly( results, 'criticalValue', gc ); 130 setReadOnlyAccessor( results, 'statistic', getStatistic ); 131 setReadOnly( results, 'df', df ); 132 setReadOnlyAccessor( results, 'mean', getMean ); 133 setReadOnlyAccessor( results, 'sd', getStDev ); 134 setReadOnlyAccessor( results, 'min', getMin ); 135 setReadOnlyAccessor( results, 'max', getMax ); 136 setReadOnly( results, 'alt', opts.alternative ); 137 setReadOnly( results, 'method', 'Grubbs\' Test' ); 138 setReadOnly( results, 'print', print ); 139 140 return accumulator; 141 142 /** 143 * If provided a value, the accumulator function returns updated Grubbs' test results. If not provided a value, the accumulator function returns the current Grubbs' test results. 144 * 145 * @private 146 * @param {number} [x] - new value 147 * @returns {(Object|null)} test results or null 148 */ 149 function accumulator( x ) { 150 var md; 151 if ( arguments.length === 0 ) { 152 if ( N < W ) { 153 return null; 154 } 155 return results; 156 } 157 N += 1; 158 159 // Update the index for managing the circular buffer: 160 i = (i+1) % W; 161 162 // Update model statistics: 163 meanstdev( x, i ); 164 minmax( x, i ); 165 166 // Insert the value into the buffer: 167 buf[ i ] = x; 168 169 if ( N < W ) { 170 return null; 171 } 172 // Compute the test statistic... 173 if ( opts.alternative === 'min' ) { 174 G = ( ms[0]-mm[0] ) / ms[ 1 ]; 175 } else if ( opts.alternative === 'max' ) { 176 G = ( mm[1]-ms[0] ) / ms[ 1 ]; 177 } else { // two-sided 178 md = max( ms[0]-mm[0], mm[1]-ms[0] ); // maximum absolute deviation 179 G = md / ms[ 1 ]; 180 } 181 return results; 182 } 183 184 /** 185 * Returns a `boolean` indicating whether the null hypothesis should be rejected. 186 * 187 * @private 188 * @returns {boolean} boolean indicating whether the null hypothesis should be rejected 189 */ 190 function getRejected() { 191 return ( G > gc ); 192 } 193 194 /** 195 * Returns the test statistic. 196 * 197 * @private 198 * @returns {number} test statistic 199 */ 200 function getStatistic() { 201 return G; 202 } 203 204 /** 205 * Returns the sample mean. 206 * 207 * @private 208 * @returns {number} sample mean 209 */ 210 function getMean() { 211 return ms[ 0 ]; 212 } 213 214 /** 215 * Returns the corrected sample standard deviation. 216 * 217 * @private 218 * @returns {number} corrected sample standard deviation 219 */ 220 function getStDev() { 221 return ms[ 1 ]; 222 } 223 224 /** 225 * Returns the sample minimum. 226 * 227 * @private 228 * @returns {number} sample minimum 229 */ 230 function getMin() { 231 return mm[ 0 ]; 232 } 233 234 /** 235 * Returns the sample maximum. 236 * 237 * @private 238 * @returns {number} sample maximum 239 */ 240 function getMax() { 241 return mm[ 1 ]; 242 } 243 244 /** 245 * Pretty-print test results. 246 * 247 * @private 248 * @param {Object} [options] - options object 249 * @param {PositiveInteger} [options.digits=4] - number of digits after the decimal point 250 * @param {boolean} [options.decision=true] - boolean indicating whether to print the test decision 251 * @throws {TypeError} options argument must be an object 252 * @throws {TypeError} must provide valid options 253 * @returns {string} formatted output 254 */ 255 function print( options ) { 256 var decision; 257 var digits; 258 var str; 259 260 digits = opts.digits; 261 decision = opts.decision; 262 if ( arguments.length > 0 ) { 263 if ( !isObject( options ) ) { 264 throw new TypeError( 'invalid argument. Must provide an object. Value: `' + options + '`.' ); 265 } 266 if ( hasOwnProp( options, 'digits' ) ) { 267 if ( !isPositiveInteger( options.digits ) ) { 268 throw new TypeError( 'invalid option. `digits` option must be a positive integer. Option: `' + options.digits + '`.' ); 269 } 270 digits = options.digits; 271 } 272 if ( hasOwnProp( options, 'decision' ) ) { 273 if ( !isBoolean( options.decision ) ) { 274 throw new TypeError( 'invalid option. `decision` option must be boolean. Option: `' + options.decision + '`.' ); 275 } 276 decision = options.decision; 277 } 278 } 279 str = ''; 280 str += results.method; 281 str += '\n\n'; 282 str += 'Alternative hypothesis: '; 283 if ( opts.alternative === 'max' ) { 284 str += 'The maximum value (' + mm[ 1 ] + ') is an outlier'; 285 } else if ( opts.alternative === 'min' ) { 286 str += 'The minimum value (' + mm[ 0 ] + ') is an outlier'; 287 } else { // two-sided 288 str += 'The '; 289 if ( ms[0]-mm[0] > mm[1]-ms[0] ) { 290 str += 'minimum value (' + mm[ 0 ] + ')'; 291 } else { 292 str += 'maximum value (' + mm[ 1 ] + ')'; 293 } 294 str += ' is an outlier'; 295 } 296 str += '\n\n'; 297 str += ' criticalValue: ' + roundn( gc, -digits ) + '\n'; 298 str += ' statistic: ' + roundn( G, -digits ) + '\n'; 299 str += ' df: ' + df + '\n'; 300 str += '\n'; 301 if ( decision ) { 302 str += 'Test Decision: '; 303 if ( G > gc ) { 304 str += 'Reject null in favor of alternative at ' + (opts.alpha*100.0) + '% significance level'; 305 } else { 306 str += 'Fail to reject null in favor of alternative at ' + (opts.alpha*100.0) + '% significance level'; 307 } 308 str += '\n'; 309 } 310 return str; 311 } 312 } 313 314 315 // EXPORTS // 316 317 module.exports = incrmgrubbs;