f6-4

6.4IncompleteBetaFunction,Student’sDistribution,F-Distribution,CumulativeBinomialDistribution 219Sample 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).return END CITED REFERENCES AND FURTHER READING: Stegun, I.A., and Zucker, R. 1974, Journal of Research of the National Bureau of Standards , vol. 78B, pp. 199–216; 1976, op. cit., vol. 80B, pp. 291–311. Amos D.E. 1980, ACM Transactions on Mathematical Software , vol. 6, pp. 365–377 [1]; also vol. 6, pp. 420–428. 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 5. Wrench J.W. 1952, Mathematical Tablesand Other Aids to Computation , vol. 6, p. 255. [2] 6.4 Incomplete Beta Function, Student’s Distribution, F-Distribution, CumulativeBinomial Distribution The incomplete beta function is defined by Ix(a, b)≡Bx(a, b) B(a, b)≡1 B(a, b)/integraldisplayx 0ta−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 ) Ifaandbare both rather greater than one, then Ix(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). The incomplete beta function has a series expansion Ix(a, b)=xa(1−x)b aB(a, b)/bracketleftBigg 1+∞/summationdisplay n=0B(a+1,n+1 ) B(a+b, n+1 )xn+1/bracketrightBigg , (6.4.4 ) butthisdoesnotprovetobeveryusefulinitsnumericalevaluation. (Note,however, that the beta functions in the coefficients can be evaluated for each value of nwith just the previousvalue and a few multiplies, using equations6.1.9 and 6.1.3.) The continued fraction representationproves to be much more useful, Ix(a, b)=xa(1−x)b aB(a, b)/bracketleftbigg1 1+d1 1+d2 1+···/bracketrightbigg (6.4.5 ) 220 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).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 .80.2.4.6.81incomplete beta function Ix(a,b) x Figure 6.4.1. The incomplete beta function Ix(a, b )forfive 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). 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(/radicalbig max(a, b))iterations. But for x> (a+1 )/(a+b+2 )we can justusethesymmetryrelation(6.4.3)toobtainanequivalentcomputationwherethe continued fraction will also converge rapidly. Hence we have FUNCTION betai(a,b,x) REAL betai,a,b,x C USES betacf,gammln Returns the incomplete beta function Ix(a,b). REAL bt,betacf,gammln if(x.lt.0..or.x.gt.1.)pause ’bad argument x in betai’ if(x.eq.0..or.x.eq.1.)then bt=0. else Factors in front of the continued fraction. bt=exp(gammln(a+b)-gammln(a)-gammln(b) * +a*log(x)+b*log(1.-x)) endif if(x.lt.(a+1.)/(a+b+2.))then Use continued fraction directly. 6.4IncompleteBetaFunction,Student’sDistribution,F-Distribution,CumulativeBinomialDistribution 221Sample 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).betai=bt*betacf(a,b,x)/a return else betai=1.-bt*betacf(b,a,1.-x)/b Use continued fraction after making the symme- try transformation. return endif END which utilizes the continued fraction evaluation routine FUNCTION betacf(a,b,x)INTEGER MAXIT REAL betacf,a,b,x,EPS,FPMIN PARAMETER (MAXIT=100,EPS=3.e-7,FPMIN=1.e-30) Used by betai : Evaluates continued fraction for incomplete beta function by modified Lentz’s method ( §5.2). INTEGER m,m2 REAL aa,c,d,del,h,qab,qam,qapqab=a+b These q’s will be used in factors that occur in the coefficients (6.4.6). qap=a+1. qam=a-1.c=1. First step of Lentz’s method. d=1.-qab*x/qap if(abs(d).lt.FPMIN)d=FPMIN d=1./dh=ddo 11m=1,MAXIT m2=2*m aa=m*(b-m)*x/((qam+m2)*(a+m2))d=1.+aa*d One step (the even one) of the recurrence. if(abs(d).lt.FPMIN)d=FPMIN c=1.+aa/cif(abs(c).lt.FPMIN)c=FPMINd=1./d h=h*d*c aa=-(a+m)*(qab+m)*x/((a+m2)*(qap+m2))d=1.+aa*d Next step of the recurrence (the odd one). if(abs(d).lt.FPMIN)d=FPMIN c=1.+aa/cif(abs(c).lt.FPMIN)c=FPMINd=1./d del=d*c h=h*delif(abs(del-1.).lt.EPS)goto 1 Are we done? enddo 11 pause ’a or b too big, or MAXIT too small in betacf’ 1 betacf=h return END Student’s DistributionProbability Function Student’s distribution, denoted A(t|ν), is useful in several statistical contexts, notablyinthetestofwhethertwoobserveddistributionshavethesamemean. 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 222 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).significantly different if, e.g., A(t|ν)>0.99. In other words, 1−A(t|ν)is the significance level at which the hypothesisthat the means are equalis disproved. The mathematical de finition of the function is A(t|ν)=1 ν1/2B(1 2,ν 2)/integraldisplayt −t/parenleftbigg 1+x2 ν/parenrightbigg−ν+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/parenleftbiggν 2,1 2/parenrightbigg (6.4.9 ) So, you can use (6.4.9)and the above routine betaito evaluate the function. F-DistributionProbability 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 probabilitythat Fwould be as largeas it is if the first sample’s underlying distribution actually has smallervariance than the second ’si s denoted Q(F|ν1,ν2), where ν1andν2are the number of degrees of freedom in the firstandsecondsamples,respectively. Inotherwords, Q(F|ν1,ν2)isthesigni ficance level at which the hypothesis “1 has smaller variance than 2 ”can be rejected. A small numerical value implies a very signi ficant 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 relationtothe incompletebetafunction Ix(a, b)as evaluatedby betaiaboveis Q(F|ν1,ν2)=I ν2 ν2+ν1F/parenleftbiggν2 2,ν1 2/parenrightbigg (6.4.11 ) CumulativeBinomialProbabilityDistribution Supposeaneventoccurswith probability ppertrial. Thenthe probability Pof itsoccurring kormoretimesin ntrialsistermeda cumulativebinomialprobability , and is related to the incomplete beta function Ix(a, b)as follows: P≡n/summationdisplay j=k/parenleftbiggn j/parenrightbigg pj(1−p)n−j=Ip(k,n−k+1 ) ( 6.4.12 ) 6.5BesselFunctionsofIntegerOrder 223Sample 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).Fornlarger than a dozen or so, betaiis a much better way to evaluate the sum in (6.4.12)than would be the straightforward sum with concurrent computationof thebinomialcoef ficients. (For nsmaller than a dozen,either methodis 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 byDover Publications, New York), Chapters 6 and 26. Pearson, E., and Johnson, N. 1968, Tablesof the Incomplete Beta Function (Cambridge: Cam- bridge University Press). 6.5 Bessel Functions of Integer Order Thissectionandthenextonepresentpracticalalgorithmsforcomputingvarious 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 underlyingrationalapproximations. Forfulldoubleprecision,itisbesttoworkwiththe routines for fractional order in §6.7. For any real ν, the Bessel function J ν(x)can be de fined by the series representation Jν(x)=/parenleftbigg1 2x/parenrightbiggν∞/summationdisplay k=0(−1 4x2)k k!Γ(ν+k+1 )(6.5.1 ) Theseries convergesfor all x, but it is not computationallyveryusefulfor x/greatermuch1. Forνnotan integer the Bessel function Yν(x)is given by Yν(x)=Jν(x)c o s ( νπ)−J−ν(x) sin(νπ)(6.5.2 ) Theright-handsidegoestothecorrectlimitingvalue Yn(x)asνgoestosomeinteger n, but this is also not computationally useful. For arguments x<ν, both Bessel functions look qualitatively like simple power laws, with the asymptotic forms for 0<x/lessmuchν Jν(x)∼1 Γ(ν+1 )/parenleftbigg1 2x/parenrightbiggν ν≥0 Y0(x)∼2 πln(x) Yν(x)∼−Γ(ν) π/parenleftbigg1 2x/parenrightbigg−ν ν> 0(6.5.3 )