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),taking in the worst case O(Vmax(a,b))iterations.But 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 <math.h> 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<(a+1.0)/(a+b+2.0)) Use continued fraction directly. Numerical Recipes 10-521 43108 return bt*betacf(a,b,x)/a; else Use continued fraction after making the sym- return 1.0-bt*betacf(b,a,1.0-x)/b; metry transformation. (outside 膜 North Software. which utilizes the continued fraction evaluation routine #include <math.h> #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 machine￾readable 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), taking in the worst case O( max(a, b)) iterations. But 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 <math.h> 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 < 0.0 || 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 < (a+1.0)/(a+b+2.0)) Use continued fraction directly. return bt*betacf(a,b,x)/a; else Use continued fraction after making the sym￾return 1.0-bt*betacf(b,a,1.0-x)/b; metry transformation. } which utilizes the continued fraction evaluation routine #include <math.h> #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[]);