Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C (523184), страница 70
Текст из файла (страница 70)
(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. It7.3 Rejection Method: Gamma, Poisson, Binomial Deviates291Afirst randomdeviate in⌠x⌡0f(x)dxreject x0f(x)f(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.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)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.
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 + y2whose 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.7.3 Rejection Method: Gamma, Poisson, Binomial Deviatesdo {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);}return x;293These 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 is&j+Prob(j) =j−px (m)dm =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:294Chapter 7.Random Numbers1inrejectaccept012345Figure 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.7.3 Rejection Method: Gamma, Poisson, Binomial Deviates295em=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.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
