compute.js (6942B)
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 /* eslint-disable max-statements */ 20 21 'use strict'; 22 23 // MODULES // 24 25 var logger = require( 'debug' ); 26 var evalpoly = require( './../../../../base/tools/evalpoly' ); 27 var gammaln = require( './../../../../base/special/gammaln' ); 28 var erfcinv = require( './../../../../base/special/erfcinv' ); 29 var gamma = require( './../../../../base/special/gamma' ); 30 var sqrt = require( './../../../../base/special/sqrt' ); 31 var abs = require( './../../../../base/special/abs' ); 32 var exp = require( './../../../../base/special/exp' ); 33 var min = require( './../../../../base/special/min' ); 34 var pow = require( './../../../../base/special/pow' ); 35 var ln = require( './../../../../base/special/ln' ); 36 var SQRT_TWO_PI = require( '@stdlib/constants/float64/sqrt-two-pi' ); 37 var MAX_FLOAT32 = require( '@stdlib/constants/float32/max' ); 38 var TWO_PI = require( '@stdlib/constants/float64/two-pi' ); 39 var higherNewton = require( './higher_newton.js' ); 40 var lambdaeta = require( './lambdaeta.js' ); 41 var gamstar = require( './gamstar.js' ); 42 var eps1 = require( './eps1.js' ); 43 var eps2 = require( './eps2.js' ); 44 var eps3 = require( './eps3.js' ); 45 46 47 // VARIABLES // 48 49 var debug = logger( 'gammaincinv:compute' ); 50 var HALF = 0.5; 51 var ONEO3 = 0.333333333333333333333333333333; 52 var ONEO4 = 0.25; 53 var ONEO5 = 0.2; 54 var ONEO6 = 0.166666666666666666666666666667; 55 var ONEO12 = 0.0833333333333333333333333333333; 56 var ONEO24 = 0.0416666666666666666666666666667; 57 58 // Coefficient workspace: 59 var CK = [ 0.0, 0.0, 0.0, 0.0, 0.0 ]; // WARNING: not thread safe 60 61 62 // MAIN // 63 64 /** 65 * Computes `x` in the equations `P(a,xr) = p` and `Q(a,xr) = q`, where `a` is a positive parameter and `p` and `q` satisfy `p+q = 1`. 66 * 67 * ## Notes 68 * 69 * - The equation is inverted with `min(p,q)`. 70 * 71 * @private 72 * @param {number} a - scale value of incomplete gamma function 73 * @param {Probability} p - probability value 74 * @param {Probability} q - probability value 75 * @returns {number} solution of the equations `P(a,xr) = p` and `Q(a,xr) = q` where `a` is a positive parameter 76 */ 77 function compute( a, p, q ) { 78 var ap1inv; 79 var invfp; 80 var lgama; 81 var pcase; 82 var porq; 83 var ainv; 84 var logr; 85 var ap22; 86 var ap14; 87 var ap13; 88 var ap12; 89 var vgam; 90 var vmin; 91 var xini; 92 var ap1; 93 var ap2; 94 var ap3; 95 var eta; 96 var p6; 97 var p5; 98 var x0; 99 var a2; 100 var L2; 101 var L3; 102 var L4; 103 var b2; 104 var b3; 105 var p3; 106 var a4; 107 var fp; 108 var p4; 109 var p2; 110 var a3; 111 var xr; 112 var ck; 113 var b; 114 var L; 115 var i; 116 var k; 117 var m; 118 var r; 119 var s; 120 var t; 121 var y; 122 123 if ( p < HALF ) { 124 pcase = true; 125 porq = p; 126 s = -1.0; 127 } else { 128 pcase = false; 129 porq = q; 130 s = 1.0; 131 } 132 k = 0; 133 if ( abs( a-1.0 ) < 1.0e-4 ) { 134 m = 0; 135 if ( pcase ) { 136 if ( p < 1.0e-3 ) { 137 p2 = p * p; 138 p3 = p2 * p; 139 p4 = p3 * p; 140 p5 = p4 * p; 141 p6 = p5 * p; 142 x0 = p + ( p2*HALF ) + ( p3*(ONEO3) ) + ( p4*ONEO4 ) + ( p5*ONEO5 ) + ( p6*(ONEO6) ); // eslint-disable-line max-len 143 } else { 144 x0 = -ln( 1.0-p ); 145 } 146 } else { 147 x0 = -ln( q ); 148 } 149 if ( a === 1.0 ) { 150 k = 2; 151 xr = x0; 152 } else { 153 lgama = gammaln( a ); 154 k = 1; 155 } 156 } 157 if ( q < 1.0e-30 && a < HALF ) { 158 m = 0; 159 x0 = -ln( q*gamma(a) ) + ( ( a-1.0 ) * ln( -ln( q*gamma(a) ) )); 160 k = 1; 161 lgama = gammaln( a ); 162 } 163 if ( a > 1.0 && a < 500.0 && p < 1.0e-80 ) { 164 m = 0; 165 ainv = 1.0 / a; 166 ap1inv = 1.0 / ( a+1.0 ); 167 x0 = ( gammaln( a+1.0 ) + ln( p ) ) * ainv; 168 x0 = exp( x0 ); 169 xini = x0; 170 for ( i = 0; i < 10; i++ ) { 171 x0 = xini * exp( x0*ainv ) * pow( 1.0-( x0*ap1inv ), ainv ); 172 } 173 k = 1; 174 lgama = gammaln( a ); 175 } 176 177 logr = (1.0/a) * ( ln(p) + gammaln( a+1.0 ) ); 178 if ( ( logr < ln( ONEO5 * ( 1.0+a ) ) ) && ( k === 0 ) ) { 179 r = exp( logr ); 180 m = 0; 181 a2 = a * a; 182 a3 = a2 * a; 183 a4 = a3 * a; 184 ap1 = a + 1.0; 185 ap12 = ap1 * ap1; 186 ap13 = ap1 * ap12; 187 ap14 = ap12 * ap12; 188 ap2 = a + 2.0; 189 ap22 = ap2 * ap2; 190 ap3 = a + 3.0; 191 CK[ 0 ] = 1.0; 192 CK[ 1 ] = 1.0 / ap1; 193 CK[ 2 ] = HALF * ( ( 3.0*a ) + 5.0 ) / ( ap12*ap2 ); 194 CK[ 3 ] = ONEO3 * ( 31.0 + (8.0*a2) + (33.0*a) ) / ( ap13*ap2*ap3 ); 195 CK[ 4 ] = ONEO24 * ( 2888.0 + (1179.0*a3) + (125.0*a4) + (3971.0*a2) + (5661.0*a) ) / ( ap14*ap22*ap3*( a+4.0 ) ); // eslint-disable-line max-len 196 x0 = r * evalpoly( CK, r ); 197 lgama = gammaln( a ); 198 k = 1; 199 } 200 if ( ( a < 10.0 ) && ( k === 0 ) ) { 201 vgam = sqrt( a ) / ( gamstar(a)*SQRT_TWO_PI ); 202 vmin = min( 0.02, vgam ); 203 if ( q < vmin ) { 204 m = 0; 205 b = 1.0 - a; 206 b2 = b * b; 207 b3 = b2 * b; 208 eta = sqrt( -2.0/a * ln( q/vgam ) ); 209 x0 = a * lambdaeta( eta ); 210 L = ln( x0 ); 211 if ( x0 > 5.0 ) { 212 L2 = L * L; 213 L3 = L2 * L; 214 L4 = L3 * L; 215 r = 1.0 / x0; 216 CK[ 0 ] = L - 1.0; 217 CK[ 1 ] = ( (3.0*b) - (2.0*b*L) + L2 - ( 2.0*L ) + 2.0 ) * HALF; 218 CK[ 2 ] =( (24.0*b*L) - (11.0*b2) - (24.0*b) - (6.0*L2) + (12.0*L) - 12.0 - (9.0*b*L2) + (6.0*b2*L) + (2.0*L3) ) * ONEO6; // eslint-disable-line max-len 219 CK[ 3 ] = ( (-12.0*b3*L) + (8.04*b*L2) - (114.0*b2*L) + (72.0+(36.0*L2)) + (((3.0*L4)-(72.0*L)+162.0) * (b-(168.0*b*L))) - ((12.0*L3)+(25.0*b3)) - ( (22.0*b*L3)+(36.0*b2*L2)+(120.0*b2) ) ) * ONEO12; // eslint-disable-line max-len 220 CK[ 4 ] = 0.0; 221 x0 = x0 - L + ( b*r*evalpoly( CK, r ) ); 222 } else { 223 r = 1.0 / x0; 224 L2 = L * L; 225 ck = L - 1.0; 226 t = L - (b*r*ck); 227 if ( t < x0 ) { 228 x0 -= t; 229 } 230 } 231 lgama = gammaln( a ); 232 k = 1; 233 } 234 } 235 if ( ( abs( porq-HALF ) < 1.0e-5 ) && ( k === 0 ) ) { 236 m = 0; 237 ainv = 1.0 / a; 238 x0 = a - ONEO3 + ( ( 0.0197530864197530864197530864198 + 239 ( 0.00721144424848128551832255535959*ainv ) ) * ainv ); 240 lgama = gammaln( a ); 241 k = 1; 242 } 243 if ( ( a < 1.0 ) && ( k === 0 ) ) { 244 m = 0; 245 if (pcase) { 246 x0 = exp( (1.0/a) * ( ln(porq) + gammaln(a+1.0) ) ); 247 } else { 248 x0 = exp( (1.0/a) * ( ln(1.0-porq) + gammaln(a+1.0) ) ); 249 } 250 lgama = gammaln( a ); 251 k = 1; 252 } 253 if ( k === 0 ) { 254 m = 1; 255 ainv = 1.0 / a; 256 r = erfcinv( 2.0 * porq ); 257 eta = s * r / sqrt( a*HALF ); 258 if ( r < MAX_FLOAT32 ) { 259 eta += ( eps1(eta) + ( (eps2(eta)+(eps3(eta)*ainv))*ainv ) ) * ainv; 260 x0 = a * lambdaeta(eta); 261 y = eta; 262 fp = -sqrt( a/TWO_PI ) * exp( -HALF*a*y*y ) / ( gamstar(a) ); 263 invfp = 1.0 / fp; 264 } else { 265 debug( 'Warning: Overflow problems in one or more steps of the computation.' ); 266 return NaN; 267 } 268 } 269 if ( k < 2 ) { 270 xr = higherNewton( x0, a, m, p, q, lgama, invfp, pcase ); 271 } 272 return xr; 273 } 274 275 276 // EXPORTS // 277 278 module.exports = compute;