time-to-botec

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

quantile.js (3284B)

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