time-to-botec

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

main.js (8201B)


      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 isNonNegativeInteger = require( '@stdlib/assert/is-nonnegative-integer' ).isPrimitive;
     24 var isCollection = require( '@stdlib/assert/is-collection' );
     25 var isndarrayLike = require( '@stdlib/assert/is-ndarray-like' );
     26 var isNumber = require( '@stdlib/assert/is-number' ).isPrimitive;
     27 var isString = require( '@stdlib/assert/is-string' ).isPrimitive;
     28 var absdiff = require( '@stdlib/math/base/utils/absolute-difference' );
     29 var FLOAT64_SQRT_EPS = require( '@stdlib/constants/float64/sqrt-eps' );
     30 var PINF = require( '@stdlib/constants/float64/pinf' );
     31 var chisqCDF = require( './../../base/dists/chisquare/cdf' );
     32 var isnan = require( '@stdlib/assert/is-nan' );
     33 var daxpy = require( '@stdlib/blas/base/daxpy' );
     34 var dscal = require( '@stdlib/blas/base/dscal' );
     35 var dsumpw = require( '@stdlib/blas/ext/base/dsumpw' );
     36 var Float64Array = require( '@stdlib/array/float64' );
     37 var defaults = require( './defaults.js' );
     38 var validate = require( './validate.js' );
     39 var getPMF = require( './get_pmf.js' );
     40 var testStatistic = require( './statistic.js' );
     41 var simulate = require( './simulate.js' );
     42 var Results = require( './results.js' );
     43 
     44 
     45 // MAIN //
     46 
     47 /**
     48 * Performs a chi-square goodness-of-fit test.
     49 *
     50 * @param {(Collection|VectorLike)} x - observation frequencies
     51 * @param {(Collection|VectorLike|string)} y - expected frequencies or a discrete probability distribution name
     52 * @param {...number} [args] - probability mass function (PMF) arguments
     53 * @param {Options} [options] - function options
     54 * @param {number} [options.alpha=0.05] - significance level
     55 * @param {NonNegativeInteger} [options.ddof=0] - degrees of freedom adjustment
     56 * @param {boolean} [options.simulate=false] - boolean indicating whether to compute p-values by Monte Carlo simulation
     57 * @param {PositiveInteger} [options.iterations=500] - number of Monte Carlo iterations
     58 * @throws {TypeError} first argument must be an array-like object or a 1-dimensional array containing nonnegative integers
     59 * @throws {TypeError} second argument must be either an array-like object (or a 1-dimensional array) of nonnegative numbers, an array-like object (or a 1-dimensional array) of probabilities summing to one, or a discrete probability distribution name
     60 * @throws {TypeError} options argument must be an object
     61 * @throws {TypeError} must provide valid options
     62 * @throws {Error} first and second arguments must have the same length
     63 * @throws {Error} first argument must contain at least one element greater than zero
     64 * @throws {RangeError} significance level must be a number on the interval `[0,1]`
     65 * @throws {TypeError} probability mass function (PMF) arguments must be number primitives
     66 * @returns {Object} test results
     67 *
     68 * @example
     69 * var x = [ 89, 37, 30, 28, 2 ];
     70 * var p = [ 0.40, 0.20, 0.20, 0.15, 0.05 ];
     71 *
     72 * var out = chi2gof( x, p );
     73 *
     74 * var o = out.toJSON();
     75 * // returns { 'pValue': ~0.0406, 'statistic': ~9.9901, ... }
     76 */
     77 function chi2gof( x, y ) {
     78 	var expected;
     79 	var nargs;
     80 	var args;
     81 	var opts;
     82 	var pval;
     83 	var stat;
     84 	var obs;
     85 	var err;
     86 	var pmf;
     87 	var sum;
     88 	var df;
     89 	var N;
     90 	var d;
     91 	var s;
     92 	var o;
     93 	var n;
     94 	var p;
     95 	var v;
     96 	var i;
     97 
     98 	if ( isndarrayLike( x ) && x.ndims === 1 && x.strides.length === 1 ) { // is ndarray-like vector?
     99 		d = x.data;
    100 		s = x.strides[ 0 ];
    101 		o = x.offset;
    102 	} else if ( isCollection( x ) ) {
    103 		d = x;
    104 		s = 1;
    105 		o = 0;
    106 	} else {
    107 		throw new TypeError( 'invalid argument. First argument must be either an array-like object or a 1-dimensional ndarray. Value: `' + x + '`.' );
    108 	}
    109 	N = x.length;
    110 
    111 	// Initialize an array for storing a copy of the observations array:
    112 	obs = new Float64Array( N+1 ); // Note: `N+1` is intentional in the event that we need to add a remaining category for all values greater than or equal to `N`
    113 
    114 	n = 0;
    115 	for ( i = 0; i < N; i++ ) {
    116 		v = d[ o+(s*i) ];
    117 		if ( !isNonNegativeInteger( v ) ) {
    118 			throw new TypeError( 'invalid argument. First argument must contain nonnegative integers. Index: `' + i + '`. Value: `' + v + '`.' );
    119 		}
    120 		obs[ i ] = v;
    121 		n += v;
    122 	}
    123 	if ( n === 0 ) {
    124 		throw new Error( 'invalid argument. First argument must contain at least one element greater than zero (i.e., the total number number of observations must be greater than zero).' );
    125 	}
    126 	// NOTE: `obs` is now a single-segment contiguous Float64Array
    127 
    128 	nargs = 0;
    129 	if ( isString( y ) ) {
    130 		pmf = getPMF( y );
    131 		if ( pmf instanceof Error ) {
    132 			throw pmf;
    133 		}
    134 		nargs += pmf.length - 1; // WARNING: this relies on PMF functions having an explicit arity
    135 		args = [ 0 ];
    136 		for ( i = 0; i < nargs; i++ ) {
    137 			v = arguments[ i+2 ];
    138 			if ( !isNumber( v ) || isnan( v ) ) {
    139 				throw new TypeError( 'invalid argument. Probability mass function (PMF) arguments must be number primitives. Argument: `' + (i+2) + '`. Value: `' + v + '`.' );
    140 			}
    141 			args.push( v );
    142 		}
    143 		expected = new Float64Array( N+1 );
    144 		sum = 0.0;
    145 		for ( i = 0; i < N; i++ ) {
    146 			args[ 0 ] = i;
    147 			if ( y === 'discrete-uniform' ) {
    148 				args[ 0 ] += args[ 1 ]; // scales the value at which to evaluate the PMF based on the minimum support of the distribution (which should have been provided as the first distribution parameter)
    149 			}
    150 			v = pmf.apply( null, args );
    151 			sum += v;
    152 			expected[ i ] = v * n;
    153 		}
    154 		// Check whether we need to add a remaining category for all values greater than or equal to `N`...
    155 		if ( sum < 1.0 ) {
    156 			expected[ N ] = (1.0-sum) * n;
    157 			N += 1;
    158 		}
    159 	} else {
    160 		if ( isndarrayLike( y ) && y.ndims === 1 && y.strides.length === 1 ) { // is ndarray-like vector?
    161 			d = y.data;
    162 			s = y.strides[ 0 ];
    163 			o = y.offset;
    164 		} else if ( isCollection( y ) ) {
    165 			d = y;
    166 			s = 1;
    167 			o = 0;
    168 		} else {
    169 			throw new TypeError( 'invalid argument. Second argument must be either an array-like object (or 1-dimensional ndarray) of probabilities summing to one, an array-like object (or 1-dimensional ndarray) of expected frequencies, or a discrete probability distribution name. Value: `' + y + '`.' );
    170 		}
    171 		if ( y.length !== N ) {
    172 			throw new Error( 'invalid arguments. First and second arguments must have the same length.' );
    173 		}
    174 		expected = new Float64Array( N );
    175 		sum = 0.0;
    176 		for ( i = 0; i < N; i++ ) {
    177 			v = d[ o+(s*i) ];
    178 			if ( !isNumber( v ) ) {
    179 				throw new TypeError( 'invalid argument. Second argument must only contain numbers. Index: `' + i + '`. Value: `' + v + '`.' );
    180 			}
    181 			if ( v < 0.0 ) {
    182 				throw new TypeError( 'invalid argument. Second argument must only contain nonnegative numbers. Index: `' + i + '`. Value: `' + v + '`.' );
    183 			} else if ( v > 1.0 ) {
    184 				sum += PINF;
    185 			} else {
    186 				sum += v;
    187 			}
    188 			expected[ i ] = v;
    189 		}
    190 		// Check if provided a unity probability array (otherwise, assume provided an expected frequencies array)...
    191 		if ( absdiff( sum, 1.0 ) <= FLOAT64_SQRT_EPS ) {
    192 			p = y; // NOTE: `y` may not be a Float64Array
    193 			expected = dscal( N, n, expected, 1 );
    194 		}
    195 	}
    196 	// NOTE: `expected` is now a single-segment contiguous Float64Array
    197 
    198 	opts = defaults();
    199 	if ( arguments.length > 2+nargs ) {
    200 		err = validate( opts, arguments[ 2+nargs ] );
    201 		if ( err ) {
    202 			throw err;
    203 		}
    204 	}
    205 	stat = testStatistic( N, obs, 1, expected, 1 ); // TODO: consider replacing with low-level double-precision strided interface
    206 	if ( opts.simulate ) {
    207 		if ( p === void 0 ) {
    208 			v = dsumpw( N, expected, 1 );
    209 			p = daxpy( N, 1.0/v, expected, 1, new Float64Array( N ), 1 );
    210 		}
    211 		pval = simulate( N, expected, p, stat, n, opts.iterations );
    212 	} else {
    213 		df = N - 1 - opts.ddof;
    214 		pval = 1.0 - chisqCDF( stat, df );
    215 	}
    216 	return new Results( pval, opts.alpha, stat, ( df === void 0 ) ? null : df );
    217 }
    218 
    219 
    220 // EXPORTS //
    221 
    222 module.exports = chi2gof;