time-to-botec

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

factory.js (3165B)

      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/math/base/assert/is-nonnegative-integer' );
     24 var constantFunction = require( '@stdlib/utils/constant-function' );
     25 var degenerate = require( './../../../../../base/dists/degenerate/quantile' ).factory;
     26 var erfcinv = require( '@stdlib/math/base/special/erfcinv' );
     27 var isnan = require( '@stdlib/math/base/assert/is-nan' );
     28 var round = require( '@stdlib/math/base/special/round' );
     29 var sqrt = require( '@stdlib/math/base/special/sqrt' );
     30 var cdf = require( './../../../../../base/dists/binomial/cdf' );
     31 var SQRT2 = require( '@stdlib/constants/float64/sqrt-two' );
     32 var PINF = require( '@stdlib/constants/float64/pinf' );
     33 var searchLeft = require( './search_left.js' );
     34 var searchRight = require( './search_right.js' );
     35 
     36 
     37 // MAIN //
     38 
     39 /**
     40 * Returns a function for evaluating the quantile function for a binomial distribution with number of trials `n` and success probability `p`.
     41 *
     42 * @param {NonNegativeInteger} n - number of trials
     43 * @param {Probability} p - success probability
     44 * @returns {Function} quantile function
     45 *
     46 * @example
     47 * var quantile = factory( 10, 0.5 );
     48 * var y = quantile( 0.1 );
     49 * // returns 3
     50 *
     51 * y = quantile( 0.9 );
     52 * // returns 7
     53 */
     54 function factory( n, p ) {
     55 	var sigmaInv;
     56 	var sigma;
     57 	var mu;
     58 
     59 	if (
     60 		isnan( n ) ||
     61 		isnan( p ) ||
     62 		!isNonNegativeInteger( n ) ||
     63 		n === PINF ||
     64 		p < 0.0 ||
     65 		p > 1.0
     66 	) {
     67 		return constantFunction( NaN );
     68 	}
     69 	if ( p === 0.0 || n === 0.0 ) {
     70 		return degenerate( 0.0 );
     71 	}
     72 	if ( p === 1.0 ) {
     73 		return degenerate( n );
     74 	}
     75 	mu = n * p;
     76 	sigma = sqrt( n * p * ( 1.0-p ) );
     77 	sigmaInv = 1.0 / sigma;
     78 	return quantile;
     79 
     80 	/**
     81 	* Evaluates the quantile function for a binomial distribution.
     82 	*
     83 	* @private
     84 	* @param {Probability} r - input value
     85 	* @returns {NonNegativeInteger} evaluated quantile function
     86 	*
     87 	* @example
     88 	* var y = quantile( 0.3 );
     89 	* // returns <number>
     90 	*/
     91 	function quantile( r ) {
     92 		var guess;
     93 		var corr;
     94 		var x2;
     95 		var x;
     96 
     97 		if ( isnan( r ) || r < 0.0 || r > 1.0 ) {
     98 			return NaN;
     99 		}
    100 		if ( r === 0.0 ) {
    101 			return 0;
    102 		}
    103 		if ( r === 1.0 ) {
    104 			return n;
    105 		}
    106 		// Cornish-Fisher expansion:
    107 		if ( r < 0.5 ) {
    108 			x = -erfcinv( 2.0 * r ) * SQRT2;
    109 		} else {
    110 			x = erfcinv( 2.0 * ( 1.0-r ) ) * SQRT2;
    111 		}
    112 		x2 = x * x;
    113 
    114 		// Skewness correction:
    115 		corr = x + ( sigmaInv * ( x2-1.0 ) / 6.0 );
    116 		guess = round( mu + (sigma * corr) );
    117 		if ( cdf( guess, n, p ) >= r ) {
    118 			return searchLeft( guess, r, n, p );
    119 		}
    120 		return searchRight( guess, r, n, p );
    121 	}
    122 }
    123 
    124 
    125 // EXPORTS //
    126 
    127 module.exports = factory;