time-to-botec

Benchmark sampling in different programming languages
Log | Files | Refs | README

anova1.js (4834B)


      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 isNumberArray = require( '@stdlib/assert/is-number-array' ).primitives;
     24 var isTypedArrayLike = require( '@stdlib/assert/is-typed-array-like' );
     25 var isArray = require( '@stdlib/assert/is-array' );
     26 var setReadOnly = require( '@stdlib/utils/define-read-only-property' );
     27 var hasOwnProp = require( '@stdlib/assert/has-own-property' );
     28 var cdf = require( './../../base/dists/f/cdf' );
     29 var copy = require( '@stdlib/utils/copy' );
     30 var defaults = require( './defaults.json' );
     31 var validate = require( './validate.js' );
     32 var unique = require( './unique.js' );
     33 var meanTable = require( './mean_table.js' );
     34 var mean = require( './mean.js' );
     35 var prettyPrint = require( './print.js' );
     36 
     37 
     38 // MAIN //
     39 
     40 /**
     41 * Perform a one-way analysis of variance (ANOVA).
     42 *
     43 * @param {NumericArray} x - measured values
     44 * @param {Array} factor - array of treatments
     45 * @param {Options} [options] - function options
     46 * @param {number} [options.alpha=0.05] - significance level
     47 * @throws {TypeError} options argument must be an object
     48 * @throws {TypeError} must provide valid options
     49 * @throws {TypeError} `x` must be a numeric array
     50 * @throws {TypeError} `factor` must be an array
     51 * @throws {RangeError} `factor` must have at least two unique elements
     52 * @throws {RangeError} length of `x` must be greater than or equal to two
     53 * @throws {RangeError} `x` and `factor` must have the same length
     54 * @returns {Object} test results
     55 */
     56 function anova1( x, factor, options ) {
     57 	var meanSumSqTreat; // Mean sum of squares
     58 	var meanSumSqError;
     59 	var ssTreatment;
     60 	var sumSqTotal;
     61 	var sumSqError;
     62 	var treatment; // Index variable
     63 	var grandMean;
     64 	var nGroups;
     65 	var fScore;
     66 	var treats;
     67 	var means;
     68 	var numDf;
     69 	var denDf;
     70 	var nobs;
     71 	var pVal;
     72 	var opts;
     73 	var err;
     74 	var out;
     75 	var sq;
     76 	var i;
     77 
     78 	if ( !isTypedArrayLike( x ) && !isNumberArray( x ) ) {
     79 		throw new TypeError( 'invalid argument. First argument must be a numeric array. Value: `' + x + '`.' );
     80 	}
     81 	opts = copy( defaults );
     82 	if ( arguments.length > 2 ) {
     83 		err = validate( opts, options );
     84 		if ( err ) {
     85 			throw err;
     86 		}
     87 	}
     88 	nobs = x.length;
     89 	if ( nobs <= 1 ) {
     90 		throw new RangeError( 'invalid argument. First argument must have at least two elements. Value: `' + x + '`.' );
     91 	}
     92 	if ( !isArray( factor ) ) {
     93 		throw new TypeError( 'invalid argument. Second argument must be an array. Value: `' + treats + '`.' );
     94 	}
     95 
     96 	treats = unique( factor );
     97 	nGroups = treats.length;
     98 	if ( nGroups <= 1 ) {
     99 		throw new RangeError( 'invalid argument. Second argument must contain at least two unique elements. Value: `' + treats + '`.' );
    100 	}
    101 	if ( nobs !== factor.length ) {
    102 		throw new RangeError( 'invalid arguments. Arguments `x` and `factor` must be arrays of the same length.' );
    103 	}
    104 
    105 	sumSqTotal = 0.0;
    106 	ssTreatment = 0.0;
    107 	means = meanTable( x, factor, treats );
    108 	grandMean = mean( x );
    109 
    110 	// Now get total ss:
    111 	for ( i = 0; i < nobs; i++ ) {
    112 		sq = ( x[i] - grandMean ) * ( x[i] - grandMean );
    113 		sumSqTotal += sq;
    114 	}
    115 
    116 	sq = 0.0;
    117 	for ( treatment in means ) {
    118 		if ( hasOwnProp( means, treatment ) ) {
    119 			// Already have sq defined
    120 			sq = ( means[treatment].mean - grandMean ) *
    121 				( means[treatment].mean - grandMean );
    122 			ssTreatment += means[treatment].sampleSize * sq;
    123 		}
    124 	}
    125 	numDf = nGroups - 1;
    126 	denDf = nobs - nGroups;
    127 	sumSqError = sumSqTotal - ssTreatment;
    128 	meanSumSqTreat = ssTreatment / numDf;
    129 	meanSumSqError = sumSqError / denDf;
    130 	fScore = meanSumSqTreat / meanSumSqError;
    131 
    132 	pVal = 1.0 - cdf( fScore, numDf, denDf );
    133 
    134 	out = {};
    135 
    136 	treatment = {};
    137 	setReadOnly( treatment, 'df', numDf );
    138 	setReadOnly( treatment, 'ss', ssTreatment );
    139 	setReadOnly( treatment, 'ms', meanSumSqTreat );
    140 	setReadOnly( out, 'treatment', treatment );
    141 
    142 	err = {};
    143 	setReadOnly( err, 'df', denDf );
    144 	setReadOnly( err, 'ss', sumSqError );
    145 	setReadOnly( err, 'ms', meanSumSqError );
    146 	setReadOnly( out, 'error', err );
    147 
    148 	setReadOnly( out, 'statistic', fScore );
    149 	setReadOnly( out, 'pValue', pVal );
    150 	setReadOnly( out, 'means', means );
    151 	setReadOnly( out, 'method', 'One-Way ANOVA' );
    152 	setReadOnly( out, 'alpha', opts.alpha );
    153 	setReadOnly( out, 'rejected', pVal <= opts.alpha );
    154 	setReadOnly( out, 'print', prettyPrint( out ) );
    155 	return out;
    156 }
    157 
    158 
    159 // EXPORTS //
    160 
    161 module.exports = anova1;