On using Newton's method for computing the noncentrality parameter of the noncentral F distribution

On: 28 April 2014, At: 06:35 Publisher: Taylor & Francis

Informa Ltd Registered in England and Wales Registered Number: 1072954 Registered office: Mortimer House, 37-41 Mortimer Street, London W1T 3JH, UK

Communications in Statistics - Simulation and

Computation

Publication details, including instructions for authors and subscription information: http://www.tandfonline.com/loi/lssp20

On using newton's method for computing the

noncentrality parameter of the noncentral f

distribution

Cherng G. Ding a a

Institute of Management Science , National Chiao-Tung University , Taipei, Taiwan, China Published online: 27 Jun 2007.

To cite this article: Cherng G. Ding (1997) On using newton's method for computing the noncentrality parameter of the noncentral f distribution, Communications in Statistics - Simulation and Computation, 26:1, 259-268, DOI: 10.1080/03610919708813377

To link to this article: http://dx.doi.org/10.1080/03610919708813377

PLEASE SCROLL DOWN FOR ARTICLE

Taylor & Francis makes every effort to ensure the accuracy of all the information (the “Content”) contained in the publications on our platform. However, Taylor & Francis, our agents, and our licensors make no

representations or warranties whatsoever as to the accuracy, completeness, or suitability for any purpose of the Content. Any opinions and views expressed in this publication are the opinions and views of the authors, and are not the views of or endorsed by Taylor & Francis. The accuracy of the Content should not be relied upon and should be independently verified with primary sources of information. Taylor and Francis shall not be liable for any losses, actions, claims, proceedings, demands, costs, expenses, damages, and other liabilities whatsoever or howsoever caused arising directly or indirectly in connection with, in relation to or arising out of the use of the Content.

This article may be used for research, teaching, and private study purposes. Any substantial or systematic reproduction, redistribution, reselling, loan, sub-licensing, systematic supply, or distribution in any

form to anyone is expressly forbidden. Terms & Conditions of access and use can be found at http:// www.tandfonline.com/page/terms-and-conditions

C O M M U N . STATIST.-SIMULA.. 26(1), 259-268 (1997)

O N USING NEWTON'S METHOD FOR COMPUTING THE NONCENTRALlTY PARAMETER O F T H E

NONCENTRAL F DISTRIBUTlON

Cherng G. Ding

Institute of Management Science, National Chiao-Tung University,

4F, 114 Chung-Hsiao W. Road, Sec. 1, Taipei, Taiwan, Republic of China

ABSTRACT

This article deals with how to enhance computational efficiency of Newton's method for computing the noncentrality parameter of the noncentral F distribution. Some relationship between recurrence fonnulas for computing the noncentral F distribution function and its derivative with respect to the noncentrality parameter is investigated. Use of the relationship discovered can combine the evaluations of the noncentral F distribution function and its derivative, and therefore can reduce the amount of computations during Newton's process. An efficient algorithm featuring this finding is provided in a step-by-step form.

259 Copyright 8 1997 by Marcel Drkker. Inc.

INTRODIJCTION

Let F(.w; r,, r?, h ) denotc the noncent~A F distribution function with r , , r2 degrees of freedom and noncentrality parameter h. Computing the noncentrality

_,.

A

parameter h such that F(.r; I . , , r,, h) = p for given s (abscissa), r,, r, (degrees of freedom), and cumulative probability p (0 < p < F ( x ; r,, r , , 0)) always relies on numcsical root-linding methods (see, e.g., Narula and Weistroffer 1986; Guirguis

1990), in which an algorithm for evaluating F(.x; r,, r,, h) is required.

It is well known that the noncentral F distribution function can be evaluated through the noncentral beta distribution function, which can be expressed as a series of central beta distribution functions with Poisson weights. More explicitly (see, e.g., Johnson and Kotz 1970b, p.192),

F ( x ; r,, r?, h) = B ( y ;

3,

2, h ) =

C

11, (h) Iy ( 5

+

i , Q),

2 2 1=0 2 2

