7.2 Transformation Method:Exponential and Normal Deviates 287 7.2 Transformation Method:Exponential and Normal Deviates In the previous section,we learned how to generate random deviates with a uniform probability distribution,so that the probability of generating a number between x and x+dr,denoted p(x)dr,is given by p(x)dx ∫dx0 (outside North Software. float expdev(long *idum) Returns an exponentially distributed,positive,random deviate of unit mean, using ran1(idum)as the source of uniform deviates. Ame float ran1(long *idum); float dum; do dum=ran1(idun）; while (dum ==0.0); return -log(dum);
7.2 Transformation Method: Exponential and Normal Deviates 287 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). 7.2 Transformation Method: Exponential and Normal Deviates In the previous section, we learned how to generate random deviates with a uniform probability distribution, so that the probability of generating a number between x and x + dx, denoted p(x)dx, is given by p(x)dx = dx 0 float expdev(long *idum) Returns an exponentially distributed, positive, random deviate of unit mean, using ran1(idum) as the source of uniform deviates. { float ran1(long *idum); float dum; do dum=ran1(idum); while (dum == 0.0); return -log(dum); }
288 Chapter 7.Random Numbers uniform deviate in F(y)=Jop(y)dy -p(y) 0 transformed 83g deviate out Figure 7.2.1.Transformation method for generating a random deviate y from a known probability 1.800 distribution p(y).The indefinite integral of p(y)must be known and invertible.A uniform deviate x is chosen between 0 and 1.Its corresponding y on the definite-integral curve is the desired deviate. Let's see what is involved in using the above transformation method to generate some arbitrary desired distribution ofy's,say one with p(y)=f(y)for some positive function f whose integral is 1.(See Figure 7.2.1.)According to (7.2.4),we need to solve the differential equation Press. dx dy =f(u) (7.2.6) But the solution of this is just x =F(y),where F(y)is the indefinite integral of f(y).The desired transformation which takes a uniform deviate into one distributed IENTIFIC as f(y)is therefore 6 y(x)=F-1(x) (7.2.7) where F-1 is the inverse function to F.Whether(7.2.7)is feasible to implement 景影 depends on whether the inverse function of the integral of fry)is itself feasible to compute,either analytically or numerically.Sometimes it is,and sometimes it isn't. Incidentally,(7.2.7)has an immediate geometric interpretation:Since F(y)is Numerica 10.621 the area under the probability curve to the left of y,(7.2.7)is just the prescription: 431 choose a uniform random then find the value y that has that fraction x of E probability area to its left,and return the value y. (outside Recipes Normal(Gaussian)Deviates North Software. Transformation methods generalize to more than one dimension.If z1,z2, are random deviates with a joint probability distribution p(1,2,...) dzidr2...,and if y,y2,...are each functions of all the x's (same number of y's as z's),then the joint probability distribution of the y's is p(1,2,)d1d2.=p(x1,2,.) ∂(x1,x2,.】 dyidy2... (7.2.8) 01,2,） where o()/o(is the Jacobian determinant of the x's with respect to the y's (or reciprocal of the Jacobian determinant of the y's with respect to the z's)
288 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). uniform deviate in 0 1 y x F(y) = 0 p(y)dy y p(y) ⌠ ⌡ transformed deviate out Figure 7.2.1. Transformation method for generating a random deviate y from a known probability distribution p(y). The indefinite integral of p(y) must be known and invertible. A uniform deviate x is chosen between 0 and 1. Its corresponding y on the definite-integral curve is the desired deviate. Let’s see what is involved in using the above transformation method to generate some arbitrary desired distribution of y’s, say one with p(y) = f(y)for some positive function f whose integral is 1. (See Figure 7.2.1.) According to (7.2.4), we need to solve the differential equation dx dy = f(y) (7.2.6) But the solution of this is just x = F(y), where F(y) is the indefinite integral of f(y). The desired transformation which takes a uniform deviate into one distributed as f(y) is therefore y(x) = F −1(x) (7.2.7) where F −1 is the inverse function to F. Whether (7.2.7) is feasible to implement depends on whether the inverse function of the integral of f(y) is itself feasible to compute, either analytically or numerically. Sometimes it is, and sometimes it isn’t. Incidentally, (7.2.7) has an immediate geometric interpretation: Since F(y) is the area under the probability curve to the left of y, (7.2.7) is just the prescription: choose a uniform random x, then find the value y that has that fraction x of probability area to its left, and return the value y. Normal (Gaussian) Deviates Transformation methods generalize to more than one dimension. If x 1, x2, ... are random deviates with a joint probability distribution p(x1, x2,...) dx1dx2 ... , and if y1, y2,... are each functions of all the x’s (same number of y’s as x’s), then the joint probability distribution of the y’s is p(y1, y2,...)dy1dy2 ... = p(x1, x2,...) ∂(x1, x2,...) ∂(y1, y2,...) dy1dy2 ... (7.2.8) where |∂( )/∂( )| is the Jacobian determinant of the x’s with respect to the y’s (or reciprocal of the Jacobian determinant of the y’s with respect to the x’s).
7.2 Transformation Method:Exponential and Normal Deviates 289 An important example of the use of (7.2.8)is the Box-Muller method for generating random deviates with a normal(Gaussian)distribution, 1 p(y)dy _e-v'/2dy V2 (7.2.9) Consider the transformation between two uniform deviates on(0,1),z1,z2,and two quantities y1,2, y1 =V-21n 1 cos 2Tx2 (7.2.10) 2=√-2lnx1sin2xx2 Equivalently we can write 3 (7.2.11) RECIPES x2= 2 arctan 1 Now the Jacobian determinant can readily be calculated (try it!) 气 Press. 1 1 0(x1,x2】 /2 OT. (7.2.12) 0(y1,y2 Programs Since this is the product of a function of y2 alone and a function of y alone,we see IENTIFIC that each y is independently distributed according to the normal distribution(7.2.9). 6 One further trick is useful in applying(7.2.10).Suppose that,instead of picking uniform deviates and 2 in the unit square,we instead pick v and v2 as the ordinate and abscissa of a random point inside the unit circle around the origin.Then the sum of their squares,R2=v is a uniform deviate,which can be used for1, 1920 COMPUTING (ISBN while the angle that(v1,v2)defines with respect to the v axis can serve as the random angle 22.What's the advantage?It's that the cosine and sine in(7.2.10)can now Numerical Recipes 10-621 be written as v/vR2 and v2/vR2,obviating the trigonometric function calls! We thus have 43106 #include (outside float gasdev(long *idum) Software. Returns a normally distributed deviate with zero mean and unit variance,using ran1(idum) North as the source of uniform deviates. float ran1(long *idum); static int iset=0; static float gset; float fac,rsq,v1,v2; 1f(*1dum<0)1set=0: Reinitialize. if (iset ==0){ We don't have an extra deviate handy,so do v1=2.0*ran1(1dum)-1.0; pick two uniform numbers in the square ex- v2=2.0*ran1(1dum)-1.0; tending from-1 to +1 in each direction, rsq=v1*v1+v2*v2; see if they are in the unit circle
7.2 Transformation Method: Exponential and Normal Deviates 289 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). An important example of the use of (7.2.8) is the Box-Muller method for generating random deviates with a normal (Gaussian) distribution, p(y)dy = 1 √ 2π e−y2/2dy (7.2.9) Consider the transformation between two uniform deviates on (0,1), x 1, x2, and two quantities y1, y2, y1 = −2 ln x1 cos 2πx2 y2 = −2 ln x1 sin 2πx2 (7.2.10) Equivalently we can write x1 = exp −1 2 (y2 1 + y2 2) x2 = 1 2π arctany2 y1 (7.2.11) Now the Jacobian determinant can readily be calculated (try it!): ∂(x1, x2) ∂(y1, y2) = ∂x1 ∂y1 ∂x1 ∂y2 ∂x2 ∂y1 ∂x2 ∂y2 = − 1 √ 2π e−y2 1/2 1 √ 2π e−y2 2/2 (7.2.12) Since this is the product of a function of y2 alone and a function of y1 alone, we see that each y is independently distributed according to the normal distribution (7.2.9). One further trick is useful in applying (7.2.10). Suppose that, instead of picking uniform deviates x1 and x2 in the unit square, we instead pick v1 and v2 as the ordinate and abscissa of a random point inside the unit circle around the origin. Then the sum of their squares, R2 ≡ v2 1 +v2 2 is a uniform deviate, which can be used for x1, while the angle that(v1, v2) defines with respect to the v1 axis can serve as the random angle 2πx2. What’s the advantage? It’s that the cosine and sine in (7.2.10) can now be written as v1/ √ R2 and v2/ √ R2, obviating the trigonometric function calls! We thus have #include float gasdev(long *idum) Returns a normally distributed deviate with zero mean and unit variance, using ran1(idum) as the source of uniform deviates. { float ran1(long *idum); static int iset=0; static float gset; float fac,rsq,v1,v2; if (*idum < 0) iset=0; Reinitialize. if (iset == 0) { We don’t have an extra deviate handy, so do { v1=2.0*ran1(idum)-1.0; pick two uniform numbers in the square exv2=2.0*ran1(idum)-1.0; tending from -1 to +1 in each direction, rsq=v1*v1+v2*v2; see if they are in the unit circle,
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
