f7-3

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 [indefinite integral of p(x)] be readily computable, much less the inverse of that function — which was required for thetransformation 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,bydefinition 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 282 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).second random deviate in 0f(x0)reject x0 accept x0first random deviate inA 0f(x)0f(x)dxx p(x) x0⌠ ⌡ Figure 7.3.1. Rejection method for generating a random deviate xfrom aknown 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 xof 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 fis found; and so on. The ratio of accepted to rejected points is the ratio of the area under pto the area between pand f. thepoint x, y). Figure7.3.1illustratestheprocedure. Then,ofcourse,thisprocedure must be repeated, on the average, Atimes beforethe final deviate is obtained. GammaDistribution The gamma distribution of integer order a> 0is the waiting time to the ath eventina Poissonrandomprocessofunitmean. Forexample,when a=1,it is just the exponential distribution of §7.2, the waiting time to the first event. A gamma deviate has probability pa(x)dxof occurring with a value between xandx+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. Sincethesumoflogarithmsisthelogarithmoftheproduct,onereallyhasonlytogenerate the product of auniform 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=aand 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 π/parenleftbigg1 1+y2/parenrightbigg dy (7.3.2 ) whose inverse inde finite integral is just the tangent function. It follows that the x-coordinateof an area-uniformrandompoint under the comparisonfunction f(x)=c0 1+( x−x0)2/a2 0(7.3.3 ) 7.3RejectionMethod: Gamma, Poisson,BinomialDeviates 283Sample 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).for any constants a0,c0,andx0, can be generated by the prescription x=a0tan(πU)+x0 (7.3.4 ) where Uisauniformdeviatebetween0and1. Thus,forsomespeci fic“bell-shaped ” p(x)probabilitydistribution,weneedonly findconstants a0,c0,x0,withtheproduct a0c0(whichdeterminesthearea)assmallaspossible,suchthat(7.3.3)iseverywhere greater than p(x). Ahrens has done this for the gamma distribution, yielding the following algorithm (as described in Knuth [1]): FUNCTION gamdev(ia,idum) INTEGER ia,idum REAL gamdev C USES ran1 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. INTEGER jREAL am,e,s,v1,v2,x,y,ran1if(ia.lt.1)pause ’bad argument in gamdev’ if(ia.lt.6)then Use direct method, adding waiting times. x=1.do 11j=1,ia x=x*ran1(idum) enddo 11 x=-log(x) else Use rejection method. 1 v1=ran1(idum) These four lines generate the tangent of a random angle, i.e., are equivalent to y = tan(3.14159265 * ran1(idum)) . v2=2.*ran1(idum)-1. if(v1**2+v2**2.gt.1.)goto 1 y=v2/v1 am=ia-1s=sqrt(2.*am+1.)x=s*y+am We decide whether to reject x: if(x.le.0.)goto 1 Reject in region of zero probability. e=(1.+y**2)*exp(am*log(x/am)-s*y) Ratio of prob. fn. to comparison fn. if(ran1(idum).gt.e)goto 1 Reject on basis of a second uniform de- viate. endif gamdev=xreturnEND Poisson Deviates The Poisson distribution is conceptually related to the gamma distribution. It gives the probability of a certain integer number mof unit rate Poisson random eventsoccurringin a givenintervaloftime x, while the gammadistributionwas the probabilityofwaitingtimebetween xandx+dxtothe mthevent. Notethat mtakes on only integer values ≥0, so that the Poisson distribution, viewed as a continuous distribution function px(m)dm, is zero everywhere except where mis an integer ≥0. At such places, it is in finite, such that the integrated probability over a region containingtheintegeris some finite number. Thetotal probabilityatan integer jis Prob (j)=/integraldisplayj+/epsilon1 j−/epsilon1px(m)dm=xje−x j!(7.3.5 ) 284 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).012345acceptrejectin1 Figure 7.3.2. Rejection method as applied to an integer-valued distribution. The method is performed on the step function shown asa dashed line, yielding a real-valued deviate. This deviate isrounded downto the next lower integer, which is output. Atfirstsightthismightseemanunlikelycandidatedistributionfortherejection method, since no continuous comparison function can be larger than the in finitely tall, but in finitely narrow, Dirac delta functions inpx(m). However, there is a trick that we can do: Spread the finite area in the spike at juniformly into the interval between jandj+1. This defines a continuousdistribution qx(m)dmgivenby 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 Figure7.3.2.) This trickis generalforanyinteger-valuedprobabilitydistribution. Forxlarge 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 ofLorentzian comparison function as was already used above. For small x, we can generateindependentexponentialdeviates(waitingtimesbetweenevents);whenthe sum of these first exceeds x, then the numberof events that would have occurredin waitingtime xbecomesknownandis oneless thanthenumberoftermsinthe sum. These ideas produce the following routine: FUNCTION poidev(xm,idum) INTEGER idumREAL poidev,xm,PI PARAMETER (PI=3.141592654) C USES gammln,ran1 Returns as a floating-point number an integer value that is a random deviate drawn from aPoisson distribution of mean xm,u s i n g ran1(idum) as a source of uniform random deviates. 7.3RejectionMethod: Gamma, Poisson,BinomialDeviates 285Sample 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).REAL alxm,em,g,oldm,sq,t,y,gammln,ran1 SAVE alxm,g,oldm,sq DATA oldm /-1./ Flag for whether xmhas changed since last call. if (xm.lt.12.)then Use direct method. if (xm.ne.oldm) then oldm=xm g=exp(-xm) Ifxmis new, compute the exponential. endifem=-1 t=1. 2 em=em+1. Instead of adding exponential deviates it is equivalent to mul- tiply uniform deviates. We never actually have to take thelog, merely compare to the pre-computed exponential.t=t*ran1(idum) if (t.gt.g) goto 2 else Use rejection method. if (xm.ne.oldm) then Ifxmhas changed since the last call, then precompute some functions that occur below. oldm=xm sq=sqrt(2.*xm) alxm=log(xm)g=xm*alxm-gammln(xm+1.) The function gammln is the natural log of the gamma function, as given in §6.1. endif 1 y=tan(PI*ran1(idum)) y is a deviate from a Lorentzian comparison function. em=sq*y+xm em isy, shifted and scaled. if (em.lt.0.) goto 1 Reject if in regime of zero probability. em=int(em) The trick for integer-valued distributions. t=0.9*(1.+y**2)*exp(em*alxm-gammln(em+1.)-g) The ratio of the desired distribu- tion to the comparison function; we accept or re-ject by comparing it to another uniform deviate. The factor 0.9 is chosen so that tnever exceeds 1.if (ran1(idum).gt.t) goto 1 endif poidev=em returnEND BinomialDeviates If an event occurs with probability q, and we make ntrials, then the number of times mthat it occurs has the binomial distribution, /integraldisplayj+/epsilon1 j−/epsilon1pn,q(m)dm=/parenleftbiggn j/parenrightbigg qj(1−q)n−j(7.3.7 ) The binomial distribution is integer valued, with mtaking on possible values from 0 to n. It depends on twoparameters, nandq, so is correspondingly a bit harder to implement than our previous examples. Nevertheless, the techniques already illustrated are suf ficiently powerful to do the job: FUNCTION bnldev(pp,n,idum) INTEGER idum,n REAL bnldev,pp,PI C USES gammln,ran1 PARAMETER (PI=3.141592654) Returns as a floating-point number an integer value that is a random deviate drawn from a binomial distribution of ntrials each of probability pp,u s i n g ran1(idum) as a source of uniform random deviates. INTEGER j,nold REAL am,em,en,g,oldg,p,pc,pclog,plog,pold,sq,t,y,gammln,ran1 286 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).SAVE nold,pold,pc,plog,pclog,en,oldg DATA nold /-1/, pold /-1./ Arguments from previous calls. if(pp.le.0.5)then The binomial distribution is invariant under changing ppto 1.-pp , if we also change the answer to nminus itself; we’ll remember to do this below.p=pp else p=1.-pp endifam=n*p T h i si st h em e a no ft h ed e v i a t et ob ep r o d u c e d . if (n.lt.25)then Use the direct method while nis not too large. This can require up to 25 calls to ran1 . bnldev=0. do 11j=1,n if(ran1(idum).lt.p)bnldev=bnldev+1. enddo 11 else if (am.lt.1.) then If fewer than one event is expected out of 25 or more tri- als, then the distribution is quite accurately Poisson. Usedirect Poisson method.g=exp(-am) t=1. do 12j=0,n t=t*ran1(idum)if (t.lt.g) goto 1 enddo 12 j=n 1 bnldev=j else Use the rejection method. if (n.ne.nold) then Ifnhas changed, then compute useful quantities. en=noldg=gammln(en+1.)nold=n endif if (p.ne.pold) then Ifphas changed, then compute useful quantities. pc=1.-p plog=log(p) pclog=log(pc)pold=p endif sq=sqrt(2.*am*pc) The following code should by now seem familiar: rejection method with a Lorentzian comparison function. 2 y=tan(PI*ran1(idum)) em=sq*y+amif (em.lt.0..or.em.ge.en+1.) goto 2 Reject. em=int(em) Trick for integer-valued distribution. t=1.2*sq*(1.+y**2)*exp(oldg-gammln(em+1.) * -gammln(en-em+1.)+em*plog+(en-em)*pclog) if (ran1(idum).gt.t) goto 2 Reject. This happens about 1.5 times per deviate, on average. bnldev=em endifif (p.ne.pp) bnldev=n-bnldev Remember to undo the symmetry transformation. return END See Devroye [2]and Bratley [3]for many additional algorithms. CITED REFERENCES AND FURTHER READING: Knuth,D.E.1981, SeminumericalAlgorithms ,2nded.,vol.2of TheArtofComputerProgramming (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: Springer- Verlag). [3]. 7.4 GenerationofRandomBits 287Sample 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.4 Generation of Random Bits This topic is not very useful for programming in high-level languages, but it can be quite useful when you have access to the machine-language level of a machine or when you are in a position to build special-purpose hardware out of readily available chips. 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 zeroand 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 muchfaster 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, suffice it to say that there are special polynomials among those whose coef ficients are zero or one. An example is x 18+x5+x2+x1+x0(7.4.1 ) which we can abbreviate by just writing the nonzero powers of x, e.g., (18,5,2,1,0) Everyprimitivepolynomialmodulo2oforder n(=18above)de finesarecurrence relation for obtaining a new random bit from the npreceding ones. The recurrence relation is guaranteed to produce a sequence of maximal length, i.e., cycle through all possible sequences of nbits (except all zeros) before it repeats. Therefore one can seed the sequence with any initial bit pattern (except all zeros), and get 2n−1 random bits before the sequence repeats. Letthebitsbenumberedfrom1(mostrecentlygenerated)through n(generated nsteps ago), and denoted a1,a2,...,a n. We want to give a formula for a new bit a0. After generating a0we will shift all the bits by one, so that the old anisfinally lost, and the new a0becomes a1. We then applythe formulaagain, and so on. “MethodI”istheeasiesttoimplementinhardware,requiringonlyasingleshift register nbits longanda few XOR ( “exclusiveor ”or bit additionmod2)gates.For the primitive polynomial given above, the recurrence formula is a0=a18XOR a5XOR a2XOR a1 (7.4.2 ) The terms that are XOR ’d together can be thought of as “taps”on the shift register, XOR’d into the register ’s input. More generally, there is precisely one term for each nonzero coef ficient in the primitive polynomial except the constant (zero bit) term. So the first term will always be anfor a primitive polynomial of degree n, while the last term might or might not be a1, depending on whether the primitive polynomial has a term in x1. Itisrathercumbersometoillustratethemethodin FORTRAN. Assumethat iand is a bitwise AND function, notis bitwise complement, ishft( ,1) is leftshift by onebit, iorisbitwiseOR.(Theseareavailableinmany FORTRAN implementations.) Then we have the following routine.