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(x<=0.0); Reject in region of zero probability. e=(1.0+y*y)*exp(am*log(x/am)-s*y); Ratio of prob.fn.to comparison fn. while (ran1(idum)>e); 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 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). do { do { do { These four lines generate the tan￾gent 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 <= 0.0); Reject in region of zero probability. e=(1.0+y*y)*exp(am*log(x/am)-s*y); Ratio of prob. fn. to comparison fn. } while (ran1(idum) > 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: