7.3 Rejection Method:Gamma,Poisson,Binomial Deviates 295 em=sq*y+xm; em is y,shifted and scaled. ]whi1e(em<0.0); Reject if in regime of zero probability. em=floor(em); The trick for integer-valued distributions t=0.9*(1.0+y*y)*exp(em*alxm-gammln(em+1.0)-g); The ratio of the desired distribution to the comparison function;we accept or reject by comparing it to another uniform deviate.The factor 0.9 is chosen so that t never exceeds 1. while (ran1(idum)>t); return em; Binomial Deviates 83 鱼 granted for 18881892 If an event occurs with probability g,and we make n trials,then the number of 11800 times m that it occurs has the binomial distribution. from NUMERICAL RECIPESI ]+6 Pn.(m)dm (1-q)n-j (7.3.7) 芝是%6sadN (Nort server The binomial distribution is integer valued,with m taking on possible values America THE from 0 to n.It depends on two parameters,n and g,so is correspondingly a ART bit harder to implement than our previous examples.Nevertheless,the techniques already illustrated are sufficiently powerful to do the job: 9 Program #include <math.h> #def1nePI3.141592654 float bnldev(float pp,int n,long *idum) Returns as a floating-point number an integer value that is a random deviate drawn from a binomial distribution of n trials each of probability pp.using ran1(idum)as a source of uniform random deviates. OF SCIENTIFIC COMPUTING(ISBN 1920 float gammln(float xx); float ran1(long *idum); int j; Numerical Recipes 2-521 static int nold=(-1); float am,em,g,angle,p,bnl,sq,t,y; 43108 static float pold=(-1.0),pc,plog,pclog,en,oldg; p=(pp<=0.5？pp:1.0-pp); The binomial distribution is invariant under changing pp to 1-pp,if we also change the (outside answer to n minus itself;we'll remember to do this below. North Software. am=n*p; This is the mean of the deviate to be produced. 1f(n<25) Use the direct method while n is not too large 6 bnl=0.0; This can require up to 25 calls to ran1. for (j=1;j<=n;j++) if (ran1(idum)<p)++bnl; 】e1seif(am<1.0)[ If fewer than one event is expected out of 25 g=exp(-am) or more trials,then the distribution is quite t=1.0: accurately Poisson.Use direct Poisson method. for (j=0;j<=n;j++){ t *ran1(idum); if (t g)break; bnl=(j <n j:n); else Use the rejection method.7.3 Rejection Method: Gamma, Poisson, Binomial Deviates 295 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). em=sq*y+xm; em is y, shifted and scaled. } while (em < 0.0); Reject if in regime of zero probability. em=floor(em); The trick for integer-valued distributions. t=0.9*(1.0+y*y)*exp(em*alxm-gammln(em+1.0)-g); The ratio of the desired distribution to the comparison function; we accept or reject by comparing it to another uniform deviate. The factor 0.9 is chosen so that t never exceeds 1. } while (ran1(idum) > t); } return em; } Binomial Deviates If an event occurs with probability q, and we make n trials, then the number of times m that it occurs has the binomial distribution, j+ j− pn,q(m)dm =
n j qj (1 − q) n−j (7.3.7) The binomial distribution is integer valued, with m taking on possible values from 0 to n. It depends on two parameters, n and q, so is correspondingly a bit harder to implement than our previous examples. Nevertheless, the techniques already illustrated are sufficiently powerful to do the job: #include <math.h> #define PI 3.141592654 float bnldev(float pp, int n, long *idum) Returns as a floating-point number an integer value that is a random deviate drawn from a binomial distribution of n trials each of probability pp, using ran1(idum) as a source of uniform random deviates. { float gammln(float xx); float ran1(long *idum); int j; static int nold=(-1); float am,em,g,angle,p,bnl,sq,t,y; static float pold=(-1.0),pc,plog,pclog,en,oldg; p=(pp <= 0.5 ? pp : 1.0-pp); The binomial distribution is invariant under changing pp to 1-pp, if we also change the answer to n minus itself; we’ll remember to do this below. am=n*p; This is the mean of the deviate to be produced. if (n < 25) { Use the direct method while n is not too large. bnl=0.0; This can require up to 25 calls to ran1. for (j=1;j<=n;j++) if (ran1(idum) < p) ++bnl; } else if (am < 1.0) { If fewer than one event is expected out of 25 or more trials, then the distribution is quite accurately Poisson. Use direct Poisson method. g=exp(-am); t=1.0; for (j=0;j<=n;j++) { t *= ran1(idum); if (t < g) break; } bnl=(j <= n ? j : n); } else { Use the rejection method.