Web And Email Hosting, Jsp Web Hosting Blog

142 RANDOM NUMBERS 3.5 3.5. WHAT IS A RANDOM SEQUENCE? A. Introductory remarks. We have seen in this chapter how to generate sequences (1) of real numbers in the range 0 2 U, < 1, and we have called them random sequences even though they are completely deterministic in character. To justify this terminology, we claimed that the numbers behave as if they are truly random. Such a statement may be satisfactory for practical purposes (at the present time), but it sidesteps a very important philosophical and theoretical question: Precisely what do we mean by random behavior ? A quantitative definition is needed. It is undesirable to talk about concepts that we do not really understand, especially since many apparently paradoxical statements can be made about random numbers. The mathematical theory of probability and statistics carefully sidesteps the question; it refrains from making absolute statements, and instead expresses everything in terms of how much probability is to be attached to statements involving random sequences of events. The axioms of probability theory are set up so that abstract probabilities can be computed readily, but nothing is said about what probability really signifies, or how this concept can be applied meaningfully to the actual world. In the book Probability, Statistics, and Truth (New York: Macmillan, 1957), R. von Mises discusses this situation in detail, and presents the view that a proper definition of probability depends on obtaining a proper definition of a random sequence. Let us paraphrase here some statements made by two of the many authors who have commented on the subject. D. H. Lehmer (1951): A random sequence is a vague notion embodying the idea of a sequence in which each term is unpredictable to the uninitiated and whose digits pass a certain number of tests, traditional with statisticians and depending somewhat on the uses to which the sequence is to be put. J. IV. Franklin (1962): The sequence (1) is random if it has every property that is shared by all infinite sequences of independent samples of random variables from the uniform distribution. Franklin s statement essentially generalizes Lehmer s to say that the se- quence must satisfy all statistical tests. His definition is not completely precise, and we will see later that a reasonable interpretation of his statement leads us to conclude that there is no such thing as a random sequence! So let us begin with Lehmer s less restrictive statement and attempt to make it precise. What we really want is a relatively short list of mathematical properties, each of which is satisfied by our intuitive notion of a random sequence; furthermore, the list is to be complete enough so that we are willing to agree that any sequence satisfying these properties is random. In this section, we will develop what seems to be

