292 Chapter 7.Random Numbers A gamma deviate has probability pa()dz of occurring with a value between x and x+dx,where Pa(z)dr =2a-le-* I(a)d x>0 (7.3.1) To generate deviates of(7.3.1)for small values of a,it is best to add up a exponentially distributed waiting times,i.e.,logarithms of uniform deviates.Since the sum of logarithms is the logarithm of the product,one really has only to generate the product of a uniform deviates,then take the log. 8 For larger values of a,the distribution(7.3.1)has a typically "bell-shaped" form,with a peak at x=a and a half-width of about va. We will be interested in several probability distributions with this same qual- 鱼a itative form.A useful comparison function in such cases is derived from the ICAL Lorentzian distribution po=(子) (7.3.2) RECIPES whose inverse indefinite integral is just the tangent function.It follows that the z-coordinate of an area-uniform random point under the comparison function Co f()=1+(z-x0/喝 (7.3.3) 9 for any constants ao,co,and xo,can be generated by the prescription x=ao tan(πU）+xo (7.3.4) 6 where U is a uniform deviate between 0 and 1.Thus,for some specific "bell-shaped" p()probability distribution,we need only find constants ao,co,o,with the product aoco(which determines the area)as small as possible,such that(7.3.3)is everywhere greater than p(x). Ahrens has done this for the gamma distribution,yielding the following o￥●N Numerica 10.621 algorithm (as described in Knuth [1]): 4312 #include <math.h> Recipes float gamdev(int ia,long *idum) Returns a deviate distributed as a gamma distribution of integer order ia,i.e.,a waiting time North Software. to the iath event in a Poisson process of unit mean,using ran1(idum)as the source of uniform deviates. float ran1(long *idum); void nrerror(char error_text[]); int j; float am,e,s,v1,v2,x,y; if (ia 1)nrerror("Error in routine gamdev"); if (ia 6){ Use direct method,adding waiting X=1.0; times. for(j=1;j<=ia;j++)x米=ran1(1dum）: x=-1og(x); else Use rejection method.292 Chapter 7. Random Numbers 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). A gamma deviate has probability pa(x)dx of occurring with a value between x and x + dx, where pa(x)dx = xa−1e−x Γ(a) dx x > 0 (7.3.1) To generate deviates of (7.3.1) for small values of a, it is best to add up a exponentially distributed waiting times, i.e., logarithms of uniform deviates. Since the sum of logarithms is the logarithm of the product, one really has only to generate the product of a uniform deviates, then take the log. For larger values of a, the distribution (7.3.1) has a typically “bell-shaped” form, with a peak at x = a and a half-width of about √a. We will be interested in several probability distributions with this same qual￾itative form. A useful comparison function in such cases is derived from the Lorentzian distribution p(y)dy = 1 π
1 1 + y2 dy (7.3.2) whose inverse indefinite integral is just the tangent function. It follows that the x-coordinate of an area-uniform random point under the comparison function f(x) = c0 1+(x − x0)2/a2 0 (7.3.3) for any constants a0, c0, and x0, can be generated by the prescription x = a0 tan(πU) + x0 (7.3.4) where U is a uniform deviate between 0 and 1. Thus, for some specific “bell-shaped” p(x) probability distribution, we need only find constants a 0, c0, x0, with the product a0c0 (which determines the area) as small as possible, such that (7.3.3) is everywhere greater than p(x). Ahrens has done this for the gamma distribution, yielding the following algorithm (as described in Knuth [1]): #include <math.h> float gamdev(int ia, long *idum) Returns a deviate distributed as a gamma distribution of integer order ia, i.e., a waiting time to the iath event in a Poisson process of unit mean, using ran1(idum) as the source of uniform deviates. { float ran1(long *idum); void nrerror(char error_text[]); int j; float am,e,s,v1,v2,x,y; if (ia < 1) nrerror("Error in routine gamdev"); if (ia < 6) { Use direct method, adding waiting x=1.0; times. for (j=1;j<=ia;j++) x *= ran1(idum); x = -log(x); } else { Use rejection method.