290 Chapter 7.Random Numbers while (rsq >1.0 II rsq =0.0); and if they are not,try again. fac=sgrt(-2.0*log(rsq)/rsq); Now make the Box-Muller transformation to get two normal deviates.Return one and save the other for next time. gset=v1*fac; 18et=1; Set flag. return v2*fac; else We have an extra deviate handy, iset=0; so unset the flag, return gset; and return it. See Devroye [1]and Bratley [2]for many additional algorithms. 83g nted for CITED REFERENCES AND FURTHER READING: Devroye,L.1986,Non-Uniform Random Variate Generation (New York:Springer-Verlag),89.1. [1] Bratley.P.,Fox,B.L.,and Schrage,E.L.1983,A Guide to Simulation (New York:Springer- Verlag).[2] Knuth,D.E.1981,Seminumerical Algorithms,2nd ed.,vol.2 of The Art of Computer Programming (Reading,MA:Addison-Wesley),pp.116ff. 为 Press. 7.3 Rejection Method:Gamma,Poisson, Binomial Deviates IENTIFIC The rejection method is a powerful,general technique for generating random deviates whose distribution function p(x)da(probability ofa value occurring between x and z+dz)is known and computable.The rejection method does not require that the cumulative distribution function [indefinite integral of p(z)]be readily computable,much less the inverse of that function-which was required for the transformation method in the previous section. 10621 The rejection method is based on a simple geometrical argument: Numerica Draw a graph of the probability distribution p()that you wish to generate,so 43106 that the area under the curve in any range ofx corresponds to the desired probability (outside Recipes of generating an x in that range.If we had some way of choosing a random point in two dimensions,with uniform probability in the area under your curve,then the x 腿 value of that random point would have the desired distribution. North Now,on the same graph,draw any other curve f()which has finite (not infinite)area and lies everywhere above your original probability distribution.(This is always possible,because your original curve encloses only unit area,by definition of probability.)We will call this f(r)the comparison function.Imagine now that you have some way of choosing a random point in two dimensions that is uniform in the area under the comparison function.Whenever that point lies outside the area under the original probability distribution,we will reject it and choose another random point.Whenever it lies inside the area under the original probability distribution,we will accept it.It should be obvious that the accepted points are uniform in the accepted area,so that their values have the desired distribution.It
290 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 machinereadable 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). } while (rsq >= 1.0 || rsq == 0.0); and if they are not, try again. fac=sqrt(-2.0*log(rsq)/rsq); Now make the Box-Muller transformation to get two normal deviates. Return one and save the other for next time. gset=v1*fac; iset=1; Set flag. return v2*fac; } else { We have an extra deviate handy, iset=0; so unset the flag, return gset; and return it. } } See Devroye [1] and Bratley [2] for many additional algorithms. CITED REFERENCES AND FURTHER READING: Devroye, L. 1986, Non-Uniform Random Variate Generation (New York: Springer-Verlag), §9.1. [1] Bratley, P., Fox, B.L., and Schrage, E.L. 1983, A Guide to Simulation (New York: SpringerVerlag). [2] Knuth, D.E. 1981, Seminumerical Algorithms, 2nd ed., vol. 2 of The Art of Computer Programming (Reading, MA: Addison-Wesley), pp. 116ff. 7.3 Rejection Method: Gamma, Poisson, Binomial Deviates The rejection method is a powerful, general technique for generating random deviates whose distribution function p(x)dx (probability of a value occurring between x and x + dx) is known and computable. The rejection method does not require that the cumulative distribution function [indefinite integral of p(x)] be readily computable, much less the inverse of that function — which was required for the transformation method in the previous section. The rejection method is based on a simple geometrical argument: Draw a graph of the probability distribution p(x) that you wish to generate, so that the area under the curve in any range of x corresponds to the desired probability of generating an x in that range. If we had some way of choosing a random point in two dimensions, with uniform probability in the area under your curve, then the x value of that random point would have the desired distribution. Now, on the same graph, draw any other curve f(x) which has finite (not infinite) area and lies everywhere above your original probability distribution. (This is always possible, because your original curve encloses only unit area, by definition of probability.) We will call this f(x) the comparison function. Imagine now that you have some way of choosing a random point in two dimensions that is uniform in the area under the comparison function. Whenever that point lies outside the area under the original probability distribution, we will reject it and choose another random point. Whenever it lies inside the area under the original probability distribution, we will accept it. It should be obvious that the accepted points are uniform in the accepted area, so that their x values have the desired distribution. It
7.3 Rejection Method:Gamma.Poisson,Binomial Deviates 291 A first random deviate in ∫6x)k reject.xo f(x) -f(xo) acceptxo second random p(x) deviate in 0 83 Figure 7.3.1.Rejection method for generating a random deviate z from a known probability distribution 鱼 p()that is everywhere less than some other function f().The transformation method is first used to generate a random deviate z of the distribution f (compare Figure 7.2.1).A second uniform deviate is used to decide whether to accept or reject that If it is rejected,a new deviate of f is found;and so on The ratio of accepted to rejected points is the ratio of the area under p to the area between p and f. should also be obvious that the fraction of points rejected just depends on the ratio of the area of the comparison function to the area of the probability distribution 令 function,not on the details of shape of either function.For example,a comparison function whose area is less than 2 will reject fewer than half the points,even if it Press. approximates the probability function very badly at some values of,e.g.,remains finite in some region where p(z)is zero. It remains only to suggest how to choose a uniform random point in two dimensions under the comparison function f(z).A variant of the transformation method (87.2)does nicely:Be sure to have chosen a comparison function whose OF SCIENTIFIC indefinite integral is known analytically,and is also analytically invertible to give z as a function of"area under the comparison function to the left of z."Now pick a uniform deviate between 0 and A,where A is the total area under f(x),and use it to get a corresponding z.Then pick a uniform deviate between 0 and f(x)as the y value for the two-dimensional point.You should be able to convince yourself that the point(,y)is uniformly distributed in the area under the comparison function f(). An equivalent procedure is to pick the second uniform deviate between zero Numerical Recipes 10-521 and one,and accept or reject according to whether it is respectively less than or greater than the ratio p()/f(x). 4310 So,to summarize,the rejection method for some given p(z)requires that one g find,once and for all,some reasonably good comparison function f(x).Thereafter, (outside each deviate generated requires two uniform random deviates,one evaluation of f (to get the coordinate y),and one evaluation of p(to decide whether to accept or reject North Software. the point y).Figure 7.3.1 illustrates the procedure.Then,of course,this procedure must be repeated,on the average,A times before the final deviate is obtained. Gamma Distribution The gamma distribution of integer order a >0 is the waiting time to the ath event in a Poisson random process of unit mean.For example,when a=1,it is just the exponential distribution of 87.2,the waiting time to the first event
7.3 Rejection Method: Gamma, Poisson, Binomial Deviates 291 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 machinereadable 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). second random deviate in 0 f(x0) reject x0 accept x0 first random deviate in A 0 f(x) 0 f(x)dx x p(x) x0 ⌠ ⌡ Figure 7.3.1. Rejection method for generating a random deviate x from a known probability distribution p(x) that is everywhere less than some other function f(x). The transformation method is first used to generate a random deviate x of the distribution f (compare Figure 7.2.1). A second uniform deviate is used to decide whether to accept or reject that x. If it is rejected, a new deviate of f is found; and so on. The ratio of accepted to rejected points is the ratio of the area under p to the area between p and f. should also be obvious that the fraction of points rejected just depends on the ratio of the area of the comparison function to the area of the probability distribution function, not on the details of shape of either function. For example, a comparison function whose area is less than 2 will reject fewer than half the points, even if it approximates the probability function very badly at some values of x, e.g., remains finite in some region where p(x) is zero. It remains only to suggest how to choose a uniform random point in two dimensions under the comparison function f(x). A variant of the transformation method (§7.2) does nicely: Be sure to have chosen a comparison function whose indefinite integral is known analytically, and is also analytically invertible to give x as a function of “area under the comparison function to the left of x.” Now pick a uniform deviate between 0 and A, where A is the total area under f(x), and use it to get a corresponding x. Then pick a uniform deviate between 0 and f(x) as the y value for the two-dimensional point. You should be able to convince yourself that the point (x, y) is uniformly distributed in the area under the comparison function f(x). An equivalent procedure is to pick the second uniform deviate between zero and one, and accept or reject according to whether it is respectively less than or greater than the ratio p(x)/f(x). So, to summarize, the rejection method for some given p(x) requires that one find, once and for all, some reasonably good comparison function f(x). Thereafter, each deviate generated requires two uniform random deviates, one evaluation of f (to get the coordinate y), and one evaluation of p (to decide whether to accept or reject the point x, y). Figure 7.3.1 illustrates the procedure. Then, of course, this procedure must be repeated, on the average, A times before the final deviate is obtained. Gamma Distribution The gamma distribution of integer order a > 0 is the waiting time to the ath event in a Poisson random process of unit mean. For example, when a = 1, it is just the exponential distribution of §7.2, the waiting time to the first event
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 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 machinereadable 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 qualitative 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 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.
7.3 Rejection Method:Gamma,Poisson,Binomial Deviates 293 do do do These four lines generate the tan- vi=ran1(idum); gent of a random angle,i.e.,they v2=2.0*ran1(1dum)-1.0: are equivalent to Jh11e(v1*v1+v2*v2>1.0); y tan(*ran1(idum)). y=v2/v1; am=1a-1; s=sqrt(2.0*am+1.0); x=s*y+am; We decide whether to reject x: Jwh11e(xe); Reject on basis of a second uniform deviate. return x; Cam Poisson Deviates The Poisson distribution is conceptually related to the gamma distribution.It gives the probability of a certain integer number m of unit rate Poisson random 令 events occurring in a given interval of time z,while the gamma distribution was the probability of waiting time between x and x+dz to the mth event.Note that m takes 未 Press. on only integer values >0,so that the Poisson distribution,viewed as a continuous distribution function p(m)dm,is zero everywhere except where m is an integer >0.At such places,it is infinite,such that the integrated probability over a region containing the integer is some finite number.The total probability at an integer j is j+e Pr(m)dm= xle-x Prob(j)= (7.3.5) 6 At first sight this might seem an unlikely candidate distribution for the rejection method,since no continuous comparison function can be larger than the infinitely tall,but infinitely narrow,Dirac delta fimnctions in p(m).However,there is a trick that we can do:Spread the finite area in the spike at j uniformly into the interval 、3p5 between j and j+1.This defines a continuous distribution g=(m)dm given by Numerica 10.621 431 4r(m)dm xmle-x lmjdm (7.3.6) Recipes where [m]represents the largest integer less than m.If we now use the rejection North method to generate a (noninteger)deviate from (7.3.6),and then take the integer part of that deviate,it will be as if drawn from the desired distribution(7.3.5).(See 、客 Figure 7.3.2.)This trick is general for any integer-valued probability distribution. For large enough,the distribution(7.3.6)is qualitatively bell-shaped(albeit with a bell made out of small,square steps),and we can use the same kind of Lorentzian comparison function as was already used above.For small we can generate independent exponential deviates(waiting times between events);when the sum of these first exceeds z,then the number of events that would have occurred in waiting time z becomes known and is one less than the number of terms in the sum. These ideas produce the following routine:
7.3 Rejection Method: Gamma, Poisson, Binomial Deviates 293 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 machinereadable 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). do { do { do { These four lines generate the tangent of a random angle, i.e., they are equivalent to v1=ran1(idum); v2=2.0*ran1(idum)-1.0; } while (v1*v1+v2*v2 > 1.0); y = tan(π * ran1(idum)). y=v2/v1; am=ia-1; s=sqrt(2.0*am+1.0); x=s*y+am; We decide whether to reject x: } while (x e); Reject on basis of a second uniform } deviate. return x; } Poisson Deviates The Poisson distribution is conceptually related to the gamma distribution. It gives the probability of a certain integer number m of unit rate Poisson random events occurring in a given interval of time x, while the gamma distribution was the probability of waiting time between x and x+dx to the mth event. Note that m takes on only integer values ≥ 0, so that the Poisson distribution, viewed as a continuous distribution function px(m)dm, is zero everywhere except where m is an integer ≥ 0. At such places, it is infinite, such that the integrated probability over a region containing the integer is some finite number. The total probability at an integer j is Prob(j) = j+ j− px(m)dm = xj e−x j! (7.3.5) At first sight this might seem an unlikely candidate distribution for the rejection method, since no continuous comparison function can be larger than the infinitely tall, but infinitely narrow, Dirac delta functions in px(m). However, there is a trick that we can do: Spread the finite area in the spike at j uniformly into the interval between j and j + 1. This defines a continuous distribution qx(m)dm given by qx(m)dm = x[m] e−x [m]! dm (7.3.6) where [m] represents the largest integer less than m. If we now use the rejection method to generate a (noninteger) deviate from (7.3.6), and then take the integer part of that deviate, it will be as if drawn from the desired distribution (7.3.5). (See Figure 7.3.2.) This trick is general for any integer-valued probability distribution. For x large enough, the distribution (7.3.6) is qualitatively bell-shaped (albeit with a bell made out of small, square steps), and we can use the same kind of Lorentzian comparison function as was already used above. For small x, we can generate independent exponential deviates (waiting times between events); when the sum of these first exceeds x, then the number of events that would have occurred in waiting time x becomes known and is one less than the number of terms in the sum. These ideas produce the following routine:
294 Chapter 7.Random Numbers 0 reject http://www.nr Permission is read able files .com or call 1-800-872- (including this one) granted fori internet 7423 1988-1992 by Cambridge University Press. from NUMERICAL RECIPES IN Figure 7.3.2.Rejection method as applied to an integer-valued distribution.The method is performed (North America to any server computer. t users to make one paper THE on the step function shown as a dashed line,yielding a real-valued deviate.This deviate is rounded down to the next lower integer,which is output. 是 ART #include #define PI3.141592654 st st Programs for their float poidev(float xm,long *idum) Returns as a floating-point number an integer value that is a random deviate drawn from a Poisson distribution of mean xm,using ran1(idum)as a source of uniform random deviates. to dir float gammln(float xx); float ran1(long *idum); ectcustser static float sq,alxm,g,oldm=(-1.0); oldm is a flag for whether xm has changed 18881920 OF SCIENTIFIC COMPUTING(ISBN float em,t,yi since last call. 1f(xm<12.0){ Use direct method. Further reproduction. Numerical Recipes 12-521 if(xm！=oldm){ oldm=xm; If xm is new,compute the exponential. 43108 g=exp(-xm) 2 em=-1: t=1.0; (outside do Instead of adding exponential deviates it is equiv- Software. ++em; alent to multiply uniform deviates.We never North t *=ran1(idum); actually have to take the log.merely com- while (t g); pare to the pre-computed exponential else Use rejection method. 1f(xm！=oldm)[ If xm has changed since the last call,then pre- visit website machine oldm=xm; compute some functions that occur below. sq=sqrt(2.0*xm); alxm=log(xm); g=xm*alxm-gammln(xm+1.0); The function gammin is the natural log of the gamma function,as given in $6.1. do do y is a deviate from a Lorentzian comparison func- y=tan(PI*ran1(idum)); tion
294 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 machinereadable 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). 0 12345 accept reject in 1 Figure 7.3.2. Rejection method as applied to an integer-valued distribution. The method is performed on the step function shown as a dashed line, yielding a real-valued deviate. This deviate is rounded down to the next lower integer, which is output. #include #define PI 3.141592654 float poidev(float xm, long *idum) Returns as a floating-point number an integer value that is a random deviate drawn from a Poisson distribution of mean xm, using ran1(idum) as a source of uniform random deviates. { float gammln(float xx); float ran1(long *idum); static float sq,alxm,g,oldm=(-1.0); oldm is a flag for whether xm has changed float em,t,y; since last call. if (xm g); } else { Use rejection method. if (xm != oldm) { If xm has changed since the last call, then preoldm=xm; compute some functions that occur below. sq=sqrt(2.0*xm); alxm=log(xm); g=xm*alxm-gammln(xm+1.0); The function gammln is the natural log of the gamma function, as given in §6.1. } do { do { y is a deviate from a Lorentzian comparison funcy=tan(PI*ran1(idum)); tion
7.3 Rejection Method:Gamma,Poisson,Binomial Deviates 295 em=sq*y+xm; em is y,shifted and scaled. ]whi1e(emt); return em; Binomial Deviates 83 鱼 granted for 18881892 If an event occurs with probability g,and we make n trials,then the number of 11800 times m that it occurs has the binomial distribution. from NUMERICAL RECIPESI ]+6 Pn.(m)dm (1-q)n-j (7.3.7) 芝是%6sadN (Nort server The binomial distribution is integer valued,with m taking on possible values America THE from 0 to n.It depends on two parameters,n and g,so is correspondingly a ART bit harder to implement than our previous examples.Nevertheless,the techniques already illustrated are sufficiently powerful to do the job: 9 Program #include #def1nePI3.141592654 float bnldev(float pp,int n,long *idum) Returns as a floating-point number an integer value that is a random deviate drawn from a binomial distribution of n trials each of probability pp.using ran1(idum)as a source of uniform random deviates. OF SCIENTIFIC COMPUTING(ISBN 1920 float gammln(float xx); float ran1(long *idum); int j; Numerical Recipes 2-521 static int nold=(-1); float am,em,g,angle,p,bnl,sq,t,y; 43108 static float pold=(-1.0),pc,plog,pclog,en,oldg; p=(pp<=0.5？pp:1.0-pp); The binomial distribution is invariant under changing pp to 1-pp,if we also change the (outside answer to n minus itself;we'll remember to do this below. North Software. am=n*p; This is the mean of the deviate to be produced. 1f(n<25) Use the direct method while n is not too large 6 bnl=0.0; This can require up to 25 calls to ran1. for (j=1;j<=n;j++) if (ran1(idum)<p)++bnl; 】e1seif(am<1.0)[ If fewer than one event is expected out of 25 g=exp(-am) or more trials,then the distribution is quite t=1.0: accurately Poisson.Use direct Poisson method. for (j=0;j<=n;j++){ t *ran1(idum); if (t g)break; bnl=(j <n j:n); else Use the rejection method
7.3 Rejection Method: Gamma, Poisson, Binomial Deviates 295 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 machinereadable 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). em=sq*y+xm; em is y, shifted and scaled. } while (em t); } return em; } Binomial Deviates If an event occurs with probability q, and we make n trials, then the number of times m that it occurs has the binomial distribution, j+ j− pn,q(m)dm = n j qj (1 − q) n−j (7.3.7) The binomial distribution is integer valued, with m taking on possible values from 0 to n. It depends on two parameters, n and q, so is correspondingly a bit harder to implement than our previous examples. Nevertheless, the techniques already illustrated are sufficiently powerful to do the job: #include #define PI 3.141592654 float bnldev(float pp, int n, long *idum) Returns as a floating-point number an integer value that is a random deviate drawn from a binomial distribution of n trials each of probability pp, using ran1(idum) as a source of uniform random deviates. { float gammln(float xx); float ran1(long *idum); int j; static int nold=(-1); float am,em,g,angle,p,bnl,sq,t,y; static float pold=(-1.0),pc,plog,pclog,en,oldg; p=(pp <= 0.5 ? pp : 1.0-pp); The binomial distribution is invariant under changing pp to 1-pp, if we also change the answer to n minus itself; we’ll remember to do this below. am=n*p; This is the mean of the deviate to be produced. if (n < 25) { Use the direct method while n is not too large. bnl=0.0; This can require up to 25 calls to ran1. for (j=1;j<=n;j++) if (ran1(idum) < p) ++bnl; } else if (am < 1.0) { If fewer than one event is expected out of 25 or more trials, then the distribution is quite accurately Poisson. Use direct Poisson method. g=exp(-am); t=1.0; for (j=0;j<=n;j++) { t *= ran1(idum); if (t < g) break; } bnl=(j <= n ? j : n); } else { Use the rejection method.
296 Chapter 7.Random Numbers if (n !nold){ If n has changed,then compute useful quanti- en=n; ties. oldg=gammln(en+1.0); nold=n; if (p !pold){ If p has changed,then compute useful quanti- pc=1.0-p; ties. plog=log(p); pclog=log(pc); pold=p; sq=sqrt(2.0*am*pc); The following code should by now seem familiar: do rejection method with a Lorentzian compar- do ison function. angle=PI*ran1(idum); y=tan(angle); em=sq*y+am; while (em 0.0 II em >(en+1.0)); Reject. em=floor(em); Trick for integer-valued distribution. t=1.2*sq*(1.0+y*y)*exp(oldg-gammln(em+1.0) -gammln(en-em+1.0)+em*plog+(en-em)*pclog); while (ran1(idum)>t); Reject.This happens about 1.5 times per devi- ate,on average if (p !pp)bnl=n-bnl; Remember to undo the symmetry transforma- A 令 return bnl; tion. Press. See Devroye [2]and Bratley [3]for many additional algorithms. CITED REFERENCES AND FURTHER READING: Knuth,D.E.1981.Seminumerical Algorithms,2nded.,vol.2 of The Art of Computer Programming IENTIFIC (Reading,MA:Addison-Wesley),pp.120ff.[1] 61 Devroye,L.1986,Non-Uniform Random Variate Generation (New York:Springer-Verlag),SX.4. 2 Bratley,P.,Fox,B.L.,and Schrage,E.L.1983,A Guide to Simulation (New York:Springer- Verlag）.[3l. COMPUTING (ISBN 188810920 Numerical Recipes 10.621 43108 7.4 Generation of Random Bits (outside The C language gives you useful access to some machine-level bitwise operations North Software. such as<<(left shift).This section will show you how to put such abilities to good use. The problem is how to generate single random bits,with 0 and 1 equally probable.Of course you can just generate uniform random deviates between zero and one and use their high-order bit(i.e.,test if they are greater than or less than 0.5).However this takes a lot of arithmetic;there are special-purpose applications, such as real-time signal processing,where you want to generate bits very much faster than that. One method for generating random bits,with two variant implementations,is based on "primitive polynomials modulo 2."The theory of these polynomials is beyond our scope (although 87.7 and $20.3 will give you small tastes of it).Here
296 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 machinereadable 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). if (n != nold) { If n has changed, then compute useful quantien=n; ties. oldg=gammln(en+1.0); nold=n; } if (p != pold) { If p has changed, then compute useful quantipc=1.0-p; ties. plog=log(p); pclog=log(pc); pold=p; } sq=sqrt(2.0*am*pc); The following code should by now seem familiar: rejection method with a Lorentzian comparison function. do { do { angle=PI*ran1(idum); y=tan(angle); em=sq*y+am; } while (em = (en+1.0)); Reject. em=floor(em); Trick for integer-valued distribution. t=1.2*sq*(1.0+y*y)*exp(oldg-gammln(em+1.0) -gammln(en-em+1.0)+em*plog+(en-em)*pclog); } while (ran1(idum) > t); Reject. This happens about 1.5 times per devibnl=em; ate, on average. } if (p != pp) bnl=n-bnl; Remember to undo the symmetry transformareturn bnl; tion. } See Devroye [2] and Bratley [3] for many additional algorithms. CITED REFERENCES AND FURTHER READING: Knuth, D.E. 1981, Seminumerical Algorithms, 2nd ed., vol. 2 of The Art of Computer Programming (Reading, MA: Addison-Wesley), pp. 120ff. [1] Devroye, L. 1986, Non-Uniform Random Variate Generation (New York: Springer-Verlag), §X.4. [2] Bratley, P., Fox, B.L., and Schrage, E.L. 1983, A Guide to Simulation (New York: SpringerVerlag). [3]. 7.4 Generation of Random Bits The C language gives you useful access to some machine-level bitwise operations such as << (left shift). This section will show you how to put such abilities to good use. The problem is how to generate single random bits, with 0 and 1 equally probable. Of course you can just generate uniform random deviates between zero and one and use their high-order bit (i.e., test if they are greater than or less than 0.5). However this takes a lot of arithmetic; there are special-purpose applications, such as real-time signal processing, where you want to generate bits very much faster than that. One method for generating random bits, with two variant implementations, is based on “primitive polynomials modulo 2.” The theory of these polynomials is beyond our scope (although §7.7 and §20.3 will give you small tastes of it). Here
