Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C (523184), страница 57
Текст из файла (страница 57)
The continued fraction —just equation (6.2.6) rewritten — converges for all x > 0:n1 n+1 21−x···(6.3.4)En (x) = ex+ 1+ x+ 1+ x+We use it in its more rapidly converging even form,1·n2(n + 1)1−x···En (x) = ex+n− x+n+2− x+n+4−(6.3.5)The continued fraction only really converges fast enough to be useful for x >∼ 1.For 0 < x <1,wecanusetheseriesrepresentation∼En (x) =(−x)n−1[− ln x + ψ(n)] −(n − 1)!∞m=0m=n−1(−x)m(m − n + 1)m!(6.3.6)The quantity ψ(n) here is the digamma function, given for integer arguments byψ(1) = −γ,ψ(n) = −γ +n−1m=11m(6.3.7)6.3 Exponential Integrals223where γ = 0.5772156649 .
. . is Euler’s constant. We evaluate the expression (6.3.6)in order of ascending powers of x:xx2(−x)n−21−+−···+(1 − n) (2 − n) · 1 (3 − n)(1 · 2)(−1)(n − 2)!(−x)n+1(−x)n−1(−x)n[− ln x + ψ(n)] −++···+(n − 1)!1 · n!2 · (n + 1)!(6.3.8)En (x) = −The first square bracket is omitted when n = 1. This method of evaluation has theadvantage that for large n the series converges before reaching the term containingψ(n).
Accordingly, one needs an algorithm for evaluating ψ(n) only for small n,n <∼ 20 – 40. We use equation (6.3.7), although a table look-up would improveefficiency slightly.Amos [1] presents a careful discussion of the truncation error in evaluatingequation (6.3.8), and gives a fairly elaborate termination criterion.
We have foundthat simply stopping when the last term added is smaller than the required toleranceworks about as well.Two special cases have to be handled separately:e−xx1,En (0) =n−1E0 (x) =(6.3.9)n>1The routine expint allows fast evaluation of En (x) to any accuracy EPSwithin the reach of your machine’s word length for floating-point numbers. Theonly modification required for increased accuracy is to supply Euler’s constant withenough significant digits. Wrench [2] can provide you with the first 328 digits ifnecessary!#include <math.h>#define MAXIT 100#define EULER 0.5772156649#define FPMIN 1.0e-30#define EPS 1.0e-7Maximum allowed number of iterations.Euler’s constant γ.Close to smallest representable floating-point number.Desired relative error, not smaller than the machine precision.float expint(int n, float x)Evaluates the exponential integral En (x).{void nrerror(char error_text[]);int i,ii,nm1;float a,b,c,d,del,fact,h,psi,ans;nm1=n-1;if (n < 0 || x < 0.0 || (x==0.0 && (n==0 || n==1)))nrerror("bad arguments in expint");else {if (n == 0) ans=exp(-x)/x;Special case.else {if (x == 0.0) ans=1.0/nm1;Another special case.else {224Chapter 6.Special Functionsif (x > 1.0) {Lentz’s algorithm (§5.2).b=x+n;c=1.0/FPMIN;d=1.0/b;h=d;for (i=1;i<=MAXIT;i++) {a = -i*(nm1+i);b += 2.0;d=1.0/(a*d+b);Denominators cannot be zero.c=b+a/c;del=c*d;h *= del;if (fabs(del-1.0) < EPS) {ans=h*exp(-x);return ans;}}nrerror("continued fraction failed in expint");} else {Evaluate series.ans = (nm1!=0 ? 1.0/nm1 : -log(x)-EULER);Set first term.fact=1.0;for (i=1;i<=MAXIT;i++) {fact *= -x/i;if (i != nm1) del = -fact/(i-nm1);else {psi = -EULER;Compute ψ(n).for (ii=1;ii<=nm1;ii++) psi += 1.0/ii;del=fact*(-log(x)+psi);}ans += del;if (fabs(del) < fabs(ans)*EPS) return ans;}nrerror("series failed in expint");}}}}return ans;}A good algorithm for evaluating Ei is to use the power series for small x andthe asymptotic series for large x.
The power series isEi(x) = γ + ln x +x2x++···1 · 1! 2 · 2!(6.3.10)where γ is Euler’s constant. The asymptotic expansion isEi(x) ∼exx1+1!2!+ 2 +···xx(6.3.11)The lower limit for the use of the asymptotic expansion is approximately | ln EPS|,where EPS is the required relative error.6.3 Exponential Integrals#include <math.h>#define EULER 0.57721566#define MAXIT 100#define FPMIN 1.0e-30#define EPS 6.0e-8225Euler’s constant γ.Maximum number of iterations allowed.Close to smallest representable floating-point number.Relative error, or absolute error near the zero of Ei atx = 0.3725.float ei(float x)Computes the exponential integral Ei(x) for x > 0.{void nrerror(char error_text[]);int k;float fact,prev,sum,term;if (x <= 0.0) nrerror("Bad argument in ei");if (x < FPMIN) return log(x)+EULER;Special case: avoid failure of convergenceif (x <= -log(EPS)) {test because of underflow.sum=0.0;Use power series.fact=1.0;for (k=1;k<=MAXIT;k++) {fact *= x/k;term=fact/k;sum += term;if (term < EPS*sum) break;}if (k > MAXIT) nrerror("Series failed in ei");return sum+log(x)+EULER;} else {Use asymptotic series.sum=0.0;Start with second term.term=1.0;for (k=1;k<=MAXIT;k++) {prev=term;term *= k/x;if (term < EPS) break;Since final sum is greater than one, term itself approximates the relative error.if (term < prev) sum += term;Still converging: add new term.else {sum -= prev;Diverging: subtract previous term andbreak;exit.}}return exp(x)*(1.0+sum)/x;}}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]; alsovol. 6, pp. 420–428.Abramowitz, M., and Stegun, I.A. 1964, Handbook of Mathematical Functions, Applied Mathematics Series, Volume 55 (Washington: National Bureau of Standards; reprinted 1968 byDover Publications, New York), Chapter 5.Wrench J.W.
1952, Mathematical Tables and Other Aids to Computation, vol. 6, p. 255. [2]226Chapter 6.Special Functions1(0.5,5.0)incomplete beta function Ix (a,b).8(8.0,10.0)(1.0,3.0).6(0.5,0.5).4.2(5.0,0.5)00.2.4.6.81xFigure 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 related by reflection symmetry around the diagonal (cf. equation 6.4.3).6.4 Incomplete Beta Function, Student’sDistribution, F-Distribution, CumulativeBinomial DistributionThe incomplete beta function is defined by& x1Bx (a, b)≡Ix (a, b) ≡ta−1 (1 − t)b−1 dtB(a, b)B(a, b) 0(a, b > 0)(6.4.1)It has the limiting valuesI0 (a, b) = 0I1 (a, b) = 1(6.4.2)and the symmetry relationIx (a, b) = 1 − I1−x(b, a)(6.4.3)If a and b are 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 functionfor several pairs (a, b).6.4 Incomplete Beta FunctionThe incomplete beta function has a series expansion !∞xa(1 − x)bB(a + 1, n + 1) n+11+x,Ix (a, b) =aB(a, b)B(a + b, n + 1)n=0227(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 withjust 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 1 d1 d2···aB(a, b) 1+ 1+ 1+(6.4.5)where(a + m)(a + b + m)x(a + 2m)(a + 2m + 1)m(b − m)x=(a + 2m − 1)(a + 2m)d2m+1 = −d2m(6.4.6)This continued fractionconverges 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 canjust use the symmetry relation (6.4.3) to obtain an equivalent computation where thecontinued 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;elseFactors 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;elseUse continued fraction after making the symreturn 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-30float betacf(float a, float b, float x)Used by betai: Evaluates continued fraction for incomplete beta function by modified Lentz’smethod (§5.2).{void nrerror(char error_text[]);228Chapter 6.Special Functionsint m,m2;float aa,c,d,del,h,qab,qam,qap;qab=a+b;These q’s will be used in factors that occurqap=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;h=d;for (m=1;m<=MAXIT;m++) {m2=2*m;aa=m*(b-m)*x/((qam+m2)*(a+m2));d=1.0+aa*d;One step (the even one) of the recurrence.if (fabs(d) < FPMIN) d=FPMIN;c=1.0+aa/c;if (fabs(c) < FPMIN) c=FPMIN;d=1.0/d;h *= d*c;aa = -(a+m)*(qab+m)*x/((a+m2)*(qap+m2));d=1.0+aa*d;Next step of the recurrence (the odd one).if (fabs(d) < FPMIN) d=FPMIN;c=1.0+aa/c;if (fabs(c) < FPMIN) c=FPMIN;d=1.0/d;del=d*c;h *= del;if (fabs(del-1.0) < EPS) break;Are we done?}if (m > MAXIT) nrerror("a or b too big, or MAXIT too small in betacf");return h;}Student’s Distribution Probability FunctionStudent’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 theobserved difference of means) would be smaller than the observed value if themeans were in fact the same. (See Chapter 14 for further details.) Two means aresignificantly different if, e.g., A(t|ν) > 0.99. In other words, 1 − A(t|ν) is thesignificance level at which the hypothesis that the means are equal is disproved.The mathematical definition of the function isA(t|ν) =1ν 1/2B( 12 , ν2 )& t1+−tx2ν− ν+12dx(6.4.7)Limiting values areA(0|ν) = 0A(∞|ν) = 1A(t|ν) is related to the incomplete beta function Ix (a, b) byν 1,A(t|ν) = 1 − I νν+t22 2(6.4.8)(6.4.9)2296.4 Incomplete Beta FunctionSo, you can use (6.4.9) and the above routine betai to evaluate the function.F-Distribution Probability FunctionThis function occurs in the statistical test of whether two observed sampleshave the same variance.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
