f7-2

7.2TransformationMethod: ExponentialandNormalDeviates 277Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Finally, the quick and dirty in-line generators ranqd1andranqd2are very fast, but they are machine dependent, nonportable, and at best only as good as a32-bitlinearcongruentialgeneratoreveris — in ourviewnot goodenoughinmany situations. We woulduse these onlyin veryspecial cases, wherespeed is critical. CITED REFERENCES AND FURTHER READING: Park, S.K., and Miller, K.W. 1988, Communications of the ACM , vol. 31, pp. 1192–1201. [1] Schrage, L. 1979, ACM Transactions on Mathematical Software , vol. 5, pp. 132–138. [2] Bratley, P., Fox, B.L., and Schrage, E.L. 1983, A Guide to Simulation (New York: Springer- Verlag). [3] Knuth,D.E.1981, SeminumericalAlgorithms ,2nded.,vol.2of TheArtofComputerProgramming (Reading, MA: Addison-Wesley), §§3.2–3.3. [4] Kahaner,D.,Moler,C.,andNash,S.1989, NumericalMethods andSoftware (EnglewoodCliffs, NJ: Prentice Hall), Chapter 10. [5] L’Ecuyer, P. 1988, Communications of the ACM , vol. 31, pp. 742–774. [6] Forsythe, G.E., Malcolm, M.A., and Moler, C.B. 1977, Computer Methods for Mathematical Computations (Englewood Cliffs, NJ: Prentice-Hall), Chapter 10. 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 xandx+dx, denoted p(x)dx, is given by p(x)dx=/braceleftBigdx 0<x< 1 0otherwise(7.2.1 ) The probability distribution p(x)is of course normalized, so that /integraldisplay∞ −∞p(x)dx=1 ( 7.2.2 ) Nowsupposethatwegenerateauniformdeviate xandthentakesomeprescribed functionofit, y(x). Theprobabilitydistributionof y,denoted p(y)dy,isdetermined by the fundamental transformation law of probabilities, which is simply |p(y)dy|=|p(x)dx| (7.2.3 ) or p(y)=p(x)/vextendsingle/vextendsingle/vextendsingle/vextendsingledx dy/vextendsingle/vextendsingle/vextendsingle/vextendsingle(7.2.4 ) 278 Chapter7. RandomNumbersSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).uniform deviate in 01 yxF(y) = 0p(y)dyy p(y)⌠ ⌡ transformed deviate out Figure 7.2.1. Transformation method for generating a random deviate yfrom a known probability distribution p(y). The inde finite integral of p(y)must be known and invertible. A uniform deviate xis chosen between 0and 1. Its corresponding yon the definite-integral curve is the desired deviate. Exponential Deviates As an example, suppose that y(x)≡− ln(x), and that p(x)is as given by equation (7.2.1) for a uniform deviate. Then p(y)dy=/vextendsingle/vextendsingle/vextendsingle/vextendsingledx dy/vextendsingle/vextendsingle/vextendsingle/vextendsingledy=e −ydy (7.2.5 ) which is distributed exponentially. This exponential distribution occurs frequently in real problems, usually as the distribution of waiting times between independent Poisson-random events, for example the radioactive decay of nuclei. You can also easily see (from7.2.4)that thequantity y/λhas theprobabilitydistribution λe−λy. So we have FUNCTION expdev(idum) INTEGER idumREAL expdev C USES ran1 Returns an exponentially distributed, positive, random deviate of unit mean, using ran1(idum) as the source of uniform deviates. REAL dum,ran1 1 dum=ran1(idum) if(dum.eq.0.)goto 1expdev=-log(dum)return END Let’sseewhatisinvolvedinusingtheabove transformationmethod togenerate somearbitrarydesireddistributionof y’s,sayonewith p(y)=f(y)forsomepositive function fwhose 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 ) 7.2TransformationMethod: ExponentialandNormalDeviates 279Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).But the solution of this is just x=F(y), where F(y)is the inde finite integral of f(y). Thedesiredtransformationwhichtakesauniformdeviateintoonedistributed asf(y)is therefore y(x)=F−1(x)( 7.2.7 ) where F−1is 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,eitheranalyticallyornumerically. Sometimesit is, andsometimesit 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, thenfind the value ythat has that fraction xof probability area to its left, and return the value y. Normal(Gaussian) Deviates Transformation methods generalize to more than one dimension. If x1,x2, ...are random deviates with a jointprobability distribution p(x1,x2,... ) dx 1dx 2..., and if y1,y2,...are each functions of all the x’s (same number of y’sa sx’s), then the joint probability distribution of the y’si s p(y1,y2,... )dy1dy2... =p(x1,x2,... )/vextendsingle/vextendsingle/vextendsingle/vextendsingle∂(x 1,x2,... ) ∂(y1,y2,... )/vextendsingle/vextendsingle/vextendsingle/vextendsingledy 1dy2... (7.2.8 ) where |∂()/∂()|is the Jacobian determinant of the x’s with respect to the y’s (or reciprocalof the Jacobian determinant of the y’s with respect to the x’s). 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), x1,x2, and two quantities y1,y2, y1=/radicalbig −2l nx1cos 2πx 2 y2=/radicalbig −2l nx1sin 2πx 2(7.2.10 ) Equivalently we can write x1=e x p/bracketleftbigg −1 2(y2 1+y2 2)/bracketrightbigg x2=1 2πarctany2 y1(7.2.11 ) Now the Jacobian determinant can readily be calculated (try it!): ∂(x1,x2) ∂(y1,y2)=/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle∂x 1 ∂y1∂x 1 ∂y2 ∂x 2 ∂y1∂x 2 ∂y2/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle=−/bracketleftbigg1 √ 2πe−y2 1/2/bracketrightbigg/bracketleftbigg1√ 2πe−y2 2/2/bracketrightbigg (7.2.12 ) 280 Chapter7. RandomNumbersSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Since this is the productofa functionof y2aloneanda functionof y1alone,we see thateach yis independentlydistributedaccordingtothe normaldistribution(7.2.9). Onefurthertrickisusefulinapplying(7.2.10). Supposethat,insteadofpicking uniform deviates x1andx2in the unit square, we instead pick v1andv2as the ordinateandabscissaofarandompointinsidetheunitcirclearoundtheorigin. Thenthesumoftheirsquares, R 2≡v2 1+v2 2isauniformdeviate,whichcanbeusedfor x1, whiletheanglethat (v1,v2)defineswithrespecttothe v1axiscanserveastherandom angle 2πx 2. What’s the advantage? It ’s that the cosine andsine in (7.2.10)can now be written as v1/√ R2andv2/√ R2, obviatingthe trigonometricfunctioncalls! We thus have FUNCTION gasdev(idum) INTEGER idum REAL gasdev C USES ran1 Returns a normally distributed deviate with zero mean and unit variance, using ran1(idum) as the source of uniform deviates. INTEGER isetREAL fac,gset,rsq,v1,v2,ran1SAVE iset,gset DATA iset/0/ if (idum.lt.0) iset=0 Reinitialize. if (iset.eq.0) then We don’t have an extra deviate handy, so 1 v1=2.*ran1(idum)-1. pick two uniform numbers in the square extend- i n gf r o m- 1t o+ 1i ne a c hd i r e c t i o n , v2=2.*ran1(idum)-1. rsq=v1**2+v2**2 see if they are in the unit circle, if(rsq.ge.1..or.rsq.eq.0.)goto 1 and if they are not, try again. fac=sqrt(-2.*log(rsq)/rsq) Now make the Box-Muller transformation to get two normal deviates. Return one and savethe other for next time.gset=v1*fac gasdev=v2*fac iset=1 Set flag. else We have an extra deviate handy, gasdev=gset so return it, iset=0 and unset the flag. endif returnEND 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, SeminumericalAlgorithms ,2nded.,vol.2of TheArtofComputerProgramming (Reading, MA: Addison-Wesley), pp. 116ff. 7.3RejectionMethod: Gamma, Poisson,BinomialDeviates 281Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).7.3 Rejection Method: Gamma, Poisson, Binomial Deviates Therejection method is a powerful, general technique for generating random deviateswhosedistributionfunction p(x)dx(probabilityofavalueoccurringbetween xandx+dx) is known and computable. The rejection method does notrequire that the cumulative distribution function [inde finite 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 probabilitydistribution p(x)that you wish to generate, so thattheareaunderthecurveinanyrangeof xcorrespondstothedesiredprobability ofgeneratingan xin that range. If we hadsomeway ofchoosinga randompoint in two dimensions , with uniform probability in the areaunder 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)areaandlies everywhere aboveyouroriginalprobabilitydistribution. (This isalwayspossible,becauseyouroriginalcurveenclosesonlyunitarea,byde finition ofprobability.) Wewillcallthis f(x)thecomparisonfunction . Imaginenowthatyou havesomewayofchoosinga randompointintwo dimensionsthatis uniforminthe areaunderthecomparisonfunction. Wheneverthatpointlies outsidetheareaunderthe original probability distribution, we will rejectit and choose another random point. Whenever it lies inside the area under the original probability distribution, we willacceptit. It should be obvious that the accepted points are uniform in the accepted area, so that their xvalues have the desired distribution. It 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 theprobability 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 Ais the total area under f(x), and use it to get a corresponding x. Then pick a uniformdeviate between 0 and f(x)as the y valueforthetwo-dimensionalpoint. Youshouldbeabletoconvinceyourselfthatthe point (x, y )is uniformlydistributedintheareaunderthecomparisonfunction 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 andfor all, some reasonablygoodcomparisonfunction f(x). Thereafter, eachdeviategeneratedrequirestwouniformrandomdeviates,oneevaluationof f(to get the coordinate y), and one evaluationof p(to decide whether to accept or reject