226 Chapter 6.Special Functions TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT (0.5,5.0) (8.010.0) /(1.0,3.0) Permission is 6 (0.5,0.5) .com or call 1-800-872- granted for (including this one) internet to any server computer, t users to make one paper 1988-1992 by Cambridge University Press. from NUMERICAL RECIPES IN C: (5.0,0.5) 0 (North America THE 0 .6 Programs Figure 6.4.1.The incomplete beta function (,b)for five different pairs of (a,b).Notice that the pairs (0.5,5.0)and(5.0,0.5)are symmetrically related as indicated in equation (6.4.3). 6.4 Incomplete Beta Function,Student's Distribution,F-Distribution,Cumulative 1788-1982 ART OF SCIENTIFIC COMPUTING (ISBN 0-521 Binomial Distribution The incomplete beta function is defined by Numerical Recipes -43108-5 I(a,b)= Ba,0=1 B(a,b) B(a,b)Jo ta-1(1-t)b-1dt (a,b>0) (6.4.1) (outside Software. It has the limiting values North Io(a,b)=0 11(a,b)=1 (6.4.2) visit website f machine and the symmetry relation Iz(a,b)=1-11-x(b,a) (6.4.3) If a and b are both rather greater than one,then I(a,b)rises from"near-zero"to "near-unity"quite sharply at about x =a/(a+b).Figure 6.4.1 plots the function for several pairs (a,b)
226 Chapter 6. Special Functions Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). 0 (5.0,0.5) (0.5,0.5) (8.0,10.0) (1.0,3.0) (0.5,5.0) .2 .4 .6 1 .8 0 .2 .4 .6 .8 1 incomplete beta function Ix(a,b) x Figure 6.4.1. The incomplete beta function Ix(a, b) for five different pairs of (a, b). Notice that the pairs (0.5, 5.0) and (5.0, 0.5) are symmetrically related as indicated in equation (6.4.3). 6.4 Incomplete Beta Function, Student’s Distribution, F-Distribution, Cumulative Binomial Distribution The incomplete beta function is defined by Ix(a, b) ≡ Bx(a, b) B(a, b) ≡ 1 B(a, b) x 0 t a−1(1 − t) b−1dt (a, b > 0) (6.4.1) It has the limiting values I0(a, b)=0 I1(a, b)=1 (6.4.2) and the symmetry relation Ix(a, b)=1 − I1−x(b, a) (6.4.3) If a and b are both rather greater than one, then I x(a, b) rises from “near-zero” to “near-unity” quite sharply at about x = a/(a + b). Figure 6.4.1 plots the function for several pairs (a, b)
6.4 Incomplete Beta Function 227 The incomplete beta function has a series expansion Iz(a,b)= xa(1-z)6 1+Be+1,n+1 rn+1 (6.4.4) aB(a,b) B(a+b,n+1) n=0 but this does not prove to be very useful in its numerical evaluation.(Note,however, that the beta functions in the coefficients can be evaluated for each value of n with just the previous value and a few multiplies,using equations 6.1.9 and 6.1.3.) The continued fraction representation proves to be much more useful. 三 Iz(a,b)= x(1-x)b「1dd2 aB(a,b) 1+1+1+ (6.4.5) where (a+m)(a+b+m)x d2m+1=- (a+2m)(a+2m+1) (6.4.6) d2m三 m(b-m)x (a+2m-1)(a+2m) This continued fraction converges rapidly for x(a +1)/(a+b+2)we can 、。 令 Press. just use the symmetry relation(6.4.3)to obtain an equivalent computation where the ART continued fraction will also converge rapidly.Hence we have Programs #include ictly proh float betai(float a,float b,float x) Returns the incomplete beta function Ix(a,b) float betacf(float a,float b,float x); to dir float gammln(float xx); void nrerror(char error_text []) 1881892 OF SCIENTIFIC COMPUTING(ISBN float bt; if (x 0.0 II x 1.0)nrerror("Bad x in routine betai"); v@cam 1f(x==0.011x==1.0)bt=0.0; else Factors in front of the continued fraction. bt=exp(gammln(a+b)-gammln(a)-gammln(b)+a*log(x)+b*log(1.0-x)); if(x #define MAXIT 100 #define EPS 3.0e-7 #define FPMIN 1.0e-30 float betacf(float a,float b,float x) Used by betai:Evaluates continued fraction for incomplete beta function by modified Lentz's method (85.2). void nrerror(char error.-text☐)；
6.4 Incomplete Beta Function 227 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). The incomplete beta function has a series expansion Ix(a, b) = xa(1 − x)b aB(a, b) 1 + ∞ n=0 B(a + 1, n + 1) B(a + b, n + 1) xn+1 , (6.4.4) but this does not prove to be very useful in its numerical evaluation. (Note, however, that the beta functions in the coefficients can be evaluated for each value of n with just the previous value and a few multiplies, using equations 6.1.9 and 6.1.3.) The continued fraction representation proves to be much more useful, Ix(a, b) = xa(1 − x)b aB(a, b) 1 1+ d1 1+ d2 1+ ··· (6.4.5) where d2m+1 = − (a + m)(a + b + m)x (a + 2m)(a + 2m + 1) d2m = m(b − m)x (a + 2m − 1)(a + 2m) (6.4.6) This continued fraction converges rapidly for x (a + 1)/(a + b + 2) we can just use the symmetry relation (6.4.3) to obtain an equivalent computation where the continued fraction will also converge rapidly. Hence we have #include float betai(float a, float b, float x) Returns the incomplete beta function Ix(a, b). { float betacf(float a, float b, float x); float gammln(float xx); void nrerror(char error_text[]); float bt; if (x 1.0) nrerror("Bad x in routine betai"); if (x == 0.0 || x == 1.0) bt=0.0; else Factors in front of the continued fraction. bt=exp(gammln(a+b)-gammln(a)-gammln(b)+a*log(x)+b*log(1.0-x)); if (x #define MAXIT 100 #define EPS 3.0e-7 #define FPMIN 1.0e-30 float betacf(float a, float b, float x) Used by betai: Evaluates continued fraction for incomplete beta function by modified Lentz’s method (§5.2). { void nrerror(char error_text[]);
228 Chapter 6. Special Functions int m,m2; float aa,c,d,del,h,qab,qam,qap; qab=a+b; These g's will be used in factors that occur 9ap=a+1.0: in the coefficients (6.4.6). qam=a-1.0; c=1.0; First step of Lentz's method d=1.0-qab*x/qap; if (fabs(d)FPMIN)d=FPMIN; d=1.0/d; hed; for (m=1;m0.99.In other words,1-A(tlv)is the Numerical Recipes 431985 significance level at which the hypothesis that the means are equal is disproved. The mathematical definition of the function is (outside T< North Software. A(tv)= /2B(2,) (6.4.7) Limiting values are A(0lw)=0 A(∞lw)=1 (6.4.8) A(t)is related to the incomplete beta function I(a,b)by A=1-1（25 (6.4.9)
228 Chapter 6. Special Functions Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). int m,m2; float aa,c,d,del,h,qab,qam,qap; qab=a+b; These q’s will be used in factors that occur qap=a+1.0; in the coefficients (6.4.6). qam=a-1.0; c=1.0; First step of Lentz’s method. d=1.0-qab*x/qap; if (fabs(d) MAXIT) nrerror("a or b too big, or MAXIT too small in betacf"); return h; } Student’s Distribution Probability Function Student’s distribution, denoted A(t|ν), is useful in several statistical contexts, notably in the test of whether two observed distributions have the same mean. A(t|ν) is the probability, for ν degrees of freedom, that a certain statistic t (measuring the observed difference of means) would be smaller than the observed value if the means were in fact the same. (See Chapter 14 for further details.) Two means are significantly different if, e.g., A(t|ν) > 0.99. In other words, 1 − A(t|ν) is the significance level at which the hypothesis that the means are equal is disproved. The mathematical definition of the function is A(t|ν) = 1 ν1/2B( 1 2 , ν 2 ) t −t 1 + x2 ν − ν+1 2 dx (6.4.7) Limiting values are A(0|ν)=0 A(∞|ν)=1 (6.4.8) A(t|ν) is related to the incomplete beta function Ix(a, b) by A(t|ν)=1 − I ν ν+t2 ν 2 , 1 2 (6.4.9)
6.4 Incomplete Beta Function 229 So,you can use(6.4.9)and the above routine betai to evaluate the function. F-Distribution Probability Function This function occurs in the statistical test of whether two observed samples have the same variance.A certain statistic F,essentially the ratio of the observed dispersion of the first sample to that of the second one,is calculated.(For further details,see Chapter 14.)The probability that F would be as large as it is if the first sample's underlying distribution actually has smaller variance than the second's is denoted (Fv1,v2),where vi and v2 are the number of degrees of freedom in the first and second samples,respectively.In other words,Q(Fv,v2)is the significance 虽 level at which the hypothesis"I has smaller variance than 2"can be rejected.A small numerical value implies a very significant rejection,in turn implying high confidence in the hypothesis"1 has variance greater or equal to 2." Q(Fv,v2)has the limiting values 袋 Q0,2)=1 Q(∞h,2)=0 (6.4.10) 2 Its relation to the incomplete beta function I(a.b)as evaluated by betai above is 、。 Press. Q(FlV1,v2)=I V2 V1 2+1F 2 (6.4.11) Programs SCIENTIFIC 6 Cumulative Binomial Probability Distribution Suppose an event occurs with probability p per trial.Then the probability P of its occurring k or more times in n trials is termed a cumulative binomial probability, and is related to the incomplete beta function I(a,b)as follows: 10.621 (1-p)n-j=I2(k,n-k+1) 43106 (6.4.12) E喜 Numerical Recipes (outside For n larger than a dozen or so,betai is a much better way to evaluate the sum in North Software. (6.4.12)than would be the straightforward sum with concurrent computation of the binomial coefficients.(For n smaller than a dozen,either method is acceptable. CITED REFERENCES AND FURTHER READING: Abramowitz,M.,and Stegun,I.A.1964.Handbook of Mathematical Functions,Applied Mathe- matics Series,Volume 55 (Washington:National Bureau of Standards;reprinted 1968 by Dover Publications,New York).Chapters 6 and 26. Pearson,E.,and Johnson,N.1968,Tables of the Incomplete Beta Function (Cambridge:Cam- bridge University Press)
6.4 Incomplete Beta Function 229 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). So, you can use (6.4.9) and the above routine betai to evaluate the function. F-Distribution Probability Function This function occurs in the statistical test of whether two observed samples have the same variance. A certain statistic F, essentially the ratio of the observed dispersion of the first sample to that of the second one, is calculated. (For further details, see Chapter 14.) The probability that F would be as large as it is if the first sample’s underlying distribution actually has smaller variance than the second’s is denoted Q(F|ν1, ν2), where ν1 and ν2 are the number of degrees of freedom in the first and second samples, respectively. In other words, Q(F|ν 1, ν2)is the significance level at which the hypothesis “1 has smaller variance than 2” can be rejected. A small numerical value implies a very significant rejection, in turn implying high confidence in the hypothesis “1 has variance greater or equal to 2.” Q(F|ν1, ν2) has the limiting values Q(0|ν1, ν2)=1 Q(∞|ν1, ν2)=0 (6.4.10) Its relation to the incomplete beta function Ix(a, b) as evaluated by betai above is Q(F|ν1, ν2) = I ν2 ν2+ν1F ν2 2 , ν1 2 (6.4.11) Cumulative Binomial Probability Distribution Suppose an event occurs with probability p per trial. Then the probability P of its occurring k or more times in n trials is termed a cumulative binomial probability, and is related to the incomplete beta function Ix(a, b) as follows: P ≡ n j=k n j pj (1 − p) n−j = Ip(k, n − k + 1) (6.4.12) For n larger than a dozen or so, betai is a much better way to evaluate the sum in (6.4.12) than would be the straightforward sum with concurrent computation of the binomial coefficients. (For n smaller than a dozen, either method is acceptable.) CITED REFERENCES AND FURTHER READING: Abramowitz, M., and Stegun, I.A. 1964, Handbook of Mathematical Functions, Applied Mathematics Series, Volume 55 (Washington: National Bureau of Standards; reprinted 1968 by Dover Publications, New York), Chapters 6 and 26. Pearson, E., and Johnson, N. 1968, Tables of the Incomplete Beta Function (Cambridge: Cambridge University Press).
230 Chapter 6.Special Functions 6.5 Bessel Functions of Integer Order This section and the next one present practical algorithms for computing various kinds of Bessel functions of integer order.In $6.7 we deal with fractional order.In fact,the more complicated routines for fractional order work fine for integer order too.For integer order,however,the routines in this section (and 86.6)are simpler and faster.Their only drawback is that they are limited by the precision of the underlying rational approximations.For full double precision,it is best to work with the routines for fractional order in 86.7. For any real v,the Bessel function J(z)can be defined by the series representation 菲 J(x) (-7x2) (6.5.1) 3 The series converges for all z,but it is not computationally very useful for >1. For v not an integer the Bessel function Y()is given by Y(z)= J(z)cos(vT)-Jv() (6.5.2) sin(wT)】 The right-hand side goes to the correct limiting value Yn()as v goes to some integer n,but this is also not computationally useful. For arguments x0 Numerica 10621 For >v,both Bessel functions look qualitatively like sine or cosine waves whose amplitude decays as x-1/2.The asymptotic forms forvare 43126 2 1 (outside (6.5.4) 2. 71 1 In the transition region wherev,the typical amplitudes of the Bessel functions are on the order 21/31 0.4473 四）~3亦~ v1/3 21/3 (6.5.5) 3VT方~-0.7748 1 Y(w) v1/3
230 Chapter 6. Special Functions Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). 6.5 Bessel Functions of Integer Order This section and the next one present practical algorithms for computing various kinds of Bessel functions of integer order. In §6.7 we deal with fractional order. In fact, the more complicated routines for fractional order work fine for integer order too. For integer order, however, the routines in this section (and §6.6) are simpler and faster. Their only drawback is that they are limited by the precision of the underlying rational approximations. For full double precision, it is best to work with the routines for fractional order in §6.7. For any real ν, the Bessel function Jν(x) can be defined by the series representation Jν (x) = 1 2 x ν ∞ k=0 (−1 4x2)k k!Γ(ν + k + 1) (6.5.1) The series converges for all x, but it is not computationally very useful for x 1. For ν not an integer the Bessel function Yν (x) is given by Yν(x) = Jν(x) cos(νπ) − J−ν(x) sin(νπ) (6.5.2) The right-hand side goes to the correct limiting value Y n(x) as ν goes to some integer n, but this is also not computationally useful. For arguments x 0 (6.5.3) For x>ν, both Bessel functions look qualitatively like sine or cosine waves whose amplitude decays as x−1/2. The asymptotic forms for x ν are Jν(x) ∼ 2 πx cos x − 1 2 νπ − 1 4 π Yν(x) ∼ 2 πx sin x − 1 2 νπ − 1 4 π (6.5.4) In the transition region where x ∼ ν, the typical amplitudes of the Bessel functions are on the order Jν (ν) ∼ 21/3 32/3Γ( 2 3 ) 1 ν1/3 ∼ 0.4473 ν1/3 Yν (ν) ∼ − 21/3 31/6Γ( 2 3 ) 1 ν1/3 ∼ −0.7748 ν1/3 (6.5.5)
