Photography web hosting - 124 RANDOM NUMBERS 3.4.1 (3) The odd-even method,
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