where B and I denote, respectively, the noncentral and central beta distribution function, y = r,.w 1 (r,x

+

r;), and u, (A) = (h/2)' exp(-h/2)l(i!). Norton (19x3)

gave a simple algorithm for evaluating (1). A bound on the e m r of truncating the series was proposed to control the computat~onal accuracy. The corresponding FORTRAN program was provided by Narula and Weistroffcr (1986). Lenth (1987) developed a recursive algorithm, whcrc only one evaluation of the central beta distribution function is required instead of several, and gave a sharper buncation error bound. Posten (1993) extracted Lenth's algorithm and provided a step-by- step listing for easy programming in any desired computer language. Smgh and Relyea (1992) proposed methods of computing (1) for integer degrees of freedom, which are based on exact expressions for l,(r,/2, rJ2) and Lenth's error hound. Ding (1994b) developed an algorithm for computing the noncentral beta dist~ibut~on function based on an alternative sesies representation. No central beta soutine is

NEWTON'S METHOD 26 1

required. The accuracy can be effectively contt.olled by using the corresponding truncation error bound.

,.

A simple approach for computing h is to use a numerical root-finding method to solve the equation F(x; r,, r?, h ) -I) = 0, where F may be evaluated by any of the algorithms mentioned above. N a ~ v l a

,.

and Weistroffer (1986) used the quadratic interpolation procedure to computed h, and claimed that it is faster than the bisection method. Guirguis (1990) introduced a modification to Newton's method, which is robust to poor starting values. He also compared the performance of Newton's method (called L-Newton), modified Newton's method (called E-

A

Newton), and quadratic interpolation for computing h , and concluded that L- Newton and E-Newton are both faster than quadratic interpolation, and E-Newton outperforms L-Newton when a good starting value is absent. Use of Newton's method, no matter with E-Newton or L-Newton, requires that w e evaluate both

F(x; r,, I.,, h ) and its derivative F'(h) with respect to

h.

T h e purpose of this

article is to investigate some relationship between recurrence formulas for evaluating

F and F' so that their evaluations can be combined at each Newton's iteration. The amount of computations during Newton's process can then be greatly reduced. T h e idea follows that given in Ding (1994a) for computing the noncentrality parameter of the noncentral chi-squared distribution, but the analytical process is

_,.

different. An algorithm for computing h based on efficient E-Newton is provided in a step-by-step form.

RELATIONSHIP BETWEEN RECURRENCE FORMULAS ,.

Computing h by using Newton's method requires the evaluations of F ( x ; r,,

r?,

A)

and F ' j h ) . It is easy to show that (see also Guirguis 1990)

Using the equation (see, e.g., Abrornowitz and Slegun 1965, p. 944) we have 1 F '( I ) = - - _u,

_{( h )}

'(('1 +

RP

+

i)

' , / 2 + i - y)r2/2 2 is0 T(r1/2

+

i

+

1) T(r.212) T((rl

+ r2)/2

+

i )

where u, = u,(h), and t, = y r1/2 + i (1 - y)'2I2. The terns T(r1/2

+

i

+

1) r(r212)

o f the series in ( 2 ) for F ' are closely related to the corresponding ones o f the series in ( 1 ) for F, which can be evaluated recursively as follows (see Posten 1993):

F(x; r,, r,, h ) = u, I,,

where l i = ly (

3

+

i, Q), i111d

2 2

NEWTON'S METHOD

The teims of thc series in (2) can be evaluated by taking advantage of the recussion given in (4).

For large values of noncentrality parameter, the lower terms of the series in

(2) and (3) make little contribution to the values of the series, and underflow error propagation may occur in recursive computations. Thus, both of the summations should be started at a higher index. A bound for the lower truncation

m - 1

error

x

u i l , (error in omitting the f i s t m terms) of the series in (3) was given by

i=0

Frick (1990) as follows:

where O denotes the standard normal distribution function. For a specified accuracy a , the starting index rn is determined by rn = [rnax(kL?

+

u ~ G

, O)].

where

U a

is the (100a)th percentage point of the standard normal distribution and

