Web And Email Hosting, Jsp Web Hosting Blog

132 RANDOM NUMBERS 3.4.1 (3) The Poisson distribution with mean p. This distribution is related to the exponential distribution as the binomial distribution is related to the geometric: It represents the number of occurrences, per unit time, of an event that can occur at any instant of time. For example, the number of alpha particles emitted by a radioactive substance in a single second has a Poisson distribution. According to this principle, we can produce a Poisson deviate N by generat- ing independent exponential deviates Xi, Xs, . . . with mean l/p, stopping as soon as Xi + . . . +X, 2 1; then N c m -1. The probability that Xl + . . . + X, > 1 is the probability that a gamma deviate of order m is 2 CL, and this comes to JWm t + eet dt/(m -l)!; hence the probability that N = n is co 1 O tn-le-t & = e-PIli ap tneet dt -n 2 0. s n! (40) If we generate exponential deviates by the logarithm method, the above recipe tells us to stop when -(ln Vi +. . + + In Um)/p 2 1. Simplifying this expression, we see that the desired Poisson deviate can be obtained by calculating e-p, converting it to a fixed point representation, then generating one or more uniform deviates U1, Us, . . . until the product satisfies Vi . . . U, 5 e-p, finally setting Ntm-1. On the average this requires the generation of p + 1 uniform deviates, so it is a very useful approach when fi is not too large. When p is large, we can obtain a method of order log p by using the fact that we know how to handle the gamma and binomial distributions for large orders: First generate X with the gamma distribution of order m = Lap], where Q is a suitable constant. (Since X is equivalent to -ln(Ui . . . Urn), we are essentially bypassing m steps of the previous method.) If X < p, set N c m + Ni, where Ni is a Poisson deviate with mean p -X; and if X > p, set N c Ni, where Ni has the binomial distribution (m -l,p/X). This method is due to J. H. Ahrens and U. Dieter, whose experiments suggest that 6 is a good choice for o. The validity of the above reduction when X 2 p is a consequence of the following important principle: Let Xi, . . . , X, be independent exponential deviates with the same mean; let S, = Xi + . . . + Xj and let Vj = Sj/Sm for 1 5 j 5 m. Then the distribution of VI, VZ, . . . , V,-I is the same as the distribution of m -1 independent uniform deviates sorted into increasing order. To establish this principle formally, we compute the probability that K 5 Ult f I&-I 5 v,-~, given the value of S, = s, for arbitrary values 0 5 Wl 5 .-. 5 21,-l < 1: Let f(vl, ws,. . . , urn-i) be the (m - 1)-fold integral VlS-h/P pevtzlw dtz . . . v,-ls-tl-…-tm–2 J0 pe ,-&Jv2S-tl ,jtmpl ; X J 0 pe-t,-l/~.pe-(S-t’–tm-l)/~ 0

