J2ee Web Hosting, Perl Scripting And Developing Blog

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.

122 RANDOM NUMBERS 3.4.1 Ml. Get U + No Wedge M3. Wedge or tail? M4. Get U Pj, set X + Yj + fZj and go to M9. Otherwise if j 5 15 (i.e., bl = 0), set X t Sj + f&j and go to M9. (This is an adaptation of Walker s alias method (3).)

3.4.1 NUMERICAL DISTRIBUTIONS 121 Fig. 11. Region of acceptance in Algorithm L. cases f(x)/cg(x) is always 0 or 1; then U need not be generated. In other cases if f(x)/cg(Z) is hard to compute, we may be able to squeeze it between two bounding functions r(x) I f(x)/cdx) 5 4×1 (18) that are much simpler, and the exact value of f(z)/cg(z) need not be calculated unless T(Z) 5 U < s(z). The following algorithm solves the wedge problem by developing the rejection method still further. Algorithm L (Nearly linear densities). This algorithm may be used to generate a random variable X for any distribution whose density f(x) satisfies the following conditions (cf. Fig. 10): f(x) = 0, for x < s and for x > s + h; (1% a -b(x -s)/h 5 f(x) 5 b - b(x -s)/h, for s < x 5 s + h. Ll. [Get U 5 V.] Generate two independent random variables U, V, uniformly distributed between zero and one. If U > V, exchange U ++ V. L2. [Easy case?] If V 5 a/b, go to L4. L3. [Try again?] If V > U + (l/b)f(s + hU), go back to step Ll. (If a/b is close to 1, this step of the algorithm will not be necessary very often.) L4. [Compute X.1 Set X c s + hU. I When step L4 is reached, the point (U, V) is a random point in the area shaded in Fig. 11, namely, 0 5 U < V < U + (l/b)f(s + hU). Conditions (19) ensure that ; 5 U + $I(s + hU) 5 1. Now the probability that X 5 s + hx, for 0 5 x 5 1, is the ratio of area to the left of the vertical line U = x in Fig. 11 to the total area, namely, oz f f(s + hu) du / 6 ;f(s + hu) du = 1 f(v) dv; .I therefore X has the correct distribution.

120 RANDOM NUMBERS 3.4.1 b . . a .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 ~ S s+h Fig. 10. Density functions for which Algorithm L may be used to generate random numbers. In order to generate such rectangular portions of the distribution, we simply compute x=QJ+s, (16) where U is uniform and S takes the value (J -1)/5 with probability pj. Since pl + . . . + p15 = .9183, we can use simple uniform deviates like this about 92 percent of the time. In the remaining 8 percent, we will usually have to generate one of the wedge-shaped distributions Frs, . . . , Fso. Typical examples of what we need to do are shown in Fig. 10. When z < 1, the curved part is concave downward, and when 2 > 1 it is concave upward, but in each case the curved part is reasonably close to a straight line, and it can be enclosed in two parallel lines as shown. To handle these wedge-shaped distributions, we will rely on yet another general technique, von Neumann s rejection method for obtaining a complicated density from another one that encloses it. The polar method described above is a simple example of such an approach: Steps Pl-P3 obtain a random point inside the unit circle by first generating a random point in a larger square, rejecting it and starting over again if the point was outside the circle. The general rejection method is even more powerful than this. To generate a random variable X with density f, let g be another probability density function such that f(t) 2 cdt) (17) for all t, where c is a constant. Now generate X according to density g, and also generate an independent uniform deviate U. If U 2 f(X)/cg(X), reject X and start again with another X and U. When the condition U < f(X)/cg(X) finally occurs, the resulting X will have density f as desired. [Proof: X 2 2 will occur with probability p(z) = Jz,(g(t)dt . f(t)/cg(t)) + qp(z), where the quantity 4 = r&w * (1 - fwm)) = 1 - l/ c is the probability of rejection; hence p(z) = sz, f(t) dt.] The rejection technique is most efficient when c is small, since there will be c iterations on the average before a value is accepted. (See exercise 6.) In some

