time-to-botec

inverse_students_t.js (5918B)

---

 /**
 * @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/detail/t_distribution_inv.hpp}. The implementation has been modified for JavaScript.
 *
 * ```text
 * Copyright John Maddock 2006.
 * Copyright Paul A. Bristow 2007.
 *
 * 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 erfcinv = require( './../../../../base/special/erfcinv' );
 var floor = require( './../../../../base/special/floor' );
 var ldexp = require( './../../../../base/special/ldexp' );
 var round = require( './../../../../base/special/round' );
 var acos = require( './../../../../base/special/acos' );
 var sqrt = require( './../../../../base/special/sqrt' );
 var abs = require( './../../../../base/special/abs' );
 var cos = require( './../../../../base/special/cos' );
 var pow = require( './../../../../base/special/pow' );
 var sin = require( './../../../../base/special/sin' );
 var SQRT2 = require( '@stdlib/constants/float64/sqrt-two' );
 var PI = require( '@stdlib/constants/float64/pi' );
 var inverseStudentsTBodySeries = require( './inverse_students_t_body_series.js' );
 var inverseStudentsTTailSeries = require( './inverse_students_t_tail_series.js' );
 var inverseStudentsTHill = require( './inverse_students_t_hill.js' );
 
 
 // VARIABLES //
 
 var DF_THRESHOLD = 0x10000000; // 2^28
 var ONE_THIRD = 1.0 / 3.0;
 var EXP = ( 2.0 * 53.0 ) / 3.0;
 var C = 0.85498797333834849467655443627193;
 
 
 // MAIN //
 
 /**
 * Evaluates Student's t quantiles.
 *
 * @private
 * @param {PositiveNumber} df - degrees of freedom
 * @param {Probability} u - input probability
 * @param {Probability} v - probability equal to `1-u`
 * @returns {number} function value
 */
 function inverseStudentsT( df, u, v ) {
 	var crossover;
 	var tolerance;
 	var rootAlpha;
 	var invert;
 	var result;
 	var alpha;
 	var tmp;
 	var p0;
 	var p2;
 	var p4;
 	var p5;
 	var p;
 	var r;
 	var x;
 	var a;
 	var b;
 
 	result = 0;
 	if ( u > v ) {
 		// Function is symmetric, so invert it:
 		tmp = v;
 		v = u;
 		u = tmp;
 		invert = true;
 	} else {
 		invert = false;
 	}
 	if ( floor(df) === df && df < 20 ) {
 		// We have integer degrees of freedom, try for the special cases first:
 		tolerance = ldexp( 1.0, EXP );
 
 		switch ( floor( df ) ) {
 		case 1:
 			// `df = 1` is the same as the Cauchy distribution, see Shaw Eq 35:
 			if ( u === 0.5 ) {
 				result = 0.0;
 			} else {
 				result = -cos( PI * u ) / sin( PI * u );
 			}
 			break;
 		case 2:
 			// `df = 2` has an exact result, see Shaw Eq 36:
 			result = ( (2.0*u) - 1.0 ) / sqrt( 2.0 * u * v );
 			break;
 		case 4:
 			// `df = 4` has an exact result, see Shaw Eq 38 & 39:
 			alpha = 4.0 * u * v;
 			rootAlpha = sqrt( alpha );
 			r = 4 * cos( acos( rootAlpha ) / 3.0 ) / rootAlpha;
 			x = sqrt( r - 4.0 );
 			result = ( u - 0.5 < 0.0 ) ? -x : x;
 			break;
 		case 6:
 			// We get numeric overflow in this area:
 			if ( u < 1.0e-150 ) {
 				return ( ( invert ) ? -1 : 1 ) * inverseStudentsTHill( df, u );
 			}
 			// Newton-Raphson iteration of a polynomial case, choice of seed value is taken from Shaw's online supplement:
 			a = 4.0 * ( u - (u*u) );// 1 - 4 * (u - 0.5f) * (u - 0.5f);
 			b = pow( a, ONE_THIRD );
 			p = 6.0 * ( 1.0 + ( C * ( (1.0/b) - 1.0 ) ) );
 			do {
 				p2 = p * p;
 				p4 = p2 * p2;
 				p5 = p * p4;
 				p0 = p;
 
 				// Next term is given by Eq 41:
 				p = 2.0 * ( (8.0*a*p5) - (270.0*p2) + 2187 ) /
 					( 5.0 * ( (4.0*a*p4) - (216.0*p) - 243.0 ) );
 			} while ( abs( (p - p0) / p ) > tolerance );
 
 			// Use Eq 45 to extract the result:
 			p = sqrt( p - df );
 			result = ( u - 0.5 < 0.0 ) ? -p : p;
 			break;
 		default:
 			if ( df > DF_THRESHOLD ) { // 2^28
 				result = erfcinv( 2.0 * u ) * SQRT2;
 			} else if ( df < 3 ) {
 				// Use a roughly linear scheme to choose between Shaw's tail series and body series:
 				crossover = 0.2742 - ( df * 0.0242143 );
 				if ( u > crossover ) {
 					result = inverseStudentsTBodySeries( df, u );
 				} else {
 					result = inverseStudentsTTailSeries( df, u );
 				}
 			} else {
 				// Use Hill's method except in the extreme tails where we use Shaw's tail series. The crossover point is roughly exponential in -df:
 				crossover = ldexp( 1.0, round( df / -0.654 ) );
 				if ( u > crossover ) {
 					result = inverseStudentsTHill( df, u );
 				} else {
 					result = inverseStudentsTTailSeries( df, u );
 				}
 			}
 		}
 	} else if ( df > DF_THRESHOLD ) {
 		result = -erfcinv( 2.0 * u ) * SQRT2;
 	} else if ( df < 3 ) {
 		// Use a roughly linear scheme to choose between Shaw's tail series and body series:
 		crossover = 0.2742 - ( df * 0.0242143 );
 		if ( u > crossover ) {
 			result = inverseStudentsTBodySeries( df, u );
 		} else {
 			result = inverseStudentsTTailSeries( df, u );
 		}
 	} else {
 		// Use Hill's method except in the extreme tails where we use Shaw's tail series. The crossover point is roughly exponential in -df:
 		crossover = ldexp( 1.0, round( df / -0.654 ) );
 		if ( u > crossover ) {
 			result = inverseStudentsTHill( df, u );
 		} else {
 			result = inverseStudentsTTailSeries( df, u );
 		}
 	}
 	return ( invert ) ? -result : result;
 }
 
 
 // EXPORTS //
 
 module.exports = inverseStudentsT;