time-to-botec

ibeta_power_terms.js (7497B)

---

 /**
 * @license Apache-2.0
 *
 * Copyright (c) 2018 The Stdlib Authors.
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *    http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *
 *
 * ## Notice
 *
 * The original C++ code and copyright notice are from the [Boost library]{@link http://www.boost.org/doc/libs/1_62_0/boost/math/special_functions/beta.hpp}. The implementation has been modified for JavaScript.
 *
 * ```text
 * (C) Copyright John Maddock 2006.
 *
 * Use, modification and distribution are subject to the
 * Boost Software License, Version 1.0. (See accompanying file
 * LICENSE or copy at http://www.boost.org/LICENSE_1_0.txt)
 * ```
 */
 
 'use strict';
 
 // MODULES //
 
 var lanczosSumExpGScaled = require( './../../../../base/special/gamma-lanczos-sum-expg-scaled' );
 var maxabs = require( './../../../../base/special/maxabs' );
 var minabs = require( './../../../../base/special/minabs' );
 var expm1 = require( './../../../../base/special/expm1' );
 var log1p = require( './../../../../base/special/log1p' );
 var sqrt = require( './../../../../base/special/sqrt' );
 var abs = require( './../../../../base/special/abs' );
 var exp = require( './../../../../base/special/exp' );
 var pow = require( './../../../../base/special/pow' );
 var min = require( './../../../../base/special/min' );
 var ln = require( './../../../../base/special/ln' );
 var MAX_LN = require( '@stdlib/constants/float64/max-ln' );
 var MIN_LN = require( '@stdlib/constants/float64/min-ln' );
 var G = require( '@stdlib/constants/float64/gamma-lanczos-g' );
 var E = require( '@stdlib/constants/float64/e' );
 
 
 // MAIN //
 
 /**
 * Computes the leading power terms in the incomplete beta function.
 *
 * When normalized,
 *
 * ```tex
 * \frac{ x^a y^b }{ \operatorname{Beta}(a,b) }
 * ```
 *
 * and otherwise
 *
 * ```tex
 * x^a y^b
 * ```
 *
 * ## Notes
 *
 * -   Almost all of the error in the incomplete beta comes from this function, particularly when \\( a \\) and \\( b \\) are large. Computing large powers are _hard_ though, and using logarithms just leads to horrendous cancellation errors.
 *
 * -   For \\( l1 * l2 > 0 \\) or \\( \operatorname{min}( a, b ) < 1 \\), the two power terms both go in the same direction (toward zero or toward infinity). In this case if either term overflows or underflows, then the product of the two must do so also. Alternatively, if one exponent is less than one, then we can't productively use it to eliminate overflow or underflow from the other term.  Problems with spurious overflow/underflow can't be ruled out. In this case, but it is _very_ unlikely since one of the power terms will evaluate to a number close to 1.
 *
 * -   If \\( \max( \abs(l1), \abs(l2) ) < 0.5 \\), both exponents are near one and both the exponents are greater than one, and, further, these two power terms tend in opposite directions (one toward zero, the other toward infinity), so we have to combine the terms to avoid any risk of overflow or underflow. We do this by moving one power term inside the other, we have:
 *
 *     ```tex
 *     (1 + l_1)^a \cdot (1 + l_2)^b \\
 *     = ((1 + l_1) \cdot (1 + l_2)^(b/a))^a \\
 *     = (1 + l_1 + l_3 + l_1*l_3)^a
 *     ```
 *
 *     and
 *
 *     ```tex
 *     l_3 = (1 + l_2)^(b/a) - 1 \\
 *     = \exp((b/a) * \ln(1 + l_2)) - 1
 *     ```
 *
 *     The tricky bit is deciding which term to move inside. By preference, we move the larger term inside, so that the size of the largest exponent is reduced.  However, that can only be done as long as l3 (see above) is also small.
 *
 * @private
 * @param {NonNegativeNumber} a - function parameter
 * @param {NonNegativeNumber} b - function parameter
 * @param {Probability} x - function parameter
 * @param {Probability} y - probability equal to `1-x`
 * @param {boolean} normalized - boolean indicating whether to evaluate the power terms of the regularized or non-regularized incomplete beta function
 * @returns {number} power terms
 */
 function ibetaPowerTerms( a, b, x, y, normalized ) {
 	var result;
 	var smallA;
 	var ratio;
 	var agh;
 	var bgh;
 	var cgh;
 	var l1;
 	var l2;
 	var l3;
 	var p1;
 	var b1;
 	var b2;
 	var c;
 	var l;
 
 	if ( !normalized ) {
 		// Can we do better here?
 		return pow( x, a ) * pow( y, b );
 	}
 	c = a + b;
 
 	// Combine power terms with Lanczos approximation:
 	agh = a + G - 0.5;
 	bgh = b + G - 0.5;
 	cgh = c + G - 0.5;
 	result = lanczosSumExpGScaled( c );
 	result /= lanczosSumExpGScaled( a ) * lanczosSumExpGScaled( b );
 
 	// Combine with the leftover terms from the Lanczos approximation:
 	result *= sqrt( bgh / E );
 	result *= sqrt( agh / cgh );
 
 	// `l1` and `l2` are the base of the exponents minus one:
 	l1 = ( ( x * b ) - ( y * agh ) ) / agh;
 	l2 = ( ( y * a ) - ( x * bgh ) ) / bgh;
 	if ( minabs( l1, l2 ) < 0.2 ) {
 		// When the base of the exponent is very near 1 we get really gross errors unless extra care is taken:
 		if ( l1 * l2 > 0 || min( a, b ) < 1 ) {
 			if ( abs(l1) < 0.1 ) {
 				result *= exp( a * log1p( l1 ) );
 			} else {
 				result *= pow( ( x*cgh ) / agh, a );
 			}
 			if ( abs(l2) < 0.1 ) {
 				result *= exp( b * log1p( l2 ) );
 			} else {
 				result *= pow((y * cgh) / bgh, b);
 			}
 		}
 		else if ( maxabs( l1, l2 ) < 0.5 ) {
 			smallA = a < b;
 			ratio = b / a;
 			if (
 				(smallA && (ratio * l2 < 0.1)) ||
 				(!smallA && (l1 / ratio > 0.1))
 			) {
 				l3 = expm1( ratio * log1p( l2 ) );
 				l3 = l1 + l3 + ( l3 * l1 );
 				l3 = a * log1p( l3 );
 				result *= exp( l3 );
 			}
 			else {
 				l3 = expm1( log1p( l1 ) / ratio );
 				l3 = l2 + l3 + ( l3 * l2 );
 				l3 = b * log1p( l3 );
 				result *= exp( l3 );
 			}
 		}
 		else if ( abs(l1) < abs(l2) ) {
 			// First base near 1 only:
 			l = ( a * log1p( l1 ) ) + ( b * ln( ( y*cgh ) / bgh ) );
 			if ( l <= MIN_LN || l >= MAX_LN ) {
 				l += ln(result);
 				if ( l >= MAX_LN ) {
 					return NaN;
 				}
 				result = exp( l );
 			} else {
 				result *= exp( l );
 			}
 		}
 		else {
 			// Second base near 1 only:
 			l = ( b * log1p( l2 ) ) + ( a * ln( (x*cgh) / agh ) );
 			if ( l <= MIN_LN || l >= MAX_LN ) {
 				l += ln(result);
 				if ( l >= MAX_LN ) {
 					return NaN;
 				}
 				result = exp( l );
 			} else {
 				result *= exp( l );
 			}
 		}
 	}
 	else {
 		// General case:
 		b1 = (x * cgh) / agh;
 		b2 = (y * cgh) / bgh;
 		l1 = a * ln(b1);
 		l2 = b * ln(b2);
 		if (
 			l1 >= MAX_LN ||
 			l1 <= MIN_LN ||
 			l2 >= MAX_LN ||
 			l2 <= MIN_LN
 		) {
 			// Oops, under/overflow, sidestep if we can:
 			if ( a < b ) {
 				p1 = pow( b2, b / a );
 				l3 = a * ( ln(b1) + ln(p1) );
 				if ( l3 < MAX_LN && l3 > MIN_LN ) {
 					result *= pow( p1 * b1, a );
 				} else {
 					l2 += l1 + ln(result);
 					if ( l2 >= MAX_LN ) {
 						return NaN;
 					}
 					result = exp( l2 );
 				}
 			}
 			else {
 				p1 = pow( b1, a / b );
 				l3 = ( ln(p1) + ln(b2) ) * b;
 				if ( l3 < MAX_LN && l3 > MIN_LN ) {
 					result *= pow( p1 * b2, b );
 				} else {
 					l2 += l1 + ln( result );
 					if (l2 >= MAX_LN) {
 						return NaN;
 					}
 					result = exp( l2 );
 				}
 			}
 		}
 		else {
 			// Finally the normal case:
 			result *= pow( b1, a ) * pow( b2, b );
 		}
 	}
 	return result;
 }
 
 
 // EXPORTS //
 
 module.exports = ibetaPowerTerms;