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.It290 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). } 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: 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. 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