3.42 RANDOM SAMPLING AND SHUFFLING 141 12. [M,E6] The gist of Algorithm P is that any permutation T can be uniquely written as a product of transpositions in the form T = (att) . . . (as3)(~2), where 1 2 a3 2 j for t 2 j > 1. Prove that there is also a unique representation of the form r = (b22)(b33). . . (btt), where 1 5 bj 2 j for 1 < j 5 t, and design an algorithm that computes the b s from the u s in O(t) steps. 13. [MB] (S. W. Golomb.) One of the most common ways to shuffle cards is to divide the deck into two parts as equal as possible, and to riffle them together. (See the discussion of card-playing etiquette in Hoyle s rules of card games; we read, A shuffle of this sort should be made about three times to mix the cards thoroughly. ) Consider adeckof2n-1 cardsXl,Xz, . . . . X2,-1; a perfect shuffle s divides this deck into X1, X2, . . . , X, and Xn+l, . . . , X2n--1, then perfectly interleaves them to obtain Xl, xn+1, x2, -%x+2, . . . , Xzn-l, X,. The cut operation c3 changes X1, X2, . . , Xzn-l into X3+1, . . . , X2n-1, Xl, ...I Xj. Show that by combining perfect shuffles and cuts, at most (2n -1)(2n -2) different arrangements of the deck are possible, if n > 1. F 14. [so] (Ole-Johan Dahl.) If Xk = k for 1 5 k 5 t at the start of Algorithm P, and if we terminate the algorithm when j reaches the value t -n, the sequence Xten+l, . . . , Xt is a random permutation of a random combination of n elements. Show how to simulate the effect of this procedure using only O(n) cells of memory. b 15. [Mz~ ] Devise a way to compute a random sample of n records from N, given N and n, based on the idea of hashing (Section 6.4). Your method should use O(n) storage locations and an average of O(n) units of time, and it should present the sample as a sorted set of integers 1 5 Xl < X2 < … < X, 5 N. 16. [.%I Discuss the problem of weighted sampling, where each subset of n elements is obtained with probability proportional to the product of the weights of the elements.

140 RANDOM NUMBERS 3.4.2 This algorithm was first published by L. E. Moses and R. V. Oakford, in Tables of Random Permutations (Stanford University Press, 1963); and by R. Durstenfeld, CACM 7 (1964), 420. It can also be modified to obtain a random permutation of a random combination (see exercise 14). For a discussion of random combinatorial objects of other kinds (e.g., parti- tions), see Section 7.2 and/or the book Combinatorial Algorithms by Nijenhuis and Wilf (New York: Academic Press, 1975). EXERCISES 1. [Mz?] Explain Eq. (1). 2. [.20] Prove that Algorithm S never tries to read more than N records of its input file. b 3. [,% I The (t+l)st item in Algorithm S is selected with probability (n–m)/(N–t), not n/N, yet the text claims that the sample is unbiased-so each item should be selected with the same probability. How can both of these statements be true? 4. [A&B] Let p(m, t) be the probability that exactly m items are selected from among the first t in the selection sampling technique. Show directly from Algorithm S that 5. [M.24] What is the average value of t when Algorithm S terminates? (In other words, how many of the N records have been passed, on the average, before the sample is complete?) 6. [M2 4] What is the standard deviation of the value computed in the previous exercise? b 7. [Mz~] Prove that any given choice of n records from the set of N is obtained by Algorithm S with probability l/(E). Th ere f ore the sample is completely unbiased. 8. [A4461 Algorithm S computes one uniform deviate for each input record it handles. Find a more efficient way to determine the number of input records to skip before the first is selected, assuming that N/n is rather large. (We could iterate this process to select the remaining n -1 records, thus reducing the number of necessary random deviates from order N to order n.) 9. [1,6] Let n = 3. If Algorithm R is applied to a file containing 20 records numbered 1 thru 20, and if the random numbers generated in step R3 are respectively 4,1,6,7,5,3,5,11,11,3,7,9,3,11,4,5,4, which records go into the reservoir? Which are in the final sample? 10. [15] Modify Algorithm R so that the reservoir is eliminated, assuming that the n records of the current sample can be held in memory. b 11. [M%] Let pm be the probability that exactly m elements are put into the reservoir during the first pass of Algorithm R. Determine the generating function G(z) = C,,, pmzm, and find the mean and standard deviation. (Use the ideas of Section 1.2.10.)

3.4.2 RANDOM SAMPLING AND SHUFFLING 139 Algorithm R is due to Alan G. Waterman. The reader may wish to work out the example of its operation that appears in exercise 9. If the records are sufficiently short, it is of course unnecessary to have a reservoir at all; we can keep the n records of the current sample in memory at all times (see exercise 10). The natural question to ask about Algorithm R is, What is the expected size of the reservoir? Exercise 11 shows that the average value of m is exactly n(l + H,~J - H,); this is approximately n(1 + ln(iV/n)). So if N/n = 1000, the reservoir will contain only about & as many items as the original file. Note that Algorithms S and R can be used to obtain samples for several independent categories simultaneously. For example, if we have a large file of names and addresses of U.S. residents, we could pick random samples of exactly 10 people from each of the 50 states without making 50 passes through the file, and without first sorting the file by state. The sampling problem can be regarded as the computation of a random combination, according to the conventional definition of combinations of N things taken n at a time (see Section 1.2.6). Now let us consider the problem of computing a random permutation of t objects; we will call this the shuffling problem, since shuffling a deck of cards is nothing more than subjecting it to a random permutation. A moment s reflection is enough to convince oneself that the approaches people traditionally use to shuffle cards are miserably inadequate; there is no hope of obtaining each of the t! permutations with anywhere near equal probability by such methods. It has been said that expert bridge players make use of this fact when deciding whether or not to finesse. If t is small, we can obtain a random permutation very quickly by generating a random integer between 1 and t!. For example, when t = 4, a random number between 1 and 24 suffices to select a random permutation from a table of all possibilities. But for large t, it is necessary to be more careful if we want to claim that each permutation is equally likely, since t! is much larger than the accuracy of individual random numbers. A suitable shuffling procedure can be obtained by recalling Algorithm 3.3.2P, which gives a simple correspondence between each of the t! possible permutations and a sequence of numbers (cl, ~2,. . . , c+~), with 0 5 c3 5 j. It is easy to compute such a set of numbers at random, and we can use the correspondence to produce a random permutation. Algorithm P (Shuffling). Let X1, X2, . . . , X, be a set oft numbers to be shuffled. Pl. [Initialize.] Set j c t. P2. [Generate V.] Generate a random number U, uniformly distributed between zero and one. P3. [Exchange.] Set k t [jU] + 1. (N ow k is a random integer, between 1 and j.) Exchange X, c) X,. P4. [Decrease j.] Decrease j by 1. If j > 1, return to step P2. 1

138 RANDOM NUMBERS 3.4.2 A problem arises if we don t know the value of N in advance, since the precise value of N is crucial in Algorithm S. Suppose we want to select n items at random from a file, without knowing exactly how many are present in that file. We could first go through and count the records, then take a second pass to select them; but it is generally better to sample m 2 n of the original items on the first pass, where m is much less than N, so that only m items must be considered on the second pass. The trick, of course, is to do this in such a way that the final result is a truly random sample of the original file. Since we don t know when the input is going to end, we must keep track of a random sample of the input records seen so far, thus always being prepared for the end. As we read the input we will construct a reservoir that contains only those m records that have appeared among the previous samples. The first n records always go into the reservoir. When the (t + 1)st record is being input, for t 2 n, we will have in memory a table of n indices pointing to those records in the reservoir that belong to the random sample we have chosen from the first t records. The problem is to maintain this situation with t increased by one, namely to find a new random sample from among the t + 1 records now known to be present. It is not hard to see that we should include the new record in the new sample with probability n/(t + l), and in such a case it should replace a random element of the previous sample. Thus, the following procedure does the job: Algorithm R (Reservoir sampling). To select n records at random from a file of unknown size 2 n, given n > 0. An auxiliary file called the reservoir contains all records that are candidates for the final sample. The algorithm uses a table of distinct indices 1[j] for 1 < j 5 n, each of which points to one of the records in the reservoir. Rl. [Initialize.] Input the first n records and copy them to the reservoir. Set 1[j] + j for 1 < j 5 n, and set t t m c n. (If the file being sampled has fewer than n records, it will of course be necessary to abort the algorithm and report failure. During this algorithm, {1[1], . . . , I[n]} point to the records in the current sample, m is the size of the reservoir, and t is the number of input records dealt with so far.) R2. [End of file?] If there are no more records to be input, go to step R6. R3. [Generate and test.] Increase t by 1, then generate a random integer M between 1 and t (inclusive). If M > n, go to R5. R4. [Add to reservoir.] Copy the next record of the input file to the reservoir, increase m by 1, and set I[M] c m. (The record previously pointed to by I[M] is being replaced in the sample by the new record.) Go back to R2. R5. [Skip.] Skip over the next record of the input file (do not include it in the reservoir), and return to step R2. R6. [Second pass.] Sort the I table entries so that 1[1] < ..-< I[n]; then go through the reservoir, copying the records with these indices into the output file that is to hold the final sample. 1

3.4.2 RANDOM SAMPLING AND SHUFFLING 137 of these select the (t + 1)st element. The idea developed in the preceding paragraph leads immediately to the following algorithm: Algorithm S (Selection sampling technique). To select n records at random from a set of N, where 0 < n 5 N. Sl. [Initialize.] Set t t 0, m c 0. (During this algorithm, m represents the number of records selected so far, and t is the total number of input records we have dealt with.) S2. [Generate U.] Generate a random number U, uniformly distributed between zero and one. S3. [Test.] If (N -t)U 2 n -m, go to step S5. S4. [Select.] Select the next record for the sample, and increase m and t by 1. If m < n, go to step S2; otherwise the sample is complete and the algorithm terminates. S5. [Skip.] Skip the next record (do not include it in the sample), increase t by 1, and go to step S2. 1 This algorithm may appear to be unreliable at first glance and, in fact, to be incorrect; but a careful analysis (see the exercises below) shows that it is completely trustworthy. It is not difficult to verify that a) At most N records are input (we never run off the end of the file before choosing n items). b) The sample is completely unbiased; in particular, the probability that any given element is selected, e.g., the last element of the file, is n/N. Statement (b) is true in spite of the fact that we are not selecting the (t+l)st item with probability n/N, we select it with the probability in Eq. (l)! This has caused some confusion in the published literature. Can the reader explain this seeming contradiction? (Note: When using Algorithm S, one should be careful to use a different source of random numbers U each time the program is run, to avoid connections between the samples obtained on different days. This can be done, for example, by choosing a different value of X0 for the linear congruential method each time; X0 could be set to the current date, or to the last X value generated on the previous run of the program.) We will usually not have to pass over all N records; in fact, since (b) above says that the last record is selected with probability n/N, we will terminate the algorithm before considering the last record exactly (1 - n/N) of the time. The average number of records considered when n = 2 is about 3N, and the general formulas are given in exercises 5 and 6. Algorithm S and a number of other sampling techniques are discussed in a paper by C. T. Fan, Mervin E. Muller, and Ivan Rezucha, J. Amer. Stat. Assoc. 57 (1962), 387-402. The method was independently discovered by T. G. Jones, CACM 5 (1962), 343.