3.4.1 NUMERICAL DISTRIBUTIONS 119 Fig. 9. The density function divided into 31 parts. The area of each part represents the average number of times a random number with that density is to be computed. probability pj is very small in this ease. Most of the distributions F,(s) will be quite easy to accommodate, since they will be trivial modifications of the uniform distribution. The resulting method yields an extremely efficient program, since its average running time is very small. It is easier to understand the method if we work with the derivatives of the distributions instead of the distributions themselves. Let f(x) = F (x), fj(X) = Q (x); these are called the density functions of the probability distributions. Equation (13) b ecomes f(z) = Plflcz) + Pz.fi(Z) + . . . + Pnfn(X). (14) Each f?(z) is > 0, and the total area under the graph of f,(z) is 1; so there is a convenient graphical way to display the relation (14): The area under f(z) is divided into n parts, with the part corresponding to fj(z) having area pj. See Fig. 9, which illustrates the situation in the case of interest to us here, with f(x) = F (x) = me-z2/2; the area under this curve has been divided into n = 31 parts. There are 15 rectangles, which represent plfl(z), . . . , p15fl5(z); there are 15 wedge-shaped pieces, which represent p,sf,s(s), . . . , psofso(z); :md the remaining part pslfsl(z) is the tail, namely the entire graph of f(z) for x 2 3. The rectangular parts fl(x), . . . , f15 ( X) represent uniform distributions. For example, fs(x) represents a random variable uniformly distributed between 3 and 3. The altitude of pjf,(x) is f(j/5), hence the area of the jth rectangle is pj = if(j/5) = & e–j2/50, for 1 < j 5 15. (15) i

