time-to-botec

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

ibeta_series.js (4600B)

      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 * ## Notice
     20 *
     21 * The original C++ code and copyright notice are from the [Boost library]{@link http://www.boost.org/doc/libs/1_64_0/boost/math/special_functions/beta.hpp}. The implementation has been modified for JavaScript.
     22 *
     23 * ```text
     24 * (C) Copyright John Maddock 2006.
     25 *
     26 * Use, modification and distribution are subject to the
     27 * Boost Software License, Version 1.0. (See accompanying file
     28 * LICENSE or copy at http://www.boost.org/LICENSE_1_0.txt)
     29 * ```
     30 */
     31 
     32 'use strict';
     33 
     34 // MODULES //
     35 
     36 var lanczosSumExpGScaled = require( './../../../../base/special/gamma-lanczos-sum-expg-scaled' );
     37 var sumSeries = require( './../../../../base/tools/sum-series' );
     38 var log1p = require( './../../../../base/special/log1p' );
     39 var sqrt = require( './../../../../base/special/sqrt' );
     40 var exp = require( './../../../../base/special/exp' );
     41 var pow = require( './../../../../base/special/pow' );
     42 var ln = require( './../../../../base/special/ln' );
     43 var MIN_VALUE = require( '@stdlib/constants/float64/smallest-normal' );
     44 var MAX_LN = require( '@stdlib/constants/float64/max-ln' );
     45 var MIN_LN = require( '@stdlib/constants/float64/min-ln' );
     46 var G = require( '@stdlib/constants/float64/gamma-lanczos-g' );
     47 var E = require( '@stdlib/constants/float64/e' );
     48 
     49 
     50 // VARIABLES //
     51 
     52 var opts = {
     53 	'maxTerms': 100
     54 };
     55 
     56 
     57 // FUNCTIONS //
     58 
     59 /**
     60 * Series approximation to the incomplete beta.
     61 *
     62 * @private
     63 * @param {NonNegativeNumber} a - function parameter
     64 * @param {NonNegativeNumber} b - function parameter
     65 * @param {Probability} x - function parameter
     66 * @param {number} result - initial result value
     67 * @returns {Function} series function
     68 */
     69 function ibetaSeriesT( a, b, x, result ) {
     70 	var poch = 1.0 - b;
     71 	var n = 1;
     72 	return next;
     73 
     74 	/**
     75 	* Calculate the next term of the series.
     76 	*
     77 	* @private
     78 	* @returns {number} series expansion term
     79 	*/
     80 	function next() {
     81 		var r = result / a;
     82 		a += 1.0;
     83 		result *= poch * x / n;
     84 		n += 1;
     85 		poch += 1.0;
     86 		return r;
     87 	}
     88 }
     89 
     90 
     91 // MAIN //
     92 
     93 /**
     94 * Incomplete beta series.
     95 *
     96 * @private
     97 * @param {NonNegativeNumber} a - function parameter
     98 * @param {NonNegativeNumber} b - function parameter
     99 * @param {Probability} x - function parameter
    100 * @param {NonNegativeInteger} s0 - initial value
    101 * @param {boolean} normalized - boolean indicating whether to evaluate the power terms of the regularized or non-regularized incomplete beta function
    102 * @param {(Array|TypedArray|Object)} out - output array holding the derivative as the second element
    103 * @param {Probability} y - probability equal to `1-x`
    104 * @returns {number} function value
    105 */
    106 function ibetaSeries( a, b, x, s0, normalized, out, y ) {
    107 	var result;
    108 	var agh;
    109 	var bgh;
    110 	var cgh;
    111 	var l1;
    112 	var l2;
    113 	var c;
    114 	var s;
    115 
    116 	if ( normalized ) {
    117 		c = a + b;
    118 
    119 		// Incomplete beta power term, combined with the Lanczos approximation:
    120 		agh = a + G - 0.5;
    121 		bgh = b + G - 0.5;
    122 		cgh = c + G - 0.5;
    123 		result = lanczosSumExpGScaled( c ) / ( lanczosSumExpGScaled( a ) * lanczosSumExpGScaled( b ) ); // eslint-disable-line max-len
    124 
    125 		l1 = ln( cgh / bgh ) * ( b - 0.5 );
    126 		l2 = ln( x * cgh / agh ) * a;
    127 
    128 		// Check for over/underflow in the power terms:
    129 		if (
    130 			l1 > MIN_LN &&
    131 			l1 < MAX_LN &&
    132 			l2 > MIN_LN &&
    133 			l2 < MAX_LN
    134 		) {
    135 			if ( a * b < bgh * 10.0 ) {
    136 				result *= exp( ( b-0.5 ) * log1p( a / bgh ) );
    137 			} else {
    138 				result *= pow( cgh / bgh, b - 0.5 );
    139 			}
    140 			result *= pow( x * cgh / agh, a );
    141 			result *= sqrt( agh / E );
    142 
    143 			if ( out ) {
    144 				out[ 1 ] = result * pow( y, b );
    145 			}
    146 		}
    147 		else {
    148 			// We need logs, and this *will* cancel:
    149 			result = ln( result ) + l1 + l2 + ( ( ln( agh ) - 1.0 ) / 2.0 );
    150 			if ( out ) {
    151 				out[ 1 ] = exp( result + ( b * ln( y ) ) );
    152 			}
    153 			result = exp( result );
    154 		}
    155 	}
    156 	else {
    157 		// Non-normalized, just compute the power:
    158 		result = pow( x, a );
    159 	}
    160 	if ( result < MIN_VALUE ) {
    161 		return s0; // Safeguard: series can't cope with denorms.
    162 	}
    163 	s = ibetaSeriesT( a, b, x, result );
    164 	opts.initialValue = s0;
    165 	return sumSeries( s, opts );
    166 }
    167 
    168 
    169 // EXPORTS //
    170 
    171 module.exports = ibetaSeries;