Efficient algorithms for the inversion of the cumulative central beta distribution
aa r X i v : . [ m a t h . NA ] M a y Efficient algorithms for the inversion of thecumulative central beta distribution
A. GilDepartamento de Matem´atica Aplicada y CC. de la Computaci´on.ETSI Caminos. Universidad de Cantabria. 39005-Santander, Spain.J. SeguraDepartamento de Matem´aticas, Estad´ıstica y Computaci´on,Univ. de Cantabria, 39005 Santander, Spain.N.M. TemmeIAA, 1391 VD 18, Abcoude, The Netherlands ∗ Abstract
Accurate and efficient algorithms for the inversion of the cumulativecentral beta distribution are described. The algorithms are based onthe combination of a fourth-order fixed point method with good non-local convergence properties (the Schwarzian-Newton method), asymp-totic inversion methods and sharp bounds in the tails of the distribu-tion function.
The cumulative central beta distribution (also known as the incomplete betafunction) is defined by I x ( p, q ) = 1 B ( p, q ) Z x t p − (1 − t ) q − dt, (1.1)where we assume that p and q are real positive parameters and 0 ≤ x ≤ B ( p, q ) is the Beta function ∗ Former address: Centrum Wiskunde & Informatica (CWI), Science Park 123, 1098XG Amsterdam, The Netherlands ( p, q ) = Γ( p )Γ( q )Γ( p + q ) . (1.2)From the integral representation in (1.1) it is easy to check the followingrelation: I x ( p, q ) = 1 − I − x ( q, p ) . (1.3)In this paper we describe algorithms for solving the equation I x ( p, q ) = α, < α < , (1.4)with p , q given fixed real positive values. In statistical terms, we are com-puting the quantile function for I x ( p, q ).The beta distribution ia a standard and widely used statistical distribu-tion which has as particular cases other important distributions like the Stu-dent’s distribution, the F-distribution and the binomial distribution. There-fore, the computational schemes for inverting the central beta distributioncan be used to compute percentiles for other distributions related to thebeta. For an example for the F-distribution see [1].The quantile function is useful, for instance, for the generation of ran-dom variables following the beta distribution density. In some Monte Carlosimulations the generation of such random variables are required and a mas-sive number of inversions of the beta cumulative distribution are needed.Therefore, it is important to construct methods as reliable and efficient aspossible.Existing algorithms use some simple initial approximations which areimproved by iterating with the Newton method. In particular, this is theapproach used in the inversion method of the statistical software package R , which is based on the algorithm of [9] and the succesive improvementsand corrections [4, 2, 3]. In [9], a simple approximation is used in terms ofthe error function together with two additional starting value approxima-tions for the tails; these initial values are refined by the Newton iteration.As discussed in [4, 2], the Newton method needs some modification to en-sure convergence inside the interval [0 , R (butsome convergence problem still remain in the present version, as we laterdiscuss).In this paper the methods for the computation of the inversion of thecumulative beta distribution are improved in several directions. In the first2lace, we introduce the Schwarzian-Newton method (SNM) as alternative toNewton’s method (NM). With respect to Newton’s method the SNM has theadvantage of having order of convergence four instead of two. In addition,as explained in [10], the SNM has good non-local properties for this typeof functions and it is possible to build an algorithm with certified conver-gence. In the second place, we analyze initial value approximations (muchsharper than those given in [9]) in terms of asymptotic approximations forlarge values of p and/or q , but which also give accurate values for moderatevalues; these approximations are given in terms of inverse error functions orthe inverse gamma distribution ([12], [6, § § We next describe the methods of computation used in the algorithms. Firstwe describe the SNM method, and discuss how a standalone algorithm withcertified convergence can be built with this method, provided an accuratemethod of computation of the beta distribution is available. In the secondplace we describe the methods for estimating the quantile function basedon asymptotics for large p and/or q . Finally, we describe sharp upper andlower bounds for the tails that can be used for solving the problem (1.4) for α close to zero or 1. The Schwarzian-Newton method (SNM) is a fourth order fixed point methodwith good non-local convergence properties for solving nonlinear equations f ( x ) = 0 [10]. The SNM has Halley’s method as limit when the Schwarzianderivative of the function f ( x ) tends to zero.Given a function f ( x ) with positive derivative (in our case f ( x ) = I x ( p, q ) − α ), it is easy to prove that Φ( x ) = f ( x ) / p f ′ ( x ) satisfies the3ifferential equationΦ ′′ ( x ) + Ω( x )Φ( x ) = 0 , Ω = 12 { f, x } , (2.1)where { f, x } is the Schwarzian derivative of f ( x ) with respect to x : { f, x } = f ′′′ f ′ − f ′′ f ′ ! . (2.2)The SNM is obtained from the approximate integration of the Riccatiequation h ′ ( x ) = 1 + Ω( x ) h ( x ) , h ( x ) = Φ( x ) / Φ ′ ( x ) under the assumptionthat Ω( x ) is approximately constant. In the case of negative Schwarzianderivative (which will be the case for the beta distribution) the iterationfunction can be written as: g ( x ) = x − p | Ω | arctanh (cid:16)p | Ω | h ( x ) (cid:17) . (2.3)We discuss two implementations: a direct implementation, which gives aconvergent algorithm for p, q > It is proved in [10] that if Ω( x ) has one and only one extremum at x e ∈ I and it is a maximum, then if Ω < f ( x ) in I starting from x = x e . We will use this result for thecumulative central beta distribution, when the parameters p and q are largerthan 1. In this case, the function Ω( x ) (the Schwarzian derivative of f ( x )with a factor 1 /
2) is given byΩ( x ) = ( p − q − x (1 − x ) − p − x − q − − x ) . (2.4)It is possible to show that for p > q >
1, the function Ω( x ) in (2.4)is negative in (0 ,
1) and has only one extremum (which is a maximum) inthat interval. The extremum of Ω( x ) is at x e = 13∆ / (3 pq + 3 p + 6 p )∆ / − ∆ / + 3 pq (cid:0) ( p + q ) + 8( p + q ) + 12 (cid:1) ( p + q ) + 2( p + q ) , (2.5)4here∆ = pq (cid:8) p − q )( p + q + 1) + 27( p q + p − q p − q ) + 3 √ p + q )( p + q + 2) p ( p + q + 2)(27 q + 54 + q p + 18 pq + 27 p + p q ) o . Then, the fixed point method is (2.3) with Ω( x ) given by (2.4) and h ( x ) = f ( x )12 (cid:18) − p − x + q − − x (cid:19) f ( x ) + x p − (1 − x ) q − B ( p, q ) , where f ( x ) = I x ( p, q ) − α . When p and/or q are smaller than 1, it is possible to make a change ofvariables in order to obtain a negative Schwarzian derivative and simplermonotonicity properties. In particular, with the change of variables z ( x ) = log (cid:18) x − x (cid:19) , (2.6)we obtain that Φ( z ) = f ( z ) / q ˙ f ( z ) (where the dot represents the derivativewith respect to z ) satisfies ¨Φ( z ) + Ω( z )Φ( z ) = 0, whereΩ( z ) = 14 (cid:0) − ( p + q )( p + q − x ( z ) + 2( p + q )( p − x ( z ) − p (cid:1) . (2.7)The function Ω( z ) has its extremum at z e = log( x e / (1 − x e )), x e =( p − / ( p + q − p and/or q smaller than 1, Ω( z ) can be either bemonotonic or it can have a minimum. Convergence of the SNM can be guar-anteed, in this case, using the following results [10]: a) if Ω( z ) is negativeand decreasing in the interval I = [ a, α ], then the SNM converges mono-tonically to α for any starting value z ∈ [ a, α ]; b) if Ω( z ) is negative andincreasing in the interval I = [ α, b ], then the SNM converges monotonicallyto α for any starting value z ∈ [ α, b ].The case p = q = 1 is of course trivial. Apart from this, there are threedifferent cases to be considered:a) p ≤ , q >
1: the function Ω( z ) is decreasing. In this case, the SNMuses as starting point a large negative value (in the z variable).5) p > , q ≤
1: the function Ω( z ) is increasing. In this case, the SNMuses as starting point a large positive value (in the z variable).c) p < , q <
1: the extremum of Ω( z ) at z e is reached and it is a mini-mum. In this case, we use the sign of the function h ( z ) at z e to selecta subinterval for application of the SNM, according to the previousresults. The function h ( z ) is given by h ( z ) = f ( x ( z )) − p − qe z e z f ( x ( z )) + e zp B ( p, q )(1 + e z ) p + q . (2.8)When this sign is negative, the SNM uses a large positive value (in the z variable), otherwise the SNM uses a large negative value.Once the SNM is applied to find the root z r in the z -variable, the corre-sponding x -value will be given by x r = e z r e z r . We have constructed two methods of order four which are proven to convergewith certainty for the initial values prescribed. The method has, in addition,good non-local properties, which means that few iterations are needed fora good estimation of the inverse (typically from 2 to 4 for 20 digits), evenwithout accurate starting values. The exceptions are the tails ( α very closeto 0 or 1), but we will discuss later how to deal with these cases.Because the convergence is guaranteed, no special precaution is neededto ensure that the interval [0 ,
1] in the original variable is not abandoned, ashappened with earlier versions of the algorithm given in [9] (see [4]) and as itis still happens for some values in the latests R version of this algorithm. Forinstance, the R command qbeta(alpha,600,1.1) does not converge properlyif alpha ∈ (6 . − , . − ). Our method avoids this type of problems.The performance of the method can be improved by considering initialapproximations, which we are discussing next. The algorithm considered in [9], which is the basis of the R implementation,uses an approximation in terms of the inverse error function, which works6easonably away from the tails. However, this simple approximation doesnot give more than two accurate digits, except by accident.Much more accurate initial approximations (some of them also in termsof error functions) can be obtained from asymptotics for large p and/or q .These are accurate approximations for large and not so large p and/or q , aswe later discuss.This section is based on the results given in [12] and [13, § We start with the following representation I x ( p, q ) = erfc (cid:16) − η p r/ (cid:17) − R r ( η ) , (2.9)where we write p = r sin θ , q = r cos θ with 0 < θ < π/ η is given by − η = sin θ log x sin θ + cos θ log 1 − x cos θ . (2.10)When we take the square root for η , we choose sign( η ) = sign( x − sin θ ),this means sign( η ) = sign( x − p/ ( p + q )). In this way, the relation between x ∈ (0 ,
1) and η ∈ ( −∞ , ∞ ) becomes one-to-one.Using this representation of I x ( p, q ), we solve the equation in (1.4) firstin terms of η . When r = p + q is a large parameter, the asymptotic methodwill provide an approximation to the requested value of η in the form η ∼ η + η r + η r + η r + . . . . (2.11)The algorithm for computing the coefficients η i , i = 0 , , , . . . can besummarized as follows1. The value η is obtained from the equation erfc (cid:16) − η p r/ (cid:17) = α. (2.12)2. With η = η , equation (2.10) is inverted to obtain a first approximationto the value of x . For inverting this equation, it seems convenient towrite it in the form x p (1 − x ) q = (cid:18) pr (cid:19) p (cid:18) qr (cid:19) q e − rη / . (2.13)7. With these values of η and x , the coefficient η is given by η = log( f ( η )) η , (2.14)where f ( η ) = η sin θ cos θ ( x − sin θ ) .4. Higher-order coefficients η j , j = 2 , , . . . can be obtained in terms of x , η , η , sin θ and cos θ . As an example, the coefficient η is given by η = 112 η c s ( s − x ) (cid:0) s η − η x − s η − η s +12 s c − s c η η x + 12 s c η η x − η s c η x +12 η s c η x + 2 η xs + 2 η xs − η s c η +12 s c η x − s c η x − η xs − η x s + η x s − s c x + 12 s c x (cid:1) , (2.15)where s = sin θ , c = cos θ .5. With these coefficients in the expansion (2.11), a value for η is ob-tained. Then, the inversion of (2.10) will provide the final asymptoticestimation of the x -value.Using (2.10) we can derive the following expansion or small values of | η | : x = s + sc η + 1 − s η + 13 s − s + 136 sc η +46 s − s + 21 s + 1270 s c η + . . . , (2.16)where s = sin θ , c = cos θ . For larger values of | η | , with η <
0, we rewrite(2.10) in the form x (1 − x ) µ = u , where µ = cot θ, u = exp h ( − η + s ln s + c ln c ) /s i , (2.17)and for small values of u we expand x = u + µu + 3 µ (3 µ + 1)3! u + 4 µ (4 µ + 1)(4 µ + 2)4! u +5 µ (5 µ + 1)(5 µ + 2)(5 µ + 3)5! u + . . . . (2.18)8 similar approach is possible for positive values of η , giving an expansionfor x near unity. In that case we have the equation x ν (1 − x ) = v , where ν = tan θ, v = exp h ( − η + s ln s + c ln c ) /c i , (2.19)and we have the expansion1 − x = v + νv + 3 ν (3 ν + 1)3! v + 4 ν (4 ν + 1)(4 ν + 2)4! v +5 ν (5 ν + 1)(5 ν + 2)(5 ν + 3)5! v + . . . , (2.20)The approximations to x obtained in this way will be used for startingthe SNM for obtaining more accurate values of x . In this case, we start from I x ( p, q ) = Q ( q, ηp ) + R p,q ( η ) , (2.21)where Q ( a, x ) is the incomplete gamma function ratio.The parameter η is given by η − µ log η + (1 + µ ) log(1 + µ ) − µ = − log x − µ log(1 − x ) , (2.22)where µ = q/p and x have the following corresponding points: x = 0 ⇐⇒ η = + ∞ x = 1 / (1 + µ ) ⇐⇒ η = µ, x = 1 ⇐⇒ η = 0 . (2.23)So, for x ∈ (0 ,
1) we have sign( η − µ ) = − sign( x − / (1 + µ )).For the representation in (2.21) we assume that p is the large parameter,and we will obtain approximations to the value of η in the form η ∼ η + η p + η p + η p + . . . , p → ∞ . (2.24)We follow similar ideas as in § η can be obtained bysolving Q ( q, η p ) = α. (2.25)The inversion of Q ( a, x ) can be done by using our inversion algorithmdescribed in [7]. 9hen, a value x is obtained by solving (2.22) for x . With x and η wecompute η = log φ ( η )1 − µ/η , (2.26)where φ ( η ) is given by φ ( η ) = η − µ − x (1 + µ ) 1 p µ .Other coefficients η j , j = 2 , , . . . can be obtained in terms of µ , x , η and η .To compute x from equation (2.22) we can use the inversion algorithm forcomputing x when in (2.10) η is given. This follows from µ = q/p = cot θ and from writing (2.22) in the formsin θ (cid:18) µ − η + µ log ηµ (cid:19) = sin θ log x sin θ + cos θ log 1 − x cos θ . (2.27)This equation can also be written as x p (1 − x ) q = (cid:18) pr (cid:19) r e q (1+log η ) − pη . (2.28) Sharp lower and upper bounds for the solution x of the equation (1.4) in thelower ( α →
0) and upper ( α →
1) tails of the distribution function can beobtained by using fixed point iterations x n +1 l = g l ( x nl ) and x n +1 u = g u ( x nu ),respectively, where the iteration functions g l ( x ) and g u ( x ) for the lower tailare given by [11] g l ( x ) = (cid:0) αB ( p, q ) ( p − ( p + q ) x ) (1 − x ) − q (cid:1) /p , (2.29)and g u ( x ) = αpB ( p, q ) p + q )( p + 1) x + ( p + q )( p + q + 1)( p + 1)( p + 2) x ! (1 − x ) q /p . (2.30)The starting value of the fixed point iterations is x = 0. The solution x of the equation (1.4) satisfies x l < x < x u . These bounds of x can10e used as starting values for the SNM. Notice that, because a lower andan upper bound is obtained we have an estimation of the error for theseapproximations that can be used to decide if they are accurate enough.We notice that the approximation for the lower tail used in [9] (and alsoin the R implementation) is just g l (0) = g u (0). Our approximation is moreaccurate and provides upper and lower bounds.For the upper tail, the iteration functions are the same as before, butwith p and q interchanged (by using (1.3)). The bounds are then given for1 − x . In this section we illustrate the use of the different methods with numericalexamples which will help in deciding how to combine the methods in orderto obtain fast reliable algorithms (which are described in section 4).
In Tables 1 and 2 we show examples of the accuracy obtained in the com-putation of | I x ( p, q ) − α | /α with x the values provided by the asymptoticapproximations (before iterating the SNM). The asymptotic methods provide relatively good initial values even forquite small values of p and q : using 10 random points in the region ( p, q, α ) ∈ (0 . , . × (0 . , . × (0 ,
1) we have tested that a relative accuracy betterthan 0 .
06 was obtained when computing | I x ( p, q ) − α | /α with x , the asymp-totic approximations obtained using the error function. With these initialvalues, not more than two iterations of the SNM are needed to obtain anaccuracy better than 5 . − . The function I x ( p, q ) is computed in theiterations of the SNM by using a continued fraction representation I x ( p, q ) = x p (1 − x ) q pB ( p, q ) d d d . . . ! , (3.1)where In Table 1 the accuracy 5 . − corresponds to the case I x (3 ,
3) = 0 .
5: because ofthe symmetry, x should be 0 . η defined in (2.10) becomes 0, the same as η in (2.12). This explains why that result in the table becomes so small. p = 4 p = 3 p = 210 − . − . − . − − . − . − . − . . − . − . − . . − . − . − . . − . − . − . . − . − . − . . − . − . − .
999 4 . − . − . − . . − . − . − Table 1 . Relative errors | I x ( p, q ) − α | /α for r = p + q = 6 using the estimates providedby the asymptotic inversion method with the error function. d m = m ( q − m ) x ( p + 2 m − p + 2 m ) ,d m +1 = − ( p + m )( p + q + m ) x ( p + 2 m )( p + 2 m + 1) . (3.2)For the computation of the Beta function B ( p, q ), it is convenient, inparticular when p and q are large, to use the following expression in termsof the scaled gamma function Γ ∗ ( x ): B ( p, q ) = √ π s p + 1 q p pp + q q qp + q p + q p + q Γ ∗ ( p )Γ ∗ ( q )Γ ∗ ( p + q ) , (3.3)where Γ ∗ ( x ) is defined asΓ ∗ ( x ) = Γ( x ) p π/x x x e − x , x > . (3.4)The function Γ ∗ ( x ) is computed using the function gamstar included ina previous package developed by the authors [8]. The interval estimations in the tails of the central beta distribution functionby using the fixed point iterations of § p and q , where the asymptotic method cannot be12 µ = 0 . µ = 0 . µ = 210 − . − . − . − − . − . − . − . . − . − . − . . − . − . − . . − . − . − . . − . − . − . . − . − . − .
999 1 . − . − . − . . − . − . − Table 2 . Relative errors | I x ( p, q ) − α | /α for p = 7 and several values of µ = q/p usingthe estimates provided by the asymptotic inversion method with the incomplete gammafunction. applied. It is important to note that when the value of the parameters p or q are close to 0, the inversion of I x ( p, q ) becomes problematic in the lower(when p →
0) or upper (when q →
0) tail of the cumulative distributionfunction, because of the particular shape of the functions.In Table 3 we show the relative errors 1 − x l /x and 1 − x u /x obtainedwith the lower and upper bounds, respectively, for the solution x of theequation (1.4) for small values of p , q and α . The bounds (computed in theexamples using Maple) are obtained with just three iterations of the fixedpoint methods of § p and q ,the bounds provide in all cases reasonable approximations for starting theSNM, no matter if the value of α is small. Besides, even for not so smallvalues of p and q , the bounds provide very accurate estimations when α isvery small. In some cases, these estimations could be even better than theestimations of the asymptotic method. We have tested that the scheme for the SNM, as described in § p and q are both small, also provides a good uniform accu-racy in the computation: using 10 random points in the region ( p, q, α ) ∈ (0 . , . × (0 . , . × (0 ,
1) we have tested that a relative accuracy betterthan 4 . − was obtained when computing | I x ( p, q ) − α | /α . The maxi-mum number of iterations of the SNM was 3.13 p = 0 . p = 0 . q = 0 . q = 0 . − (LB) − . − − . − (UB) 7 . − . − − (LB) − . − − . − (UB) 1 . − . − − (LB) − . − − . − (UB) 8 . − . − Table 3 . Relative errors 1 − x l /x and 1 − x u /x obtained with the lower (LB) and upper(UB) bounds for the solution x of the equation (1.4). As we have shown in § p and q . But apart from accuracy, an important feature of any computationalscheme is also efficiency. So, we have compared whether the combined use ofthe asymptotic approximations plus iterations of the SNM is more efficientor not than the sole use of the SNM. In Table 4 we show CPU times spent by20000 runs of the inversion algorithm for different values of p , q and α usingthree methods of computation: a) asymptotic inversion method using theerror function for estimating the initial value plus iterations of the SNM;b) asymptotic inversion method using the gamma function for estimatingthe initial value plus iteration of the SNM; c) iterations of the SNM withstarting values obtained as discussed in § α are neither very small nor near to 1.A different scenario arises when the solution of equation (1.4) is com-puted with an aimed accuracy better than 10 − (single precision). In thiscase, just using the approximations provided by the asymptotic expansionswill be enough to obtain such an accuracy for a wide range of parameters.14 Method p = 4 p = 50 , p = 100 , p = 150 , p = 300 ,q = 3 q = 60 q = 80 q = 1 . q = 40010 − M1 0 .
22 0 .
29 0 .
22 0 . . .
31 0 .
39 0 .
31 0 .
22 0 . .
32 0 .
34 0 .
34 0 .
39 0 . − M1 0 .
22 0 .
23 0 .
19 0 .
36 0 . .
34 0 .
34 0 .
31 0 .
33 0 . .
25 0 .
27 0 .
28 0 .
27 0 . . .
19 0 .
16 0 .
17 0 .
20 0 . .
31 0 .
28 0 .
30 0 .
19 0 . .
11 0 .
19 0 .
20 0 .
19 0 . . .
20 0 .
17 0 .
16 0 .
23 0 . .
33 0 .
28 0 .
30 0 .
23 0 . .
16 0 .
19 0 .
19 0 .
16 0 . .
999 M1 0 .
17 0 .
17 0 .
17 0 .
14 0 . .
25 0 .
28 0 .
28 0 .
16 0 . .
25 0 .
27 0 .
28 0 .
16 0 . Table 4 . CPU times (in seconds) for 20000 runs of the inversion algorithm using dif-ferent methods. M1: Asymptotic inversion method using the error function +SNM; M2:Asymptotic inversion method using the gamma function +SNM; M3: SNM. The SNM isiterated until the solution of equation (1.4) is obtained with an accuracy near full doubleprecision. Proposed algorithms
Based on the previous numerical experiments of §
3, we conclude that ifthe precision required is not very high, the initial approximations given byasymptotics or by the tail estimations could be sufficient in a large range ofthe parameters. However, for higher precision the use of the SNM must beprevalent.This leads us to suggest two different schemes for computing the solution x of the equation (1.4).SCHEME 1. Algorithm for the inversion of the cumulative central betadistribution with an accuracy near double precision:If α ≤ . p < .
3, use the upper bound of § . < p <
1, use the SNM as described in § § < p <
30 and q <
1, use the SNM, using as starting valuesthe bounds of § p >
30 and q < .
5, use the SNM, using as starting valuesthe bounds of § p >
30, 0 . < q <
5: a) if α > . α < . § . < α ≤ . < q < p >
50, use the SNM, using as starting val-ues the approximations provided by the uniform asymptoticexpansion in terms of the incomplete gamma function.For p >
30 and q >
30, use the SNM, using as starting val-ues the approximations provided by the uniform asymptoticexpansion in terms of the error function.In other cases, use the SNM as described in § . < α <
1, use the relation (1.3) and apply the previous stepsto solve 1 − x in I − x ( q, p ) = 1 − α .SCHEME 2. Algorithm for the inversion of the cumulative central betadistribution with an accuracy near single precision:If α ≤ . p < .
5, use the upper bound of § . < p <
1, use the SNM as described in § § < p <
30 and q <
1, use the SNM, using as starting valuesthe bounds of § p >
30 and q < .
5, use the SNM, using as starting valuesthe bounds of § p >
30, 0 . < q <
5: a) if α > . α < . § § . < α ≤ . < q < p >
160 and α > .
1, use the approximationprovided by the uniform asymptotic expansion in terms ofthe incomplete gamma function as solution of the equation(1.4).For p >
30 and q >
30, use the approximation provided by theuniform asymptotic expansion in terms of the error functionas solution of the equation (1.4).In other cases, use the SNM as described in § . < α <
1, use the relation (1.3) and apply the previous stepsto solve 1 − x in I − x ( q, p ) = 1 − α . We have presented methods for the inversion of cumulative beta distribu-tions which improve previously existing methods. We have described how17he Schwarzian-Newton method provides a standalone method with certifiedconvergence which is in itself an efficient method, even without accurate ini-tial estimations (except at the tails). In addition, we have discussed how toimprove the efficiency by estimating the quantile function using asymptoticsfor large p and/or q and by considering sharp upper and lower bounds forthe tails. These initial estimations are considerably more accurate than thesimple approximations used in some standard mathematical software pack-ages (like R ) and, combined with the fourth order SNM, provide efficientand reliable algorithms for the inversion of cumulative beta distributions. The authors acknowledge financial support from
Ministerio de Econom´ıa yCompetitividad , project MTM2012-34787. NMT thanks CWI, Amsterdam,for scientific support.
References [1] R.W. Abernathy and R.P. Smith. Algorithm 724: Program to CalculateF-Percentiles.
ACM Trans Math Soft , 19(4):481–483, 1993.[2] K. J. Berry, P. W. Mielke, and G. W. Cran. Algorithm as r83: Aremark on algorithm as 109: Inverse of the incomplete beta functionratio.
Appl. Statist. , 39(2):309–310, 1990.[3] K. J. Berry, P. W. Mielke, and G. W. Cran. Correction to algorithmas r83-a remark on algorithm as 109: Inverse of the incomplete betafunction ratio.
Appl. Statist. , 40(1):236, 1991.[4] G. W. Cran, K. J. Martin, and G. E. Thomas. Remark as r19 andalgorithm as 109: A remark on algorithms: As 63: The incompletebeta integral as 64: Inverse of the incomplete beta function ratio.
Appl.Statist. , 26(1):111–114, 1977.[5] A.R. Didonato and A.H. Morris. Algorithm 708: Significant digit com-putation of the incomplete beta function ratios.
ACM Trans. Math.Softw. , 18(3):360–373, 1992.[6] A. Gil, J. Segura, and N. M. Temme.
Numerical Methods for SpecialFunctions . Society for Industrial and Applied Mathematics (SIAM),Philadelphia, PA, 2007. 187] A. Gil, N. M. Temme, and J. Segura. Efficient and accurate algorithmsfor the computation and inversion of the incomplete gamma functionratios.
SIAM J. Sci. Comput. , 34(6):A2965–A2981, 2012.[8] A. Gil, N. M. Temme, and J. Segura. Gammachi: a package for theinversion and computation of the gamma and chi-square cumulative dis-tribution functions (central and noncentral).
Comput. Phys. Commun. ,191:132–139, 2015.[9] K. L. Majumder and G. P. Bhattacharjee. Algorithm as 64: Inverse ofthe incomplete beta function ratio.
Appl. Statist. , 2(3):411–414, 1973.[10] J. Segura. The Schwarzian-Newton method for solving nonlinear equa-tions, with applications.
Math. Comput. (to appear) .[11] Javier Segura. Sharp bounds for cumulative distribution functions.
J.Math. Anal. Appl. , 436(2):748–763, 2016.[12] N. M. Temme. Asymptotic inversion of the incomplete beta function.
J. Comput. Appl. Math. , 41(1-2):145–157, 1992.[13] N. M. Temme.