f6-2

6.2IncompleteGammaFunction 209Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).If your problem requires a series of related binomial coefficients, a good idea is to use recurrence relations, for example /parenleftbiggn+1 k/parenrightbigg =n+1 n−k+1/parenleftbiggn k/parenrightbigg =/parenleftbiggn k/parenrightbigg +/parenleftbiggn k−1/parenrightbigg /parenleftbiggn k+1/parenrightbigg =n−k k+1/parenleftbiggn k/parenrightbigg (6.1.7 ) Finally, turning away from the combinatorial functions with integer valued arguments, we come to the beta function, B(z,w)=B(w, z)=/integraldisplay1 0tz−1(1−t)w−1dt (6.1.8 ) which is related to the gamma function by B(z,w)=Γ(z)Γ(w) Γ(z+w)(6.1.9 ) hence FUNCTION beta(z,w) REAL beta,w,z C USES gammln Returns the value of the beta function B(z, w). REAL gammlnbeta=exp(gammln(z)+gammln(w)-gammln(z+w))return END 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 byDover Publications, New York), Chapter 6. Lanczos, C. 1964, SIAM Journal on Numerical Analysis , ser. B, vol. 1, pp. 86–96. [1] 6.2 Incomplete Gamma Function, Error Function, Chi-Square Probability Function, Cumulative Poisson Function The incomplete gamma function is defined by P(a, x)≡γ(a, x) Γ(a)≡1 Γ(a)/integraldisplayx 0e−tta−1dt (a>0) ( 6.2.1 ) 210 Chapter6. SpecialFunctionsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).024681 0 1 2 1 40.2.4.6.81.0 a = 3.01.00.5incomplete gamma function P(a,x) xa = 10 Figure 6.2.1. The incomplete gamma function P(a, x)for four values of a. It has the limiting values P(a,0) = 0 and P(a,∞)=1 ( 6.2.2 ) Theincompletegammafunction P(a, x)ismonotonicand(for agreaterthanoneor so) rises from “near-zero ”to“near-unity ”in a range of xcentered on about a−1, and of width about√a(see Figure 6.2.1). The complement of P(a, x)is also confusingly called an incomplete gamma function, Q(a, x)≡1−P(a, x)≡Γ(a, x) Γ(a)≡1 Γ(a)/integraldisplay∞ xe−tta−1dt (a>0) (6.2.3 ) It has the limiting values Q(a,0) = 1 and Q(a,∞)=0 ( 6.2.4 ) The notations P(a, x),γ(a, x), and Γ(a, x)are standard; the notation Q(a, x)is specific to this book. There is a series development for γ(a, x)as follows: γ(a, x)=e−xxa∞/summationdisplay n=0Γ(a) Γ(a+1+ n)xn(6.2.5 ) One does not actually need to compute a new Γ(a+1+ n)for each n; one rather uses equation (6.1.3) and the previous coef ficient. 6.2IncompleteGammaFunction 211Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).A continued fraction development for Γ(a, x)is Γ(a, x)=e−xxa/parenleftbigg1 x+1−a 1+1 x+2−a 1+2 x+···/parenrightbigg (x>0) (6.2.6 ) It is computationally better to use the even part of (6.2.6), which converges twice as fast (see §5.2): Γ(a, x)=e−xxa/parenleftbigg1 x+1−a−1·(1−a) x+3−a−2·(2−a) x+5−a−···/parenrightbigg (x>0) (6.2.7 ) It turns out that (6.2.5) converges rapidly for xless than about a+1, while (6.2.6)or(6.2.7)convergesrapidlyfor xgreaterthanabout a+1. Intheserespective regimes each requires at most a few times√aterms to converge, and this many only near x=a, where the incomplete gamma functions are varying most rapidly. Thus (6.2.5) and (6.2.7) together allow evaluation of the function for all positive aandx. An extra dividend is that we never need compute a function value near zerobysubtractingtwonearlyequalnumbers. Thehigher-levelfunctionsthatreturn P(a, x)andQ(a, x)are FUNCTION gammp(a,x) REAL a,gammp,x C USES gcf,gser Returns the incomplete gamma function P(a, x). REAL gammcf,gamser,gln if(x.lt.0..or.a.le.0.)pause ’bad arguments in gammp’if(x.lt.a+1.)then Use the series representation. call gser(gamser,a,x,gln) gammp=gamser else Use the continued fraction representation call gcf(gammcf,a,x,gln) gammp=1.-gammcf and take its complement. endifreturn END FUNCTION gammq(a,x) REAL a,gammq,x C USES gcf,gser Returns the incomplete gamma function Q(a, x)≡1−P(a, x). REAL gammcf,gamser,gln if(x.lt.0..or.a.le.0.)pause ’bad arguments in gammq’if(x.lt.a+1.)then Use the series representation call gser(gamser,a,x,gln) gammq=1.-gamser and take its complement. else Use the continued fraction representation. call gcf(gammcf,a,x,gln) gammq=gammcf endifreturn END 212 Chapter6. SpecialFunctionsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).The argument glnis returned by both the series and continued fraction procedures containing the value ln Γ(a); the reason for this is so that it is available to you if you want to modifythe above two procedures to give γ(a, x)andΓ(a, x), in addition to P(a, x)andQ(a, x)(cf. equations 6.2.1 and 6.2.3). The procedures gserandgcfwhich implement (6.2.5)and (6.2.7)are SUBROUTINE gser(gamser,a,x,gln) INTEGER ITMAX REAL a,gamser,gln,x,EPSPARAMETER (ITMAX=100,EPS=3.e-7) C USES gammln Returns the incomplete gamma function P(a, x)evaluated by its series representation as gamser . Also returns ln Γ( a)asgln . INTEGER n REAL ap,del,sum,gammln gln=gammln(a)if(x.le.0.)then if(x.lt.0.)pause ’x < 0 in gser’ gamser=0.return endif ap=a sum=1./adel=sumdo 11n=1,ITMAX ap=ap+1. del=del*x/apsum=sum+del if(abs(del).lt.abs(sum)*EPS)goto 1 enddo 11 pause ’a too large, ITMAX too small in gser’ 1 gamser=sum*exp(-x+a*log(x)-gln) return END SUBROUTINE gcf(gammcf,a,x,gln) INTEGER ITMAXREAL a,gammcf,gln,x,EPS,FPMIN PARAMETER (ITMAX=100,EPS=3.e-7,FPMIN=1.e-30) C USES gammln Returns the incomplete gamma function Q(a, x)evaluated by its continued fraction repre- sentation as gammcf . Also returns ln Γ( a)asgln . Parameters: ITMAX is the maximum allowed number of iterations; EPS is the relative accu- racy; FPMIN is a number near the smallest representable floating-point number. INTEGER i REAL an,b,c,d,del,h,gammln gln=gammln(a)b=x+1.-a Set up for evaluating continued fraction by modified Lentz’s method ( §5.2) with b 0=0. c=1./FPMIN d=1./b h=ddo 11i=1,ITMAX Iterate to convergence. an=-i*(i-a) b=b+2.d=an*d+bif(abs(d).lt.FPMIN)d=FPMIN c=b+an/c if(abs(c).lt.FPMIN)c=FPMINd=1./d del=d*c 6.2IncompleteGamma Function 213Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).h=h*del if(abs(del-1.).lt.EPS)goto 1 enddo 11 pause ’a too large, ITMAX too small in gcf’ 1 gammcf=exp(-x+a*log(x)-gln)*h Put factors in front. return END ErrorFunction The error function and complementary error function are special cases of the incomplete gamma function, and are obtained moderately ef ficiently by the above procedures. Their de finitions are erf(x)=2√π/integraldisplayx 0e−t2dt (6.2.8 ) and erfc(x)≡1−erf(x)=2√π/integraldisplay∞ xe−t2dt (6.2.9 ) The functions have the following limiting values and symmetries: erf(0) = 0 erf(∞)=1 erf(−x)=−erf(x)(6.2.10) erfc(0) = 1 erfc(∞)=0 erfc(−x)=2 −erfc(x)(6.2.11) They are related to the incomplete gamma functions by erf(x)=P/parenleftbigg1 2,x2/parenrightbigg (x≥0) ( 6.2.12 ) and erfc(x)=Q/parenleftbigg1 2,x2/parenrightbigg (x≥0) ( 6.2.13 ) Hence we have FUNCTION erf(x) REAL erf,x C USES gammp Returns the error function erf (x). REAL gammpif(x.lt.0.)then erf=-gammp(.5,x**2) else erf=gammp(.5,x**2) endif return END 214 Chapter6. SpecialFunctionsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).FUNCTION erfc(x) REAL erfc,x C USES gammp,gammq Returns the complementary error function erfc (x). REAL gammp,gammq if(x.lt.0.)then erfc=1.+gammp(.5,x**2) else erfc=gammq(.5,x**2) endif returnEND If you care to do so, you can easily remedy the minor inef ficiency in erfand erfc, namely that Γ(0.5) =√πis computed unnecessarily when gammporgammq is called. Before you do that, however, you might wish to consider the following routine,based onChebyshev fitting to an inspiredguess as to the functionalform: FUNCTION erfcc(x) REAL erfcc,x Returns the complementary error function erfc (x)with fractional error everywhere less than 1.2×10−7. REAL t,zz=abs(x) t=1./(1.+0.5*z) erfcc=t*exp(-z*z-1.26551223+t*(1.00002368+t*(.37409196+ * t*(.09678418+t*(-.18628806+t*(.27886807+t*(-1.13520398+* t*(1.48851587+t*(-.82215223+t*.17087277))))))))) if (x.lt.0.) erfcc=2.-erfcc returnEND There are also some functions of twovariables that are special cases of the incomplete gamma function: CumulativePoisson ProbabilityFunction Px(<k), for positive xand integer k≥1, denotes the cumulative Poisson probability function. It is de fined as the probability that the number of Poisson randomeventsoccurringwillbebetween0and k−1inclusive,iftheexpectedmean number is x. It has the limiting values Px(<1) = e−xPx(<∞)=1 ( 6.2.14 ) Its relation to the incomplete gamma function is simply Px(<k)=Q(k,x)=gammq (k,x)( 6.2.15 ) 6.3ExponentialIntegrals 215Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Chi-Square ProbabilityFunction P(χ2|ν)is definedas the probabilitythat the observedchi-squarefor a correct model should be less than a value χ2. (We will discuss the use of this function in Chapter15.) Itscomplement Q(χ2|ν)istheprobabilitythattheobservedchi-square will exceed the value χ2by chance evenfor a correct model. In both cases νis an integer,the numberof degreesof freedom. Thefunctionshavethe limitingvalues P(0|ν)=0 P(∞|ν)=1 ( 6.2.16 ) Q(0|ν)=1 Q(∞|ν)=0 ( 6.2.17 ) and the following relation to the incomplete gamma functions, P(χ2|ν)=P/parenleftbiggν 2,χ2 2/parenrightbigg =gammp/parenleftbiggν 2,χ2 2/parenrightbigg (6.2.18 ) Q(χ2|ν)=Q/parenleftbiggν 2,χ2 2/parenrightbigg =gammq/parenleftbiggν 2,χ2 2/parenrightbigg (6.2.19 ) 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 byDover Publications, New York), Chapters 6, 7, and 26. Pearson, K. (ed.) 1951, Tables of the Incomplete Gamma Function (Cambridge: Cambridge University Press). 6.3 Exponential Integrals The standard de finition of the exponential integral is En(x)=/integraldisplay∞ 1e−xt tndt, x > 0,n =0,1,... (6.3.1 ) The function de fined by the principal value of the integral Ei(x)=−/integraldisplay∞ −xe−t tdt=/integraldisplayx −∞et tdt, x > 0( 6.3.2 ) is also called an exponential integral. Note that Ei(−x)is related to −E1(x)by analytic continuation. The function En(x)is a special case of the incomplete gamma function En(x)=xn−1Γ(1−n, x)( 6.3.3 )