3.4.1 NUMERICAL DISTRIBUTIONS 131 1956), p. 302. To get a random point inside the n-sphere, R. P. Brent suggests taking a point on the surface and multiplying it by U1ln. In three dimensions a significantly simpler method can be used, since each individual coordinate is uniformly distributed between -1 and 1: Find VI, Vi, and S by steps Pl-P3 of Algorithm P; then the desired random point on the surface of a globe is (&VI, aV2,2S -l), where CY= 2J1=s. [Robert E. Knop, CACM 13 (1970), 326.1 F. Important integer-valued distributions. A probability distribution that con- sists only of integer values can essentially be handled by the techniques described at the beginning of this section; but some of these distributions are so important in practice, they deserve special mention here. (1) The geometric distribution. If some event occurs with probability p, the number N of independent trials needed until the event first occurs (or between occurrences of the event) has the geometric distribution. We have N = 1 with probability p, N = 2 with probability (1 -p)p, . . . , N = n with probability (1 -p)+ p. Th is is essentially the situation we have already considered in the gap test of Section 3.3.2; it is also directly related to the number of times certain loops in the algorithms of this section are executed, e.g., steps Pl-P3 of the polar method. A convenient way to generate a variable with this distribution is to set N c [lnU/ ln(l -p)l. (39) To check this formula, we observe that [ln U / ln(1 -p)l = n if and only if n-l < lnU/ln(l-p) 5 n, thatis, (l--~~p)~- > U 2 (l-~p)~, and this happens with probability p(1 -p) - as required. Note that 1nU can be replaced by -Y, where Y has the exponential distribution with mean 1. The special case p = 4 can be handled more easily on a binary computer, since formula (39) becomes N +-I-log, Ul; that is, N is one more than the number of leading zero bits in the binary representation of U. (2) The binomial distribution (t, p). If some event occurs with probability p, and if we carry out t independent trials, the total number N of occurrences equals n with probability (A) pn(l -p) - . (See Section 1.2.10.) In other words if we generate U1, . . . , U,, we want to count how many of these are < p. For small t we can obtain N in exactly this way. For large t, we can generate a beta variate X with integer parameters a and b where a + b - 1 = t; this effectively gives us the bth largest of t elements, without bothering to generate the other elements. Now if X 2 p, we set N c N1 where N1 has the binomial distribution (a -1, p/X), since this tells us how many of u-l random numbers in the range [0,X) are < p; and if X < p, we set N t a + N1 where Nl has the binomial distribution (b -1, (p - X)/(1 -X)), since N1 tells us how many of b- 1 random numbers in the range [X, 1) are < p. By choosing a = 1+ Lt/2], the parameter t will be reduced to a reasonable size after about lg t reductions of this kind. (This approach is due to J. H. Ahrens, who has also suggested an alternative for medium-sized t; see exercise 27.)

130 RANDOM NUMBERS 3.4.1 K + u, la and Yz t U:/b repeatedly until YI +Yz 5 1; then X c Yl/(Yl +Yz). [See M. D. J6hnk, Metrika 8 (1964), 5-15.1 Still another approach, if a and b are integers (not too large), is to set X to the bth largest of a + b - 1 independent uniform deviates (cf. exercise 7 at the beginning of Chapter 5). See also the direct method described by R. C. H. Cheng, CACA4 21 (1978), 317-322. (3) The chi-square distribution with v degrees of freedom (Eq. 3.3.1-22) is obtained by setting X c 2Y, where Y is a random variable having the gamma distribution of order v 12. (4) The F-distribution (variance-ratio distribution) with VI and u2 degrees of freedom is defined by F(x) = v ;li2 @2 r((vl + v2)/2) z tvl/2-l(V2 + vlt)-~1/2–v2/2 dt q4qr(v2/2) s 0 (3;) where x 2 0. Let YI and Y2 be independent, having the chi-square distribution with v1 and 212 degrees of freedom, respectively; set X t YIV~/Y~V~. Or set X c vzY/v~(l -Y), where Y is a beta variate with parameters v1/2, v2/2. (5) The t-distribution with v degrees of freedom is defined by z (1 + t2/v)-(V+1) 2 dt. (37) @r(v/2) s –oo Let YI be a normal deviate (mean 0, variance 1) and let Y2 be independent of YI, having the chi-square distribution with v degrees of freedom; set X t Y~/fl. Alternatively, when v > 2, let Y, be a normal deviate and let Y2 independently have the exponential distribution with mean Z/(Y -2); set 2 +- Y~/(v -2) and reject (Yl,Y2) if e-fieZ 2 1 - 2, otherwise set X t Y~/d(l -2v)(l -z). The latter method is due to George Marsaglia, Math. Comp. 34 (1980), 235-236. [See also A. J. Kinderman, J. F. Monahan, and J. G. Ramage, Math. Comp. 31 (1977), 1009-1018.1 (6) Random point on n-dimensional sphere with radius one. Let X1, X2, . . . , X, be independent normal deviates (mean 0, variance 1); the desired point on the unit sphere is (X,/r, X2/r, . . . , L/r>, where r = XT + X$ + . . . + XZ,. (38) Note that if the X s are calculated using the polar method, Algorithm P, we compute two independent X s each time, and X: + Xg = -2 In S (in the notation of that algorithm); this saves a little of the time needed to evaluate r. The validity of this method comes from the fact that the distribution function for the point (XI, . . . , Xn) has a density that depends only on its distance from the origin, so when it is projected onto the unit sphere it has the uniform distribution. This method was first suggested by G. W. Brown, in Modern Mathematics for the Engineer, First series, ed. by E. F. Beckenbach (New York: McGraw-Hill,

3.4.1 NUMERICAL DISTRIBUTIONS 129 E. Other continuous distributions. Let, us now consider briefly how to handle some other distributions that arise reasonably often in practice. (1) The gamma distribution of order a > 0 is defined by F(x) = & 1 f-1e-t dt, x 2 0. When a = 1, this is the exponential distribution with mean 1; when a = 4, it is the distribution of iZ2, where 2 has the normal distribution (mean 0, variance 1). If X and Y are independent gamma-distributed random variables, of order a and b, respectively, then X + Y has the gamma distribution of order a + b. Thus, for example, the sum of k independent exponential deviates with mean 1 has the gamma distribution of order k. If the logarithm method (32) is being used to generate these exponential deviates, we need compute only one logarithm: X t -ln(U1. . . Uk), where Ul, . . . , Uk are nonzero uniform deviates. This technique handles all integer orders a; to complete the picture, a suitable method for 0 < a < 1 appears in exercise 16. The simple logarithm method is much too slow when a is large, since it requires [al uniform deviates. For large a, the following algorithm due to J. H. Ahrens is reasonably efficient, and it is easy to write in terms of standard subroutines. Algorithm A (Gamma distribution of order a > 1). Al. [Generate candidate.] Set Y e tan(rU), where U is a uniform deviate, and set X t J-Y + a -1. (In place of tan(TU) we could use a polar method, e.g., &/VI in step P4 of Algorithm P.) A2. [Accept?] If X 5 0, return to Al. Otherwise generate a uniform deviate V, and return to Al if V > (1 + Y2) exp((a -1) ln(X/(a -1)) -d-Y>. Otherwise accept X. 1 The average number of times step Al is performed is < 1.902 when a 2 3. For further discussion, proof, and a slightly more complex method that is two to three times faster, see J. H. Ahrens and U. Dieter, Computing 12 (1974), 223-246. There is also an attractive approach for large a based on the remarkable fact that gamma deviates are approximately equal to aX3, where X is normally dis- tributed with mean 1 - 1/(9a) and standard deviation l/a; see G. Marsaglia, Computers and Math. 3 (1977), 321-325. (.2) The beta distribution with positive parameters a and b is defined by r(a+b) F(x) s05 =~ t+l(1 -Q-1 dt, O

128 RANDOM NUMBERS 3.4.1 D. The exponential distribution. After uniform deviates and normal deviates, the next most important random quantity is an exponential deviate. Such numbers occur in arrival time situations; for example, if a radioactive substance emits alpha particles at a rate so that one particle is emitted every p seconds on the average, then the time between two successive emissions has the exponential distribution with mean p. This distribution is defined by the formula F(x) = 1 - e–s/P, x 2 0. (31) (1) Logarithm method. Clearly, if y = F(z) = 1 -e+/p, then z = F- (y) = –I-L ln(1 -y). Therefore by Eq. (7), -p ln(1 -U) has the exponential distribu- tion. Since 1 - U is uniformly distributed when U is, we conclude that X = -plnU (32) is exponentially distributed with mean CL. (The case U = 0 must be avoided.) (2) Random minimization method. We saw in Algorithm F that there are simple and fast alternatives to calculating the logarithm of a uniform deviate. The following especially efficient approach has been developed by G. Marsaglia, M. Sibuya, and J. H. Ahrens. Algorithm S (Exponential distribution with mean CL). This algorithm produces exponential deviates on a binary computer, using uniform deviates with t-bit accuracy. The constants should be precomputed, extending until Q[k] > 1 - 21mt. Sl. [Get U and shift.] Generate a t-bit uniform random binary fraction U = (.blbz . . . bt)2; locate the first zero bit bj, and shift off the leading j bits, setting U + (.bj+l . . . bt)2. (As in Algorithm F, the average value of j is 2.) 52. [Immediate acceptance?] If U < In 2, set X t p(j In 2 + U) and terminate the algorithm. (Note that Q[l] = ln2.) 53. [Minimize. ] Find the least k 2 2 such that U < Q[k]. Generate k new uniform deviates VI, . . . , uk and set V t min(Ui, . . . , uk). S4. [Deliver the answer.] Set X c p(j + V) In 2. 1 Alternative ways to generate exponential deviates (e.g., a ratio of uniforms as in Algorithm R) might also be used.

3.4.1 NUMERICAL DISTRIBUTIONS 127 for some nonnegative integrable function g. If we set X c V/U, the probability that X 5 x can be calculated by integrating du dv over the region defined by the two relations in (26) plus the auxiliary condition v/u 5 x, then dividing by the same integral without this extra condition. Letting w = tu so that dv = udt, the integral becomes z dtm 15 udu=-2 J_ g(t) dt. s–cc I 0 co Hence the probability that X 5 x is /-~mS(tPt /~~+J7W. (27) The normal distribution comes out when g(t) = e-t2/2; and the condition u2 5 g(v/u) simplifies in this case to (w/u) 5 -4 lnu. It is easy to see that the set of all (u, v) satisfying this relation is entirely contained in the rectangle of Fig. 13. The bounds in steps R2 and R3 define interior and exterior regions with simpler boundary equations. The well-known inequality which holds for all real numbers x, can be used to show that l+lnc-cu 5 -1nu 5 l/cu-l+lnc (28) for any constant c > 0. Exercise 21 proves that c = e1j4 is the best possible constant to use in step R2. The situation is more complicated in step R3, and there doesn t seem to be a simple expression for the optimum c in that case, but computational experiments show that the best value for R3 is approximately e1.35. The approximating curves (28) are tangent to the true boundary when u = l/c. It is possible to obtain a faster method by partitioning the region into subregions, most of which can be handled more quickly. Of course, this means that auxiliary tables will be needed, as in Algorithms M and F. (5) Variations of the normal distribution. So far we have considered the normal distribution with mean zero and standard deviation one. If X has this distribu- tion, then y=p++x (2% has the normal distribution with mean p and standard deviation u. Furthermore, if Xl and Xz are independent normal deviates with mean zero and standard deviation one, and if K =Pl+alXl, 35 = P2 + c72(PXl + &=-7X2), (30) then Yi and Y2 are dependent random variables, normally distributed with means ~1, ~2 and standard deviations ui, us, and with correlation coefficient p. (For a generalization to rz variables, see exercise 13.)

126 RANDOM NUMBERS 3.4.1 Fig. 13. Region of acceptance in the ratio -of -uniforms method for normal deviates. Lengths of lines with coor- dinate ratio x have the normal distribu- tion. - (0, - 42/e) Exercises 20 and 21 work out the timing analysis; four different algorithms are analyzed, since steps R2 and R3 can be included or omitted depending on one s preference. The following table shows how many times each step will be performed, on the average, depending on which of the optional tests is applied: Step Neither R2 only R3 only Both Rl 1.369 1.369 1.369 1.369 R2 0 1.369 0 1.369 (25) R3 0 0 1.369 0.467 R4 1.369 0.467 1.134 0.232 Thus it pays to omit the optional tests if there is a very fast logarithm operation, but if the log routine is rather slow it pays to include them. But why does it work? One reason is that we can calculate the probability that X < Z, and it turns out to be the correct value (10). But such a calculation isn t very easy unless one happens to hit on the right trick, and anyway it is better to understand how the algorithm might have been discovered in the first place. Kinderman and Monahan derived it by working out the following theory that can be used with any well-behaved density function f(x) [cf. ACM Trans. Math. Software 3 (1977), 257-2601. In general, suppose that a point (U, V) has been generated uniformly over the region of the (u, v)-plane defined by

3.4.1 NUMERICAL DISTRIBUTIONS 125 the rejection method described above, with h(z) = x2/2–a2/2 = y2/2+ay where y = 2 - a. Exercise 12 proves that h(x) 2 1 as required in (21).) Set Y c dj times (.bj+i . . . bt)s and V e ($Y + a)Y. (Since the average value of j is 2, there will usually be enough significant bits in (.bj+i . . . bt)s to provide decent accuracy. The calculations are readily done in fixed point arithmetic.) F4. [Reject?] Generate a uniform deviate U. If V < U, go on to step F5. Otherwise set V to a new uniform deviate; and if now U < V (i.e., if K is even, in the discussion above), go back to F3, otherwise repeat step F4. F5. [Return X.1 Set X +a+Y. If@=l,setXc-X. 1 Values of dj for 1 2 j 5 47 appear in a paper by Ahrens and Dieter, Math. Camp. 27 (1973), 927-937; their paper discusses refinements of the algorithm that improve its speed at the expense of more tables. Algorithm F is attractive since it is almost as fast as Algorithm M and it is easier to implement. The average number of uniform deviates per normal deviate is 2.53947; R. P. Brent [CACM 17 (1974), 704-7051 has shown how to reduce this number to 1.37446 at the expense of two subtractions and one division per uniform deviate saved. (4) Ratios of uniform deviates. There is yet another good way to generate normal deviates, discovered by A. J. Kinderman and J. F. Monahan in 1976. Their idea is to generate a random point (U, V) in the region defined by O

124 RANDOM NUMBERS 3.4.1 (3) The odd-even method, due to G. E. Forsythe. An amazingly simple technique for generating random deviates with a density of the general exponential form f(x) = Ce-+ ), for a 5 z < b, f(z) = 0 otherwise, ( JO> when 0 2 h(z) 5 1 for a 5 2 < b, (21) was discovered by John von Neumann and G. E. Forsythe about 1950. The idea is based on the rejection method described earlier, letting g(z) be the uniform distribution on [a, b): We set X c a+(&-a)U, where U is a uniform deviate, and then we want to accept X with probability e- cX). The latter operation could be done by testing e--h(X) vs. V, or h(X) vs. -In V, when V is another uniform deviate, but the job can be done without applying any transcendental functions in the following interesting way. Set VO t h(X), then generate uniform deviates K, E?, **a until finding some K 2 1 with VK-~ < VK. For fixed X and k, the probability that h(X) 2 VI 2 .. . 2 Vj is l/k! times the probability that max(V1, . .* tvk) 2 h(X), namely h(X)k/k!; hence the probability that K = k is h(X) - /(k -l)! -h(X) /k!, and the probability that K is odd is (22) Therefore we reject X and try again if K is even; we accept X as a random variable with density (20) if K is odd. Note that we usually won t have to gen- erate many V s in order to determine K, since the average value of K (given X) is Ck2,,probability that(K > k) = Ck,,, h(X) /k! = ehcX) 2 e. Forsythe realized some years later that this approach leads to an efficient method for calculating normal deviates, without the need for any auxiliary routines to calculate square roots or logarithms as in Algorithms P and M. His procedure, with an improved choice of intervals [a, b) due to J. H. Ahrens and U. Dieter, can be summarized as follows. Algorithm F (Odd-even method for normal deviates). This algorithm generates normal deviates on a binary computer, assuming approximately t + 1 bits of accuracy. The algorithm requires a table of values dj = aj -u+~, for 1 5 j <_ t + 1, where aj is defined by the relation e-x2/2dz = _. 1 2-7 (23) Fl. [Get U.] Generate a uniform random number U = (.bobl . . . bt)2, where bo, . 7 bt denote the bits in binary notation. Set + +- bo, j c 1, and a t 0. F2. [Find first zero bj.] If bj = 1, set a + a + dj, j + j + 1, and repeat this step. (If j = t + 1, treat bj as zero.) F3. [Generate candidate.] (Now a = aj-1, and the current value of j occurs with probability NN 2-j. We will generate X in the range [a+~, a?), using

3.4.1 NUMERICAL DISTRIBUTIONS 123 Table 1 EXAMPLE OF TABLES USED WITH ALGORITHM M* j p3 P,+ls QJ 5 y3+= 2, zj+l6 %+1 L&+15 J-$+15 0 ,000 ,067 0.00 0.59 0.20 0.21 0.0 1 ,849 ,161 ,236 -0.92 0.96 1.32 0.24 0.2 ,505 25.00 2 ,970 ,236 ,206 -5.86 -0.06 6.66 0.26 0.4 ,773 12.50 3 .855 .285 .234 -0.58 0.12 1.38 0.28 0.6 ,876 8.33 4 ,994 ,308 ,201 -33.13 1.31 34.93 0.29 0.8 ,939 6.25 5 ,995 ,304 ,201 -39.55 0.31 41.35 0.29 1.0 ,986 5.00 6 ,933 ,280 ,214 -2.57 1.12 2.97 0.28 1.2 ,995 4.06 7 ,923 ,241 ,217 -1.61 0.54 2.61 0.26 1.4 ,987 3.37 8 ,727 ,197 ,275 0.67 0.75 0.73 0.25 1.6 ,979 2.86 9 1.000 ,152 ,200 0.00 0.56 0.00 0.24 1.8 ,972 2.47 10 ,691 ,112 ,289 0.35 0.17 0.65 0.23 2.0 ,966 2.16 11 ,454 ,079 ,440 -0.17 0.38 0.37 0.22 2.2 ,960 1.92 12 ,287 ,052 ,698 0.92 -0.01 0.28 0.21 2.4 ,954 1.71 13 ,174 ,033 1.150 0.36 0.39 0.24 0.21 2.6 ,948 1.54 14 ,101 ,020 1.974 -0.02 0.20 0.22 0.20 2.8 ,942 1.40 15 ,057 ,086 3.526 0.19 0.78 0.21 0.22 3.0 ,936 1.27 *In practice, this data would be given with much greater precision; the table shows only enough figures so that interested readers may test their own algorithms for computing the values more accurately. M3. [Wedge or tail?] (Now 15 < j 5 31, and each particular value j occurs with probability p3.) If j = 31, go to M7. M4. [Get U 2 V.] Generate two new uniform deviates, U, V; if U > V, exchange U * V. (We are now performing Algorithm L.) Set X c S,-is+ 5 U. M5. [Easy case?] If V 5 Dj, go to M9. M6. [Another try?] If V > U + E~(e(s~~~~-xz)~z -l), go back to step M4; otherwise go to M9. (This step is executed with low probability.) M7. [Get supertail deviate.] Generate two new independent uniform deviates, U, V,andsetXcdS-21nV. M8. [Reject?] If UX 2 3, go back to step M7. (This will occur only about one-twelfth as often as we reach step M8.) M9. [Attach sign.] If q!~ = 1, set X e -X. 1 This algorithm is a very pretty example of mathematical theory intimately interwoven with programming ingenuity-a fine illustration of the art of com- puter programming! Only steps Ml, M2, and M9 need to be performed most of the time, and the other steps aren t terribly slow either. The first publications of the rectangle-wedge-tail method were by G. Marsaglia, Annals Math. Stat. 32 (1961), 894-899; G. Marsaglia, M. D. MacLaren, and T. A. Bray, CACA4 7 (1964), 4-10. Further refinements of Algorithm M have been developed by G. Marsaglia, K. Ananthanarayanan, and N. J. Paul, hf. Proc. Letters 5 (1976), 27-m.