c7-3 (779517)
Текст из файла
290Chapter 7.Random Numbers}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: SpringerVerlag). [2]Knuth, D.E. 1981, Seminumerical Algorithms, 2nd ed., vol. 2 of The Art of Computer Programming(Reading, MA: Addison-Wesley), pp. 116ff.7.3 Rejection Method: Gamma, Poisson,Binomial DeviatesThe rejection method is a powerful, general technique for generating randomdeviates whose distribution function p(x)dx (probability of a value occurring betweenx and x + dx) is known and computable. The rejection method does not requirethat the cumulative distribution function [indefinite integral of p(x)] be readilycomputable, 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 probability distribution p(x) that you wish to generate, sothat the area under the curve in any range of x corresponds to the desired probabilityof generating an x in that range.
If we had some way of choosing a random point intwo dimensions, with uniform probability in the area under your curve, then the xvalue of that random point would have the desired distribution.Now, on the same graph, draw any other curve f(x) which has finite (notinfinite) area and lies everywhere above your original probability distribution. (Thisis always possible, because your original curve encloses only unit area, by definitionof probability.) We will call this f(x) the comparison function. Imagine nowthat you have some way of choosing a random point in two dimensions that isuniform in the area under the comparison function.
Whenever that point lies outsidethe area under the original probability distribution, we will reject it and chooseanother random point. Whenever it lies inside the area under the original probabilitydistribution, we will accept it. It should be obvious that the accepted points areuniform in the accepted area, so that their x values have the desired distribution. ItSample 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).} while (rsq >= 1.0 || rsq == 0.0);and if they are not, try again.fac=sqrt(-2.0*log(rsq)/rsq);Now make the Box-Muller transformation to get two normal deviates. Return one andsave the other for next time.gset=v1*fac;iset=1;Set flag.return v2*fac;} else {We have an extra deviate handy,iset=0;so unset the flag,return gset;and return it.}7.3 Rejection Method: Gamma, Poisson, Binomial Deviates291Afirst randomdeviate in⌠x⌡0f(x)dxf(x0 )accept x0second randomdeviate inp(x)0x00Figure 7.3.1. Rejection method for generating a random deviate x from a known probability distributionp(x) that is everywhere less than some other function f (x).
The transformation method is first used togenerate a random deviate x of the distribution f (compare Figure 7.2.1). A second uniform deviate isused to decide whether to accept or reject that x. If it is rejected, a new deviate of f is found; and so on.The ratio of accepted to rejected points is the ratio of the area under p to the area between p and f .should also be obvious that the fraction of points rejected just depends on the ratioof the area of the comparison function to the area of the probability distributionfunction, not on the details of shape of either function.
For example, a comparisonfunction whose area is less than 2 will reject fewer than half the points, even if itapproximates the probability function very badly at some values of x, e.g., remainsfinite in some region where x is zero.It remains only to suggest how to choose a uniform random point in twodimensions under the comparison function f(x). A variant of the transformationmethod (§7.2) does nicely: Be sure to have chosen a comparison function whoseindefinite integral is known analytically, and is also analytically invertible to give xas a function of “area under the comparison function to the left of x.” Now pick auniform deviate between 0 and A, where A is the total area under f(x), and use itto get a corresponding x.
Then pick a uniform deviate between 0 and f(x) as the yvalue for the two-dimensional point. You should be able to convince yourself that thepoint (x, y) is uniformly distributed in the area under the comparison function f(x).An equivalent procedure is to pick the second uniform deviate between zeroand one, and accept or reject according to whether it is respectively less than orgreater than the ratio p(x)/f(x).So, to summarize, the rejection method for some given p(x) requires that onefind, once and for all, some reasonably good comparison function f(x).
Thereafter,each deviate generated requires two uniform random deviates, one evaluation of f (toget the coordinate y), and one evaluation of p (to decide whether to accept or rejectthe point x, y). Figure 7.3.1 illustrates the procedure. Then, of course, this proceduremust be repeated, on the average, A times before the final deviate is obtained.Gamma DistributionThe gamma distribution of integer order a > 0 is the waiting time to the athevent in a Poisson random process of unit mean.
For example, when a = 1, it is justthe exponential distribution of §7.2, the waiting time to the first event.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).reject x0f (x)292Chapter 7.Random NumbersA gamma deviate has probability pa (x)dx of occurring with a value betweenx and x + dx, wherepa (x)dx =xa−1 e−xdxΓ(a)x>0(7.3.1)whose inverse indefinite integral is just the tangent function.
It follows that thex-coordinate of an area-uniform random point under the comparison functionf(x) =c01 + (x − x0 )2 /a20(7.3.3)for any constants a0 , c0 , and x0 , can be generated by the prescriptionx = a0 tan(πU ) + x0(7.3.4)where U is a uniform deviate between 0 and 1. Thus, for some specific “bell-shaped”p(x) probability distribution, we need only find constants a0 , c0 , x0 , with the producta0 c0 (which determines the area) as small as possible, such that (7.3.3) is everywheregreater than p(x).Ahrens has done this for the gamma distribution, yielding the followingalgorithm (as described in Knuth [1]):#include <math.h>float gamdev(int ia, long *idum)Returns a deviate distributed as a gamma distribution of integer order ia, i.e., a waiting timeto the iath event in a Poisson process of unit mean, using ran1(idum) as the source ofuniform deviates.{float ran1(long *idum);void nrerror(char error_text[]);int j;float am,e,s,v1,v2,x,y;if (ia < 1) nrerror("Error in routine gamdev");if (ia < 6) {Use direct method, adding waitingx=1.0;times.for (j=1;j<=ia;j++) x *= ran1(idum);x = -log(x);} else {Use 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).To generate deviates of (7.3.1) for small values of a, it is best to add up aexponentially distributed waiting times, i.e., logarithms of uniform deviates.
Характеристики
Тип файла PDF
PDF-формат наиболее широко используется для просмотра любого типа файлов на любом устройстве. В него можно сохранить документ, таблицы, презентацию, текст, чертежи, вычисления, графики и всё остальное, что можно показать на экране любого устройства. Именно его лучше всего использовать для печати.
Например, если Вам нужно распечатать чертёж из автокада, Вы сохраните чертёж на флешку, но будет ли автокад в пункте печати? А если будет, то нужная версия с нужными библиотеками? Именно для этого и нужен формат PDF - в нём точно будет показано верно вне зависимости от того, в какой программе создали PDF-файл и есть ли нужная программа для его просмотра.