118 RANDOM NUMBERS 3.4.1 To prove the validity of this method, we use elementary analytic geometry and calculus: If 5 < 1 in step P3, the point in the plane with Cartesian coordinates (VI, VZ) is a random point uniformly distributed inside the unit circle. Transforming to polar coordinates VI = R cos 0, V2 = R sin 0, we find S = R2, Xi = &%n??cos 0, Xz = dm sin 0. Using also the polar coordinates X1 = R cos Q , X2 = R sin O , we find that 0 = 0 and R = dn. It is clear that R and 0 are independent, since R and 0 are independent inside the unit circle. Also, 0 is uniformly distributed between 0 and 215 and the probability that R < T is the probability that -2lnS < r2, i.e., the probability that S 2 ~ 27~. This equals 1 -e-r2/2, since S = E2 is uniformly distributed between zero and one. The probability that R lies between T and r+dr is therefore the derivative of l-eerai2, namely, ~e+ /~ dr. Similarly, the probability that 0 lies between 6 and 0 + d0 is (1/27r) &9. The joint probability that Xr 5 x1 and that Xz 2 52 now can be computed, it is -e --72/2 r dr de J 1 {(r,e)Ircoselsl,rsinB~zz}2~ 1 dz dy =---- ,-(~2+Y2)/2 27r J{(~,Y)I~l~l,Yl~Z) = 1 21 -x=/2 da: 1 x2 -Y=/2& . V-J5 -me WJ -me > 277 This proves that Xi and X2 are independent and normally distributed, as desired. (L?) The rectangle-wedge-tail method, introduced by G. Marsaglia. In this method we use the distribution (12) so that F(z) gives the distribution of the absolute value of a normal deviate. After X has been computed according to distribution (12), we will attach a random sign to its value, and this will make it a true normal deviate. The rectangle-wedge-tail approach is based on several important general techniques that we shall explore as we develop the algorithm. The first key idea is to regard F(z) as a mixture of several other functions, namely to write F(z) = PlFl(Z) + P2J72(Z) + . . . + P?zFn(Z>, (13) where Fi, F., . . . . F, are appropriate distributions and pl, pz, . . . , p, are nonnegative probabilities that sum to 1. If we generate a random variable X by choosing distribution Fj with probability pj, it is easy to see that X will have distribution F overall. Some of the distributions Fj(Z) may be rather difficult to handle, even harder than F itself, but we can usually arrange things so that the

3.4.1 NUMERICAL DISTRIBUTIONS 117 max(Ul, UZ,. . . , Ut) has the distribution function F(z) = xt, for 0 < 2 5 1. This is the basis of the maximum-of-t test given in Section 3.3.2. Note that the inverse function in this case is F-l(y) = fl. In the special case t = 2, we see therefore that the two formulas x=&7 and X = max( U1, UZ) (9) will give equivalent distributions to the random variable X, although this is not obvious at first glance. We need not take the square root of a uniform deviate. The number of tricks like this is endless: any algorithm that employs random numbers as input will give a random quantity with some distribution as output. The problem is to find general methods for constructing the algorithm, given the distribution function of the output. Instead of discussing such methods in purely abstract terms, we shall study how they can be applied in important cases. C. The normal distribution. Perhaps the most important nonuniform, continuous distribution is the so-called normal distribution with mean zero and standard deviation one: z F(z) = -!—t2/2 dt. (10) & s -tee The significance of this distribution was indicated in Section 1.2.10. Note that the inverse function F-l is not especially easy to compute; but we shall see that several other techniques are available. (1) The polar method, due to G. E. P. Box, M. E. Muller, and G. Marsaglia. (See Annals Math. Stat. 28 (1958), 610; and Boeing Scientific Res. Lab. report Dl-82-0203 (1962).) Algorithm P (Polar method for normal deviates). This algorithm calculates two independent normally distributed variables, X1 and X2. Pl. [Get uniform variables.] Generate two independent random variables, UI, U2, uniformly distributed between zero and one. Set VI c 2Ul -1, V2 c 2U2 -1. (Now VI and V2 are uniformly distributed between -1 and +l. On most computers it will be preferable to have VI and V2 represented in floating point form at this point.) P2. [Compute S.] Set S t VT + Vi. P3. [Is S 2 I?] If S 2 1, return to step Pl. (Steps Pl through P3 are executed 1.27 times on the average, with a standard deviation of 0.587; see exercise 6.) P4. [Compute Xl, X2.1 Set X1, X2 according to the following two equations: -21nS -21nS X2=& (11) Xl=v1 -, r ~ /—– s S These are the normally distributed variables desired. 1

116 RANDOM NUMBERS 3.4.1 We can do this using (3), if k = 16 and x3 = j for 0 5 j < 16, and if the P and Y tables are set up as follows: pj=oo 4 8 1~111~~~~000 Yj=59 7 4 * 6 * * * 8 4 7 10 6 7 8 (When Pj = 1, y3 is not used.) For example, the value 7 occurs with probability & . ((1 -Pz) + P7 + (1 -PII) + (1 -PI,)) = & as required. It is a peculiar way to throw dice, but the results are indistinguishable from the real thing. B. General methods for continuous distributions. The most general real-valued distribution may be expressed in terms of its distribution function F(x), which specifies the probability that a random quantity X will not exceed x: F(z) = probability that (X 5 x). (4 This function always increases monotonically from zero to one; i.e., F(G) I F(x2), if xl I x2; F(-co) = 0, F(+cm) = 1. (5) Examples of distribution functions are given in Section 3.3.1, Fig. 3. If F(x) is continuous and strictly increasing (so that F(xl) < F(Q) when x1 < x2), it takes on all values between zero and one, and there is an inverse function F-l(y) such that, for 0 < y < 1, Y = F(x) if and only if x = F- (y). (6) In general we can compute a random quantity X with the continuous, strictly increasing distribution F(z) by setting x = F- (U), (7) where U is uniform; this works because the probability that X 5 x is the prob- ability that F-l(U) 2 x, i.e., the probability that U 5 F(x), and this is F(x). The problem now reduces to one of numerical analysis, namely to find good methods for evaluating F- (U) to the desired accuracy. Numerical analysis lies outside the scope of this seminumerical book; yet there are a number of important shortcuts available to speed up this general approach, and we will consider them here. In the first place, if Xr is a random variable having the distribution Fl(x) and if Xz is an independent random variable with the distribution Fz(x), then m&Xl, X2) has the distribution J ~(x)~z(x), min(Xl, X2) has the distribution Fl(x)+Fz(x)–l(x)F2(~). (8) (See exercise 4.) For example, a uniform deviate U has the distribution F(x) = x, for 0 5 x 5 1; if VI, U2, . . . . Ut are independent uniform deviates, then