[.I

denotes 'largest integer less than or equal to '. Setting

U a

= - 5 would be sufficient since 0 ( - 5 ) < 3 x lo-'. The starting index m determined above also

m - 1 m - 1

applies for computing the series in (2) because u i ti is dominated by u i I, .

This can be seen by noting that Ii =

x

t, (see equation (4) in Ding 1994b), and

k = i

-

-

therefore ti < li for i 2 0. The upper truncation errors,

z

u,

Ii

and

z

u, ti (

i=n i=n

errors j n truncating the summations at i = n (n > m ) ) , arc both bounded above by

n. 1

I -

C

u , since

l=nl

T o approximate the series in (2) and (3), sum u p the t e r m from i = m until the error bound given in ( 6 ) is less than a predetermined accuracy E (e.g.,

COMPUTING T H E NONCENTRALITY PARAMETER

A

Newton's method is an efficient approach for computing

h.

The process is to repeat thc iteration

until I

hi+,

- h, I 5 6

A,+,,

whcrc

6

dcnotrs ,I convergence criterion (c g , 1 0 ' ~ ) . However, when we lack a good starting v;lluc, convergence may not be achieved. Guirguis (1990) proposed a modificat~on, called E-Newton, to dcal with this problem. It is given by

No matter using E-Newton or the regular Newton in (7) (called L-Newton by Guirguis), both F and F ' n e e d to be evaluated at each iterat~on. Since the factors u, and t, , computed through (4), are needed Sol- hoth F and F', and the hounds in

( 5 ) and (6) for the lower and upper trunca[ion errors are commonly used for

NEWTON'S METHOD 265

computing F and F ' , the evaluations of F and F' should be combined, rather than separated, to avoid redundant computations. In addition, the factors t! and I, d o not depend upon

h,

and hence their values can be repeatedly used for computing new iterates. Beforc Newton's iteration is performed, t, and I, are computed recursively and saved until I, 5 E. The efficiency of Newton's method can be greatly enhanced by the treatments mentioned above. Note that the evaluations of F and F' should be precise enough to ensure the accuracy of Newton's solution.

T o carry out Newton's iteration, we need a starting value h,. Cox and Reid

(1987) proposed a simple formula for approximating the percentage point of the noncentral F distribution, which can be used to obtain

h,.

The formula implies that

where Frl, (p) is the (100p)th percentage point of the Fdistribution with r., and r, degrees of freedom. Since h, is just a guess, a rough estimate of F,,, r.2 @)

would be sufficient, and could be obtained by xZ,,(p) 1 r,, where x2,,@) denote the (100p)th percentage point of the

x2

distribution with r, degrees of freedom. An approximation for x2,-,(y) is given by (see, e.g., Johnson and Kotz 1970a, p.176)

where U p could be approximated by Hastings' formula (see, e.g., Abramowitz and Stcgun 1965, p.933). In case the starting value h, is nonpositive, set h, = 6.

If the new llerate

hi+,

1s nonpositive, set

A,+,

= h,/2.

In this section, we provide an algorithm, based on E-Newton with the

propcrtics mentioned ahovc, for computing the nonccntsaliry parameter of the nonccntral F distribution, assuming that accuratc routincs for computing (the logarithm of) the gamma function

I-(.)

(see, e.g., Macleod 1989; Pikc and Hill 1966) and the central beta distribution function (see, c.g., Majunder and Bhattacharjee 1985) are available. The algorithm, named NCPF, is psesentcd in a step-by-stcp form. It can be easily convestcd to a co~nputes psogsam in any desired language. The algorithm involvcs two user-specified stopping pa~~arneters,

E and

6.

The input of the algorithm includes r,, r, (degrees of freedom), x (the abscissa), and y (the cumulative probability). Note that there will be no solution if p > F(x; rl , r2, 0). Algorithm NCPF: Step I. Set i t I, y t r,x I ( r , ~

+

r J , to t ( ( 1 + 2 ) r1/2( 1 . y p 1 2 , r(r112

+

I ) r(r2/2) and lo t 1 ~ ( 5 ,

9).

2 2

Step 2. Do Steps 2.1-2.2 until I, 5 E.

2.1 Compute Ii t Ii., - t ,.,, and

r ;

t t,., y ((r,

+

rJ2

+

i - l ) / ( r , / 2

+

i).

Step 3. Set N t i - 1. Then obtain an initial guess h, and set h t h, Step 4. Do Steps 4.1-4.7 until IDIFA I 61.

4.1 Set In t [max (XI2 - 5 6 , O)].

4.2 Set u t exp(1n ln(hl2) - 112 - In T(m+ I)), v t u , CDF t 11 I,,, , and DER t 14 t,,, .

4.3.1 Set u t uhl(2i), v t v

+

u , CDF t CDF

+

u

r , ,

and DER t

DER

+

uI, .

NEWTON'S METHOD

4.5 Do Steps 4.5.1 -4.5.2 until 1 - 5

e.

4.5.1 Cornpule I* t I* - I * , I* t- t * y ( ( r ,

+

r2)/2+i - l ) / ( r , / 2 + i ) ,

11 t uh l(2i), v t 11

+

1 1 , CDF t CDF

+

u t*, and DER t DER

+

u I*.

4.6 Set DIFF t 2 CDF 111(1,/ CDF) I DER.

4.7 If h

+

DIFF 5 0 , set h t h 12; othe~wise set h t h

+

DIFF A

Step 5 . Stop with h t h.

B I B L I O G R A P H Y

Abramowitz, M. and Stegun, I. A. (1965). Handbook of MathematicalFunctions. New York: Dover Publications, Inc.

Cox, D. R. and Reid, N. (1987). "Approximations to noncentral distributions,"

Canad. J. Statist., 15, 105-1 14.

Ding, C. G. (1994a). "An efficient algorithm for computing the noncentrality parameters of chi-squared tests," Colnlnun. Statist.-Sirnula., 23, 861-870. Ding, C. G. (1994b). "On the computation of the noncentral beta distribution,"

Cornp. Statist. & DataAnal., 18, 449-455.

Frick, H. (1990). "A remark on Algorithm AS 226: Computing noncentral beta probabilities," Appl. Statist., 39, 3 11-312.

Guirguis, G. H. (1990). "A note on computing the noncentrality parameter of the noncentral F-distribution," Coinmun. Statist. -Sirnula., 19, 1497- 15 1 1.

Johnson, N. L. and Kotz, S. (1970a). Distribution in Statistics: Continuous

Univariate Distributions

-

1. New York: John Wiley & Sons, Inc.

Johnson, N. L. and Kotz, S. (1970b). Distribution in Statistics: Continuous

Univarinte Distributions - 2 . New Yosk: John Wiley & Sons, Inc.

Lenth, R. V . (1987). "Computing noncentral beta probabilities," Appl. Statist., 36, 241-244.

Macleod, A. J. (1989). "A robust and reliable algorithm for the logarithm of the gamma function." Appl. Stotist., 38, 397-402.

Majundcr, K. L. and Bhattacharlee, C. P. (1985). "The incomplete beta integral,"

Applied Stcltistics Algorithms (P. Griffiths and I. D. Hill, Eds.). Chichester: Ellis Honvood, 1 17- 120.

Narula, S. C. and Weistroffer, H. R. (1986). "Computation of probability and non-centrality parameter of a non-central F-distribution," Commun. Statist. -Sirnula., 15, 871-878.

Norton, V. (1983). "A simple algorithm for computing the non-central F

distribution," Appl. Stotist., 32, 84-85.

Pike, M. C. and Hill, I. D. (1966). "Logarithm of gamma function," Commun. Ass. Comput. Mach.. 9, 684.

Posten,

H.

0. (1993). "An effective algorithm for the noncentral beta distribution function," Arner. Statistn., 47, 129- 13 1.

Singh, K. P. and Relyea, G. E. (1992). "Computation of nonccntral f p~ohahilitics:

A computer program," Cntnp. Str~tixt. & Data A n d . , 13, 95- 102.

R e c e i v e d A u g u s t , 1994; R e v i s e d May, 1996.