136 RANDOM NUMBERS 3.4.1 25. [AL%] Let XI, X2, . . . , X5 be binary words each of whose bits is independ- ently 0 or 1 with probability f. What is the probability that a given bit position of Xl V (X2 A (Xs V (X4 A X5))) contains a l? Generalize. 26. [MI81 Let Nl and NZ be independent Poisson deviates with respective means ~1 and ~2, where PI > ~2 2 0. Prove or disprove: (a) NI + NZ has the Poisson distribution with mean ~1 + ~2. (b) NI -NZ has the Poisson distribution with mean /h–pz. 27. [ZZ] (J. H. Ahrens.) On most binary computers there is an efficient way to count the number of l s in a binary word (cf. Section 7.1). Hence there is a nice way to obtain the binomial distribution (t,p) when p = f, simply by generating t random bits and counting the number of 1 s. Design an algorithm that produces the binomial distribution (t, p) for arbitrary p, using only a subroutine for the special case p = 3 as a source of random data. [Hint: Simulate a process that first looks at the most significant bits of t uniform deviates, then at the second bit of those deviates whose leading bit is not sufficient to determine whether or not their value is < p, etc.] 28. [HAL%] (R. P. Brent.) Develop a method to generate a random point on the surface of the ellipsoid defined by c ak$ = 1, where al 2 e . . 2 a, > 0. 29. [M20] (J. L. Bentley and J. D. Saxe.) Fint a simple way to generate n numbers Xl,…, X, that are uniform between 0 and 1 except for the fact that they are sorted: Xl 5 ..* 5 X,. Your algorithm should take only O(n) steps. 3.4.2. Random Sampling and Shuffling Many data processing applications call for an unbiased choice of n records at random from a file containing N records. This problem arises, for example, in quality control or other statistical calculations where sampling is needed. Usually N is very large, so that it is impossible to contain all the data in memory at once; and the individual records themselves are often very large, so that we can t even hold n records in memory. Therefore we seek an efficient procedure for selecting n records by deciding either to accept or to reject each record as it comes along, writing the accepted records onto an output file. Several methods have been devised for this problem. The most obvious approach is to select each record with probability n/N; this may sometimes be appropriate, but it gives only an average of n records in the sample. The standard deviation is n(1 - n/N), and it is possible that the sample will be either too large for the desired application, or too small to give the necessary results. A simple modification of the obvious procedure gives what we want: The (t + 1)st record should be selected with probability (n - m)/(N -t), if m items have already been selected. This is the appropriate probability, since of all the possible ways to choose n things from N such that m values occur in the first t, exactly (1)

3.4.1 NUMERICAL DISTRIBUTIONS 135 b 16. [HA&Z?] (J. H. Ahrens.) Develop an algorithm for gamma deviates of order a when 0 < a 5 1, using the rejection method with cg(t) = Y- /I(a) for 0 < t < 1, e- /r(a) for t 2 1. b 17. [M.%?4] What is the distribution function F(s) for the geometric distribution with probability p? What is the generating function G(Z)? What are the mean and standard deviation of this distribution? 18. [M.Z4] Suggest a method to compute a random integer N for which N takes the value n with probability np (1 -p) - , n 2 0. (The case of particular interest is when p is rather small.) 19. [ZZ] The negative binomial distribution (t,p) has integer values N = n with probability ( -L+,f )p (l -p) . (Unlike the ordinary binomial distribution, t need not be an integer, since this quantity is nonnegative for all n whenever t > 0.) Generalizing exercise 18, explain how to generate integers N with this distribution when t is a small positive integer. What method would you suggest if t = p = f? 20. [A&Q] Let A be the area of the shaded region in Fig. 13, and let R be the area of the enclosing rectangle. Let 1 be the area of the interior region recognized by step R2, and let E be the area between the exterior region rejected in step R3 and the outer rectangle. Determine the number of times each step of Algorithm R is performed, for each of its four variants as in (25), in terms of A, R, I, and E. 21. [HA&?9] Derive formulas for the quantities A, R, I, and E defined in exercise 20. (For 1 and especially E you may wish to use an interactive computer algebra system.) Show that c = e1j4 is the best possible constant in step R2 for tests of the form X2 5 4(1 + In c) -4cU. 22. [HM40] Can the exact Poisson distribution for large p be obtained by generating an appropriate normal deviate, converting it to an integer in some convenient way, and applying a (possibly complicated) correction a small percent of the time? 23. [HA4.23] (J. von Neumann.) Are the following two ways to generate a random quantity X equivalent (i.e., does the quantity X have the same distribution) ? Method 1: Set X + sin((r/2)U), where U is uniform. Method 8: Generate two uniform deviates, U, V, and if U2 + V2 2 1, repeat until U2 + V2 < 1. Then set X + IU2 -V l/(U + V ). 24. [HM40] (S. Ulam, J. von Neumann.) Let VO be a randomly selected real number between 0 and 1, and define the sequence (Vn) by the rule V,+l = 4V,(l -Vn). If this computation is done with perfect accuracy, the result should be a sequence with the distribution sin2 TU, where U is uniform, i.e., with distribution function F(x) = s; d+/%$-=-kj. F or i f we write V, = sin2 rU,, we find that U,+l = (2ii,) mod 1; and by the fact that almost all real numbers have a random binary expansion (see Section 3.5), this sequence U, is equidistributed. But if the computation of V, is done with only finite accuracy, the above argument breaks down because we soon are dealing with noise from the roundoff error. [Reference: von Neumann s Collected Works 5, 768-770.1 Analyze the sequence (V,) defined above when only finite accuracy is present, both empirically (for various different choices of VO) and theoretically. Does the sequence have a distribution resembling the expected distribution?

134 RANDOM NUMBERS 3.4.1 v 7. [,80] (A. J. Walker.) Suppose we have a bunch of cubes of k different colors, say nj cubes of color C, for 1 2 j 5 k, and we also have k boxes {Bl, . . . , Bk} each of which can hold exactly n cubes. Furthermore n1 + +. . + nk = kn, so the cubes will just fit in the boxes. Prove (constructively) that there is always a way to put the cubes into the boxes so that each box contains at most two different colors of cubes; in fact, there is a way to do it so that, whenever box Bj contains two colors, one of those colors is Cj. Show how to use this principle to compute the P and Y tables required in (3), given a probability distribution (~1, . . . , pk). 8. [ML51 Show that operation (3) could be changed to if U < PK then X + zK+l otherwise X t YK (i.e., using the original value of U instead of V) if this were more convenient, by suitably modifying PO, PI, . . , P&-l. 9. [HA&O] Why is the curve f(z) of Fig. 9 concave downward for 2 < 1, concave upward for x > l? b 10. [HM24] Explain how to calculate auxiliary constants Pj, Q3, Yj, Zj, Sj, Dj, Ej so that Algorithm M delivers answers with the correct distribution. b 11. [HM27] Prove that steps M7-M8 of Algorithm M generate a random variable with the appropriate tail of the normal distribution; i.e., the probability that X 5 x should be JI e-t 12dtl~ue- ./2dt, x 2 3. [Hint: Show that it is a special case of the rejection method, with g(t) = Ctept2/ for some C.] 12. [HMZI] (R. P. Brent.) Prove that the numbers oj defined in (23) satisfy the relation a;/2 -a;-,/2 < In 2 for all j 2 1. [Hint: If f(x) = ez2/ s, e–t2/2 dt, show that f(x) < f(y) for 0 I x < y.] 13. [HM.%] Given a set of n independent normal deviates, Xi, X2, . , X,, with mean 0 and variance 1, show how to find constants bj and az3, 1 5 j 2 i 5 n, so that if Yl = bl + al&l, Yz = bz + ~21×1 + azzX2, . . . , Yn = b, + ~1×1 + an2X2 +. . . + annXn, thenYi,Yz, . . . . Y, are dependent normally distributed variables, Yj has mean pj, and the Y s have a given covariance matrix (ci3). (The covariance, cij, of Yi and Yj is defined to be the average value of (Yz - pt)(Yj -pLj). In particular, Cjj is the variance of Yj, the square of its standard deviation. Not all matrices (cij) can be covariance matrices, and your construction is, of course, only supposed to work whenever a solution to the given conditions is possible.) 14. [MZ] If X is a random variable with continuous distribution F(x), and if c is a constant, what is the distribution of cX? 15. [HM.U] If Xl and X2 are independent random variables with the respective distributions Ft(x) and Fz(x), and with densities jr(x) = Fl (x), f2(x) = Fz (x), what are the distribution and density functions of the quantity X1 +X2?

3.4.1 NUMERICAL DISTRIBUTIONS 133 by making the substitution tl = sul, tl + t2 = ~~12, . . . , tl + +. . + t,-1 = XL,-1. The latter ratio is the corresponding probability that uniform deviates Vl, ** , U,-I satisfy VI 2 WI, . . . , V,-l 5 ~~-1, given that they also satisfy Ul 2 … 5 Urn-l. A more efficient but somewhat more complicated technique for binomial and Poisson deviates is sketched in exercise 22. G. For further reading. The forthcoming book Non-Uniform Random Numbers by J. H. Ahrens and U. Dieter discusses many more algorithms for the genera- tion of random variables with nonuniform distributions, together with a careful consideration of the efficiency of each technique on typical computers. From a theoretical point of view it is interesting to consider optimal methods for generating random variables with a given distribution, in the sense that the method produces the desired result from the minimum possible number of random bits. For the beginnings of a theory dealing with such questions, see D. E. Knuth and A. C. Yao, Algorithms and Complexity, ed. by J. F. Traub (New York: Academic Press, 1976), 357-428. Exercise 16 is recommended as a review of many of the techniques in this section. EXERCISES 1. [IO] If QI and ,B are real numbers with CY < ,B, how would you generate a random real number uniformly distributed between CY and p? 2. [MI61 Assuming that mU is a random integer between 0 and m -1, what is the exact probability that LkuJ = r, if 0 < r < k? Compare this with the desired probability l/k. b 3. [&I Discuss treating U as an integer and dividing by k to get a random integer between 0 and k -1, instead of multiplying as suggested in the text. Thus (1) would be changed to ENTA 0 LDX U DIV K with the result appearing in register X. Is this a good method? 4. [M20] Prove the two relations in (8). b 5. [21] Suggest an efficient method to compute a random variable with the distribu- tion px + qx2 + rx3, where p 2 0, q 2 0, r 2 0, and p + q + r = 1. 6. [HA4211 A quantity X is computed by the following method: Step 1. Generate two independent uniform deviates U, V. Step 2. If U2 + V2 2 1, return to step 1; otherwise set X +- U. What is the distribution function of X? How many times will step 1 be performed? (Give the mean and standard deviation.)