temme3.js (6127B)
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/detail/ibeta_inverse.hpp}. The implementation has been modified for JavaScript. 22 * 23 * ```text 24 * Copyright John Maddock 2006. 25 * Copyright Paul A. Bristow 2007. 26 * 27 * Use, modification and distribution are subject to the 28 * Boost Software License, Version 1.0. (See accompanying file 29 * LICENSE or copy at http://www.boost.org/LICENSE_1_0.txt) 30 * ``` 31 */ 32 33 /* eslint-disable no-mixed-operators, max-len */ 34 35 'use strict'; 36 37 // MODULES // 38 39 var gammaincinv = require( './../../../../base/special/gammaincinv' ); 40 var ln = require( './../../../../base/special/ln' ); 41 var sqrt = require( './../../../../base/special/sqrt' ); 42 var SMALLEST_SUBNORMAL = require( '@stdlib/constants/float64/smallest-subnormal' ); 43 var temmeRootFinder = require( './root_finder.js' ); 44 var newtonRaphsonIterate = require( './newton_raphson.js' ); 45 46 47 // MAIN // 48 49 /** 50 * Carries out the third method by Temme (described in section 4). 51 * 52 * ## References 53 * 54 * - Temme, N. M. 1992. "Incomplete Laplace Integrals: Uniform Asymptotic Expansion with Application to the Incomplete Beta Function." _Journal of Computational and Applied Mathematics_ 41 (1–2): 1638–63. doi:[10.1016/0377-0427(92)90244-R](https://doi.org/10.1016/0377-0427(92)90244-R). 55 * 56 * @private 57 * @param {PositiveNumber} a - function parameter 58 * @param {PositiveNumber} b - function parameter 59 * @param {Probability} p - function parameter 60 * @param {Probability} q - probability equal to `1-p` 61 * @returns {number} function value 62 */ 63 function temme3( a, b, p, q ) { 64 var cross; 65 var roots; 66 var lower; 67 var upper; 68 var eta0; 69 var eta; 70 var w10; 71 var w12; 72 var w13; 73 var w14; 74 var e1; 75 var e2; 76 var e3; 77 var mu; 78 var d2; 79 var d3; 80 var d4; 81 var w2; 82 var w3; 83 var w4; 84 var w5; 85 var w6; 86 var w7; 87 var w8; 88 var w9; 89 var w1; 90 var d; 91 var w; 92 var u; 93 var x; 94 95 // Begin by getting an initial approximation for the quantity eta from the dominant part of the incomplete beta: 96 if ( p < q ) { 97 eta0 = gammaincinv( p, b, true ); 98 } else { 99 eta0 = gammaincinv( q, b, false ); 100 } 101 eta0 /= a; 102 103 // Define the variables and powers we'll need later on: 104 mu = b / a; 105 w = sqrt( 1.0+mu ); 106 w2 = w * w; 107 w3 = w2 * w; 108 w4 = w2 * w2; 109 w5 = w3 * w2; 110 w6 = w3 * w3; 111 w7 = w4 * w3; 112 w8 = w4 * w4; 113 w9 = w5 * w4; 114 w10 = w5 * w5; 115 d = eta0 - mu; 116 d2 = d * d; 117 d3 = d2 * d; 118 d4 = d2 * d2; 119 w1 = w + 1.0; 120 w12 = w1 * w1; 121 w13 = w1 * w12; 122 w14 = w12 * w12; 123 124 // Now we need to compute the perturbation error terms that convert eta0 to eta, these are all polynomials of polynomials. Probably these should be re-written to use tabulated data (see examples above), but it's less of a win in this case as we need to calculate the individual powers for the denominator terms anyway, so we might as well use them for the numerator-polynomials as well. Refer to p154-p155 for the details of these expansions: 125 e1 = (w+2.0) * (w-1.0) / (3.0*w); 126 e1 += (w3 + 9.0*w2 + 21.0*w + 5.0) * d / (36.0*w2*w1); 127 e1 -= (w4 - 13.0*w3 + 69.0*w2 + 167.0*w + 46.0) * d2 / (1620.0*w12*w3); 128 e1 -= (7.0*w5 + 21.0*w4 + 70.0*w3 + 26.0*w2 - 93.0*w - 31.0) * d3 / (6480.0*w13*w4); 129 e1 -= (75.0*w6 + 202.0*w5 + 188.0*w4 - 888.0*w3 - 1345.0*w2 + 118.0*w + 138.0) * d4 / (272160.0*w14*w5); 130 131 e2 = (28.0*w4 + 131.0*w3 + 402.0*w2 + 581.0*w + 208.0) * (w-1.0) / (1620.0*w1*w3); 132 e2 -= (35.0*w6 - 154.0*w5 - 623.0*w4 - 1636.0*w3 - 3983.0*w2 - 3514.0*w - 925.0) * d / (12960.0*w12*w4); 133 e2 -= (2132.0*w7 + 7915.0*w6 + 16821.0*w5 + 35066.0*w4 + 87490.0*w3 + 141183.0*w2 + 95993.0*w + 21640.0) * d2 / (816480.0*w5*w13); 134 e2 -= (11053.0*w8 + 53308.0*w7 + 117010.0*w6 + 163924.0*w5 + 116188.0*w4 - 258428.0*w3 - 677042.0*w2 - 481940.0*w - 105497.0) * d3 / (14696640.0*w14*w6); 135 136 e3 = -((3592.0*w7 + 8375.0*w6 - 1323.0*w5 - 29198.0*w4 - 89578.0*w3 - 154413.0*w2 - 116063.0*w - 29632.0) * (w-1.0)) / (816480.0*w5*w12); 137 e3 -= (442043.0*w9 + 2054169.0*w8 + 3803094.0*w7 + 3470754.0*w6 + 2141568.0*w5 - 2393568.0*w4 - 19904934.0*w3 - 34714674.0*w2 - 23128299.0*w - 5253353.0) * d / (146966400.0*w6*w13); 138 e3 -= (116932.0*w10 + 819281.0*w9 + 2378172.0*w8 + 4341330.0*w7 + 6806004.0*w6 + 10622748.0*w5 + 18739500.0*w4 + 30651894.0*w3 + 30869976.0*w2 + 15431867.0*w + 2919016.0) * d2 / (146966400.0*w14*w7); 139 140 // Combine eta0 and the error terms to compute eta (Second equation p155): 141 eta = eta0 + (e1/a) + (e2/(a*a)) + (e3/(a*a*a)); 142 143 /* 144 Now we need to solve Eq 4.2 to obtain x. For any given value of 145 eta there are two solutions to this equation, and since the distribution 146 may be very skewed, these are not related by x ~ 1-x we used when 147 implementing section 3 above. However we know that: 148 149 cross < x <= 1 ; iff eta < mu 150 x == cross ; iff eta == mu 151 0 <= x < cross ; iff eta > mu 152 153 Where cross == 1 / (1 + mu) 154 Many thanks to Prof Temme for clarifying this point. Therefore we'll just jump straight into Newton iterations to solve Eq 4.2 using these bounds, and simple bisection as the first guess, in practice this converges pretty quickly and we only need a few digits correct anyway: 155 */ 156 if ( eta <= 0 ) { 157 eta = SMALLEST_SUBNORMAL; 158 } 159 u = eta - ( mu*ln(eta) ) + ( ( 1.0+mu ) * ln( 1.0+mu ) ) - mu; 160 cross = 1.0 / ( 1.0+mu ); 161 lower = (eta < mu) ? cross : 0.0; 162 upper = (eta < mu) ? 1.0 : cross; 163 x = (lower+upper) / 2.0; 164 roots = temmeRootFinder( u, mu ); 165 return newtonRaphsonIterate( roots, x, lower, upper, 32, 100 ); 166 } 167 168 169 // EXPORTS // 170 171 module.exports = temme3;