c7-3 (779517), страница 2
Текст из файла (страница 2)
Sincethe sum of logarithms is the logarithm of the product, one really has only to generatethe product of a uniform 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 = a and a half-width of about a.We will be interested in several probability distributions with this same qualitative form. A useful comparison function in such cases is derived from theLorentzian distribution11p(y)dy =dy(7.3.2)π 1 + y27.3 Rejection Method: Gamma, Poisson, Binomial Deviates}return x;These four lines generate the tangent of a random angle, i.e., theyare equivalent toy = tan(π * ran1(idum)).We decide whether to reject x:Reject in region of zero probability.Ratio of prob.
fn. to comparison fn.Reject on basis of a second uniformdeviate.}Poisson DeviatesThe Poisson distribution is conceptually related to the gamma distribution. Itgives the probability of a certain integer number m of unit rate Poisson randomevents occurring in a given interval of time x, while the gamma distribution was theprobability of waiting time between x and x + dx to the mth event.
Note that m takeson only integer values ≥ 0, so that the Poisson distribution, viewed as a continuousdistribution 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 regioncontaining the integer is some finite number. The total probability at an integer j isZj+Prob(j) =px (m)dm =j−xj e−xj!(7.3.5)At first sight this might seem an unlikely candidate distribution for the rejectionmethod, since no continuous comparison function can be larger than the infinitelytall, but infinitely narrow, Dirac delta functions in px (m). However, there is a trickthat we can do: Spread the finite area in the spike at j uniformly into the intervalbetween j and j + 1.
This defines a continuous distribution qx (m)dm given byqx (m)dm =x[m] e−xdm[m]!(7.3.6)where [m] represents the largest integer less than m. If we now use the rejectionmethod to generate a (noninteger) deviate from (7.3.6), and then take the integerpart of that deviate, it will be as if drawn from the desired distribution (7.3.5). (SeeFigure 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 (albeitwith 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 cangenerate independent exponential deviates (waiting times between events); when thesum of these first exceeds x, then the number of events that would have occurred inwaiting time x becomes known and is one less than the number of terms in the sum.These ideas produce the following routine:Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).do {do {do {v1=ran1(idum);v2=2.0*ran1(idum)-1.0;} while (v1*v1+v2*v2 > 1.0);y=v2/v1;am=ia-1;s=sqrt(2.0*am+1.0);x=s*y+am;} while (x <= 0.0);e=(1.0+y*y)*exp(am*log(x/am)-s*y);} while (ran1(idum) > e);293294Chapter 7.Random Numbers1inaccept012345Figure 7.3.2. Rejection method as applied to an integer-valued distribution.
The method is performedon the step function shown as a dashed line, yielding a real-valued deviate. This deviate is roundeddown to the next lower integer, which is output.#include <math.h>#define PI 3.141592654float poidev(float xm, long *idum)Returns as a floating-point number an integer value that is a random deviate drawn from aPoisson distribution of mean xm, using ran1(idum) as a source of uniform random deviates.{float gammln(float xx);float ran1(long *idum);static float sq,alxm,g,oldm=(-1.0);oldm is a flag for whether xm has changedfloat em,t,y;since last call.if (xm < 12.0) {Use direct method.if (xm != oldm) {oldm=xm;g=exp(-xm);If xm is new, compute the exponential.}em = -1;t=1.0;do {Instead of adding exponential deviates it is equiv++em;alent to multiply uniform deviates.
We nevert *= ran1(idum);actually have to take the log, merely com} while (t > g);pare to the pre-computed exponential.} else {Use rejection method.if (xm != oldm) {If xm has changed since the last call, then preoldm=xm;compute some functions that occur below.sq=sqrt(2.0*xm);alxm=log(xm);g=xm*alxm-gammln(xm+1.0);The function gammln is the natural log of the gamma function, as given in §6.1.}do {do {y is a deviate from a Lorentzian comparison funcy=tan(PI*ran1(idum));tion.Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).reject7.3 Rejection Method: Gamma, Poisson, Binomial Deviates295}return em;}Binomial DeviatesIf an event occurs with probability q, and we make n trials, then the number oftimes m that it occurs has the binomial distribution,Zj+j− n jpn,q (m)dm =q (1 − q)n−jj(7.3.7)The binomial distribution is integer valued, with m taking on possible valuesfrom 0 to n.
It depends on two parameters, n and q, so is correspondingly abit harder to implement than our previous examples. Nevertheless, the techniquesalready illustrated are sufficiently powerful to do the job:#include <math.h>#define PI 3.141592654float bnldev(float pp, int n, long *idum)Returns as a floating-point number an integer value that is a random deviate drawn froma binomial distribution of n trials each of probability pp, using ran1(idum) as a source ofuniform 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 theanswer 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 25g=exp(-am);or more trials, then the distribution is quitet=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.Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (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 orreject by comparing it to another uniform deviate.
The factor 0.9 is chosen sothat t never exceeds 1.} while (ran1(idum) > t);296Chapter 7.Random Numbers}if (p != pp) bnl=n-bnl;return bnl;Remember to undo the symmetry transformation.}See Devroye [2] and Bratley [3] for many additional algorithms.CITED REFERENCES AND FURTHER READING:Knuth, D.E. 1981, Seminumerical Algorithms, 2nd ed., vol. 2 of The Art of Computer Programming(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: SpringerVerlag). [3].7.4 Generation of Random BitsThe C language gives you useful access to some machine-level bitwise operationssuch as << (left shift). This section will show you how to put such abilities to good use.The problem is how to generate single random bits, with 0 and 1 equallyprobable. 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 than0.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, isbased on “primitive polynomials modulo 2.” The theory of these polynomials isbeyond our scope (although §7.7 and §20.3 will give you small tastes of it).
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
