7 Random-Number Generation

Random numbers are a necessary basic ingredient in the simulation of almost all discrete systems. Most computer languages have a subroutine, object, or function that will generate a random number. Random numbers are used to generate event times and other random variables in simulation languages. In this chapter, the generation of random numbers and their subsequent testing for randomness is described.

1 Properties of Random Numbers

A sequence ofrandom numbers,$R_{1} , R_{2} , \ldots ,$, must have two important statistical properties: uniformity and independence. Each random number $R_{i}$ must be an independent sample drawn from a continuous uniform distribution between zero and 1—that is, the pdf is given by

This density function is shown in Figure 1. The expected value of each $R_{i}$ is given by

Figure 1 pdf for random numbers.

and the variance is given by

Some consequences of the uniformity and independence properties are the following:

1. If the interval [0, 1] is divided into n classes, or subintervals of equal length, the expected number of observations in each interval is $N / n ,$ where N is the total number of observations.

2. The probability of observing a value in a particular interval is independent of the previous values drawn.

In the next section we describe desirable properties of methods that produce random numbers to drive simulations.

2 Generation of Pseudo-Random Numbers

Notice that the title of this section has the word “pseudo” in it. “Pseudo” means false, so false random numbers are being generated! In this instance, “pseudo” is used to imply that the very act of generating random numbers by a known method removes the potential for true randomness. I the method is known, the set of random numbers can be repeated. Then an argument can be made that the numbers are not truly random. The goal of any generation scheme, however, is to produce a sequence of numbers between 0 and 1 that simulates, or imitates, the ideal properties of uniform distribution and independence as closely as possible.

To be sure, in the generation of pseudo-random numbers, certain problems or errors can occur. These errors, or departures from ideal randomness, are all related to the properties stated previously in Section 1. Some examples of such departures include the following:

1. The generated numbers might not be uniformly distributed.

2. The generated numbers might be discrete-valued instead of continuous-valued.

3. The mean of the generated numbers might be too high or too low.

4. The variance of the generated numbers might be too high or too low.

5. There might be dependence as measured, for instance, by autocorrelation.

Departures from uniformity and independence for a particular generation scheme often can be detected by tests such as those described in Section 4. If departures are detected, the generation scheme should be dropped in favor of an acceptable generator. Generators that pass the tests in Section 4, and tests even more stringent, have been developed; thus, there is no excuse for using a generator that has been found to be defective.

Usually, random numbers are generated by a digital computer as part of the simulation. There are numerous methods that can be used to generate the values. Before we describe some of these methods, there are several important considerations that we should mention:

1. The method should be fast. Individual computations are inexpensive, but simulation could require many millions of random numbers. The total cost can be managed by selecting a computationally efficient method of random-number generation.

2. The method should be portable to different computers—and, ideally, to different programming languages. This is desirable so that the simulation program will produce the same results wherever it is executed.

3. The method should have a sufficiently long cycle. The cycle length, or period, represents the length of the random number sequence before previous numbers begin to repeat themselves in an earlier order. Thus, if 10,000 events are to be generated, the period should be many times that long.

A special case of cycling is degenerating. A method degenerates when the same random numbers appear repeatedly. Such an occurrence is certainly unacceptable. This can happen rapidly with some methods.

4. The random numbers should be repeatable. Given the starting point (or conditions) it should be possible to generate the same set ofrandom numbers, completely independent ofthe system that is being simulated. This is helpful for debugging purposes and facilitates comparisons between systems. For the same reasons, it should be possible to easily specify differen starting points, widely separated, within the sequence.

5. Most important, the generated random numbers should closely approximate the ideal statistical properties of uniformity and independence.

Inventing techniques that seem to generate random numbers is easy; inventing techniques tha really do produce sequences that appear to be independent, uniformly distributed random numbers is incredibly difficult. There is now a vast literature and rich theory on the topic, and many hours of testing have been devoted to establishing the properties of various generators. Even when a technique is known to be theoretically sound, it is seldom easy to implement it in a way that will be fast and portable. The goal of this chapter is to make the reader aware of the central issues in random-number generation, to enhance understanding, and to show some of the techniques that are used by those working in this area.

3 Techniques for Generating Random Numbers

The linear congruential method of Section 3.1 is the most widely used technique for generating random numbers, so we describe it in detail. We also report on an extension of this method that yields sequences with a longer period. Many other methods have been proposed, and they are reviewed in Bratley, Fox, and Schrage [1996], Law [2007], and Ripley [1987].

3.1 Linear Congruential Method

The linear congruential method, initially proposed by Lehmer [1951], first produces a sequence of integers,$X_{1} , X_{2} , . . .$. between zero and $m - 1$ by following a recursive relationship:

The initial value $X_{0}$ is called the seed, a is called the multiplier, c is the increment, and m is the modulus. If $c \neq 0$ in Equation (1), then the form is called the mixed congruential method. When $c = 0 ,$, the form is known as the multiplicative congruential method. The selection of the values for $a , c , m ,$, and $X_{0}$ drastically affects the statistical properties and the cycle length. Variations of Equation (1) are quite common in the computer generation of random numbers.

Notice that random integers are being generated rather than random numbers. These random integers should appear to be uniformly distributed on the integers zero to $m .$. Random numbers $R_{i}$ between zero and 1 can then be generated by setting

An example will illustrate how this technique operates.

Example 1

Use the linear congruential method to generate a sequence of random numbers with $X_{0} = 27 , a = 17$$c = 43$, and $m = 100$ The sequence of seeds $X_{i}$ and subsequent random numbers $R_{i}$, is computed as follows:

Recall that $a = b$ mod m provided that $b - a$ is divisible by m with no remainder; otherwise, the mod operation returns the remainder of $b / m$. Thus,$X_{1} = 502$ mod $100 = 2$ because 502/100 equals 5 with a remainder of 2.

The ultimate test of the linear congruential method, as of any generation scheme, is how closely the generated numbers $R_{1} , R_{2} , \ldots$. approximate uniformity and independence. There are, however, several secondary properties that must be considered. These include maximum density and maximum period.

First, notice that the numbers generated from Equation (2) assume values only from the set $I = \{0 , 1 / m , 2 / m , \ldots , ( m - 1 ) / m \}$, because each X is an integer in the set $\{0 , 1 , 2 , \ldots , m - 1 \}$. Thus, each $R_{i}$ is discrete on I, instead of continuous on the interval [0, 1]. This approximation appears to be of little consequence if the modulus m is a very large integer. (Values such as $m = 2^{31} - 1$ and $m = 2^{48}$ are in common use in generators appearing in many simulation languages.) By maximum density is meant that the values assumed by $R_{i} , i = 1 , 2 , \ldots$, leave no large gaps on [0, 1].

Second, to help achieve maximum density, and to avoid cycling in practical applications (i.e., recurrence ofthe same sequence ofgenerated numbers), the generator should have the largest possible period. Maximal period can be achieved by the proper choice of $a , c , m .$, and $X_{0}$[Law 2007].

For m a power of $^ {2 ,}$, say $m = 2^{b}$, and $c \neq 0 .$, the longest possible period is $P = m = 2^{b}$ which is achieved whenever c is relatively prime to m (that is, the greatest common factor of c and m is 1) and $a = 1 + 4 k$, where k is an integer.

For m a power of $2 , \operatorname{s a y} m = 2^{b}$, and $c = 0$, the longest possible period is $P = m / 4 = 2^{b - 2}$ which is achieved if the seed $X_{0}$ is odd and if the multiplier a is given by $a = 3 + 8 k$ or $a = 5 + 8 k$, for some $k = 0 ,$, 1, . . . .

For m a prime number and $c = 0 .$, the longest possible period is $P = m - 1$, which is achieved whenever the multiplier a has the property that the smallest integer k such that $a^{k} - 1$ is divisible by m is $k = m - 1$

Example 2

Using the multiplicative congruential method, find the period of the generator for $a = 13 , m = 2^{6} =$64, and $X_{0} = 1 , 2 , 3$, and 4. The solution is given in Table 1. When the seed is 1 or 3, the sequence has period 16. However, a period of length 8 is achieved when the seed is 2 and a period of length 4 occurs when the seed is 4.

In Example 2,$m = 2^{6} = 64$ and $c = 0$. The maximal period is therefore $P = m / 4 = 16 .$. Notice that this period is achieved by using odd seeds,$X_{0} = 1$ and $X_{0} = 3 ;$; even seeds,$X_{0} = 2$ and $X_{0} = 4$ yield the periods 8 and 4, respectively, both less than the maximum. Notice that $a = 13$ is of the form $5 + 8 k$ with $k = 1$, as is required to achieve maximal period.

When $X_{0} = 1$, the generated sequence assumes values from the set $\{1 , 5 , 9 , 13 , \dots , 53 , 57 , 61 \}$ The gaps in the sequence of generated random numbers $R_{i}$ are quite large (i.e., the gap is $5 / 64 - 1 / 64$ or 0.0625). Such a gap gives rise to concern about the density of the generated sequence.

The generator in Example 2 is not viable for any application—its period is too short and its density is insufficient. However, the example shows the importance of properly choosing $a , c , m ,$ and $X_{0}$ Table 1 Period Determination Using Various Seeds

Speed and efficiency in using the generator on a digital computer is also a selection consideration. Speed and efficiency are aided by use of a modulus, m, which is either a power of 2 or close to a power of 2. Since most digital computers use a binary representation of numbers, the modulo, or remaindering, operation of Equation (1) can be conducted efficiently when the modulo is a power of $2 ( {\mathrm{i . e . ,}} m = 2^{b} )$. After ordinary arithmetic yields a value for $a X_{i} + c , X_{i + 1}$ is obtained by dropping the leftmost binary digits in $a X_{i} + c$ and then using only the b rightmost binary digits. The following example illustrates, by analogy, this operation using $m = 10^{b}$, because most human beings think in decimal representation.

Let $X_{0} = 63 , a = 19 , c = 0 , \mathrm{a n d} m = 10^{2} = 100$, and generate a sequence of random integers using Equation (1).

When m is a power of 10, say $m = 10^{b} .$, the modulo operation is accomplished by saving the b rightmost (decimal) digits. By analogy, the modulo operation is most efficient for binary computers when $m = 2^{b}$ for some $b > 0 ,$

Example 4

The last example in this section is in actual use. It has been extensively tested [Learmonth and Lewis, 1973; Lewis et $a l . .$, 1969]. The values for $a , c ,$ and m have been selected to ensure that the characteristics desired in a generator are most likely to be achieved.

Let $a = 7^{5} = 16 {,} 807 , c = 0$, and $m = 2^{31} - 1 = 2 , 147 , 483 , 647$(a prime number). These choices satisfy the conditions that ensure a period of $P = m - 1$(well over 2 billion). Further, specify the seed $X_{0} = 123 , 457$. The first few numbers generated are as follows:

Notice that this method divides by $m + 1$ instead of $m ;$ however, for such a large value of $m ,$ the effect is negligible.

3.2 Combined Linear Congruential Generators

As computing power has increased, the complexity ofthe systems that we are able to simulate has also increased. A random-number generator with period $2^{31} - \overset {\cdot} {1} \approx 2 \times 10^{9}$, such as the popular generator described in Example 4, is no longer adequate for all applications. Examples include the simulation of highly reliable systems, in which hundreds of thousands of elementary events must be simulated in order to observe even a single failure event, and the simulation of complex computer networks, in which thousands of users are executing hundreds of programs. An area of current research is deriving generators with substantially longer periods.

One fruitful approach is to combine two or more multiplicative congruential generators in such a way that the combined generator has good statistical properties and a longer period. The following result from L’Ecuyer [1988] suggests how this can be done:

If $W_{i , 1} , W_{i , 2} , \ldots , W_{i , k}$ are any independent, discrete-valued random variables (not necessarily identically distributed), but one of them, say $W_{i , 1}$, is uniformly distributed on the integers from 0 to $m_{1} - 2$ then

is uniformly distributed on the integers from 0 to $m_{1} - 2 .$ To see how this result can be used to form combined generators, let $X_{i , 1} , X_{i , 2} , \ldots , X_{i , k}$ be the ith output from k different multiplicative congruential generators, where thejth generator has prime modulus $m_{j}$ and the multiplier $a_{j}$ is chosen so that the period is $m_{j} - 1$. Then the jth generator is producing integers $X_{i , j}$ that are approximately uniformly distributed on the integers from 1 to $m_{j} - 1$ and $W_{i , j} = X_{i , j} - 1$ is approximately uniformly distributed on the integers from 0 to $m_{j} - 2 .$. L’Ecuyer therefore suggests combined generators of the form

with

Notice that the $^ {66} ( - 1 ) ^ {j - 1 , 5}$ coefficient implicitly performs the subtraction $X_{i , 1} - 1 ;$ for example, if $k = 2$ then $\begin{array} {r} {( - 1 ) ^ {0} ( X_{i , 1} - 1 ) - ( - 1 ) ^ {1} ( X_{i , 2} - 1 ) = \sum_{i = 1} ^ {2} ( - 1 ) ^ {j - 1} X_{i , j}} \end{array}$ The maximum possible period for such a generator is

which is achieved by the generator described in the next example.

Example 5

For 32-bit computers, L’Ecuyer [1988] suggests combining $k = 2$ generators with $a_{1} = 40 , 014 ,$$m_{1} = 2 , 147 , 483 , 563 , a_{2} = 40 , 692$, and $m_{2} = 2 , 147 , 483 , 399$. This leads to the following algorithm:

1. Select seed $X_{1 , 0}$ in the range [1, 2,147,483,562] for the first generator, and seed $X_{2 , 0}$ in the range [1, 2,147,483,398] for the second.

2. Evaluate each individual generator.

3. Set.

4. Return.

5. Set $j = j + 1$ and go to Step 2.

This combined generator has period $( m_{1} - 1 ) ( m_{2} - 1 ) / 2 \approx 2 \times 10^{18}$. Perhaps surprisingly, even such a long period might not be adequate for all applications. See L’Ecuyer [1996, 1999] and L’Ecuyer et al. [2002] for combined generators with periods as long as $2 ^ {\mathsf {I} \mathsf {\bar {9}} 1} \approx 3 \times 10^{57}$. One of these generators has been translated into Visual Basic for Applications and can be found in the spreadsheet RandomNumberTools.xls at www.bcnn.net.

3.3 Random-Number Streams

The initial seed for a linear congruential random-number generator (initial seeds, in the case of a combined linear congruential generator) is the integer value $X_{0}$ that initializes the random-number sequence. Since the sequence of integers X ,$X_{1} , \dots , X_{P} , X_{0} , X_{1} , \dots .$ produced by a generator repeats, any value in the sequence could be used to “seed” the generator.

For a linear congruential generator, a random-number stream is nothing more than a convenien way to refer to a starting seed taken from the sequence $X_{0} , X_{1} , \ldots , X_{P}$(for a combined generator, the stream is a reference to the starting seeds for all of the basic generators); typically these starting seeds are far apart in the sequence. For instance, if the streams are b values apart, then stream i could be defined by starting seed

for $i = 1 , 2 , \dots , \lfloor P / b \rfloor$. Values of $b = 100 , 000$ were common in older generators, but values as large as $b = 10^{37}$ are in use in modern combined linear congruential generators. (See, for instance, L’Ecuyer et al. [2002] for the implementation of such a generator.) Thus, a single random-number generator with k streams acts like k distinct virtual random-number generators, provided that the current value of seed for each stream is maintained. Exercise 21 illustrates one way to create streams that are widely separated in the random-number sequence.

There is the problem of comparing two or more alternative systems via simulation; there are advantages to dedicating portions of the pseudo-random number sequence to the same purpose in each of the simulated systems. For instance, in comparing the efficiency of several queueing systems, a fairer comparison will be achieved if all of the simulated systems experience exactly the same sequence of customer arrivals. Such synchronization can be achieved by assigning a specific stream to generate arrivals in each of the queueing simulations. If the starting seeds for these streams are

spaced far enough apart, then this has the same effect as having a distinct random-number generator whose only purpose is to generate customer arrivals.

4 Tests for Random Numbers

The desirable properties of random numbers—uniformity and independence—were discussed in Section 1. To check on whether these desirable properties have been achieved, a number of tests can be performed. (Fortunately, the appropriate tests have already been conducted for most commercia simulation software.) The tests can be placed in two categories, according to the properties of interest: uniformity and independence. A brief description of two types of tests is given in this section:

1. Frequency test. Uses the Kolmogorov–Smirnov or the chi-square test to compare the distribution of the set of numbers generated to a uniform distribution.

2. Autocorrelation test. Tests the correlation between numbers and compares the sample correlation to the desired correlation, zero

In testing for uniformity, the hypotheses are as follows:

The null hypothesis,$H_{0} ,$, reads that the numbers are distributed uniformly on the interval [0, 1]. Failure to reject the null hypothesis means that evidence of nonuniformity has not been detected by this test. This does not imply that further testing of the generator for uniformity is unnecessary.

In testing for independence, the hypotheses are as follows:

This null hypothesis,$H_{0} ,$, reads that the numbers are independent. Failure to reject the null hypothesis means that evidence of dependence has not been detected by this test. This does not imply that further testing of the generator for independence is unnecessary.

For each test, a level of significance α must be stated. The level α is the probability of rejecting the null hypothesis when the null hypothesis is true:

The decision maker sets the value of α for any test. Frequently, α is set to 0.01 or 0.05.

If several tests are conducted on the same set of numbers, the probability of rejecting the nul hypothesis on at least one test, by chance alone (i.e., making a Type I (α) error), increases. Say that $\alpha = 0 . 05$ and that five different tests are conducted on a sequence of numbers. The probability o rejecting the null hypothesis on at least one test, by chance alone, could be as large as 0.25.

Similarly, if one test is conducted on many sets of numbers from a generator, the probability of rejecting the null hypothesis on at least one test by chance alone (i.e., making a Type I (α) error), increases as more sets of numbers are tested. For instance, if 100 sets of numbers were subjected to the test, with $\alpha = 0 . 05$, it would be expected that five of those tests would be rejected by chance alone. If the number of rejections in 100 tests is close to 100α, then there is no compelling reason to discard the generator. The concept discussed in this and the preceding paragraph is discussed further at the conclusion of Example 8.

If one of the well-known simulation languages or random-number generators is used, it is probably unnecessary to apply the tests just mentioned or those described in Sections 4.1 and 4.2. However, random-number generators frequently are added to software that is not specifically developed for simulation, such as spreadsheet programs, symbolic/numerical calculators, and programming languages. If the generator that is at hand is not explicitly known or documented, then the tests in this chapter should be applied to many samples of numbers from the generator. Some additiona tests that are commonly used, but are not covered here, are Good’s serial test for sampling numbers [1953, 1967], the median-spectrum test [Cox and Lewis, 1966; Durbin, 1967], the runs test [Law 2007], and a variance heterogeneity test [Cox and Lewis, 1966]. Even if a set of numbers passes al the tests, there is no guarantee of randomness; it is always possible that some underlying pattern has gone undetected.

In this text, we emphasize empirical tests that are applied to actual sequences of numbers produced by a generator. Because of the extremely long period of modern pseudo-random-number generators, as described in Section 3.2, it is no longer possible to apply these tests to a significan portion of the period of such generators. The tests can be used as a check if one encounters a generator with completely unknown properties (perhaps one that is undocumented and buried deep in a software package), but they cannot be used to establish the quality of a generator throughou its period. Fortunately, there are also families of theoretical tests that evaluate the choices for m, a, and c without actually generating any numbers, the most common being the spectral test. Many o these tests assess how k-tuples of random numbers fill up a k-dimensional unit cube. These tests are beyond the scope of this text; see, for instance, Ripley [1987].

In the examples of tests that follow, the hypotheses are not restated. The hypotheses are as indicated in the foregoing paragraphs. Although few simulation analysts will need to perform these tests, every simulation user should be aware of the qualities of a good random-number generator.

4.1 Frequency Tests

A basic test that should always be performed to validate a new generator is the test of uniformity. Two different methods of testing are the Kolmogorov–Smirnov and the chi-square test. Both of these tests measure the degree of agreement between the distribution of a sample of generated random numbers and the theoretical uniform distribution. Both tests are based on the null hypothesis of no significant difference between the sample distribution and the theoretical distribution.

1. The Kolmogorov–Smirnov test. This test compares the continuous cdf $F ( x )$ of the uniform distribution with the empirical cdf $S_{N} ( x )$ of the sample of N observations. By definition,

If the sample from the random-number generator is $R_{1} , R_{2} , \ldots , R_{N}$, then the empirical cdf $S_{N} ( x )$ is defined by

As N becomes larger,$S_{N} ( x )$ should become a better approximation to $F ( x )$, provided that the null hypothesis is true.

The cdf of an empirical distribution is a step function with jumps at each observed value.

The Kolmogorov–Smirnov test is based on the largest absolute deviation between $F ( x )$ and $S_{N} ( x )$ over the range of the random variable—that ${\mathrm{i s}} ,$ it is based on the statistic

The sampling distribution of D is known; it is tabulated as a function of N in Table A.8. For testing against a uniform cdf, the test procedure follows these steps:

Step 1. Rank the data from smallest to largest. Let $R_{( i )}$ denote the ith smallest observation, so that

Step 2. Compute

Step 3. Compute $D = \operatorname* {m a x} ( D^{+} , D^{-} )$ Step 4. Locate the critical value $D_{\alpha}$ in Table A.8 for the specified significance level α and the given sample size N.

Step 5. If the sample statistic D is greater than the critical value $D_{\alpha} .$, the null hypothesis that the data are a sample from a uniform distribution is rejected. I $: D \le D_{\alpha}$, conclude that no difference has been detected between the true distribution of $\{R_{1} , R_{2} , \ldots , R_{N} \}$ and the uniform distribution.

Example 6

Suppose that the five numbers 0.44, 0.81, 0.14, 0.05, 0.93 were generated, and it is desired to perform a test for uniformity by using the Kolmogorov–Smirnov test with the level of significance $\alpha = 0 . 05$ First, the numbers must be sorted from smallest to largest. The calculations can be facilitated by use of Table 2. The top row lists the numbers from smallest $( R_{( 1 )} )$ to largest $( R_{( 5 )} )$. The computations for $D^{+}$, namely $i / N - R_{( i )}$, and for $D^{-}$, namely $R_{( i )} - ( i - 1 ) / N$, are easily accomplished by using Table 2. The statistics are computed as $D^{+} = 0 . 26$ and $D^{-} = 0 . 21$. Therefore,$D = \operatorname* {m a x} \{0 . 26 , 0 . 21 \} = 0 . 26$ Table 2 Calculations for Kolmogorov–Smirnov Test

The critical value of $D ,$ obtained from Table A.8 in the Appendix for $\alpha = 0 . 05$ and $N = 5 ,$, is 0.565. Since the computed value, 0.26, is less than the tabulated critical value, 0.565, the hypothesis tha the distribution of the generated numbers is the uniform distribution is not rejected.

The calculations in Table 2 are illustrated in Figure 2, where the empirical cdf $S_{N} ( x )$ is compared to the uniform cdf $F ( x )$. It can be seen that $D^{+}$ is the largest deviation of $S_{N} ( x )$ above $F ( x )$), and that $D^{-}$ is the largest deviation of $S_{N} ( x )$ below $F ( x )$. For example, at $R_{( 3 )}$, the value of $D^{+}$ is given by $3 / 5 - R_{( 3 )} = 0 . 60 - 0 . 44 = 0 . 16 {\mathrm{.}}$, and that of $D^{-}$ is given by $R_{( 3 )} - 2 / 5 = 0 . 44 - 0 . 40 = 0 . 04$ Although the test statistic D is defined by Equation (3) as the maximum deviation over all x, it can be seen from Figure 2 that the maximum deviation will always occur at one of the jump points $R_{( 1 )} , R_{( 2 )} , \ldots , R_{( N )} ;$; thus, the deviation at other values of x need not be considered.

2. The chi-square test. The chi-square test uses the sample statistic

where $O_{i}$ is the observed number in the ith class,$E_{i}$ is the expected number in the ith class, and n is the number of classes. For the uniform distribution $E_{i} ,$, the expected number in each class is given by

for equally spaced classes, where N is the total number of observations. It can be shown that the sampling distribution of $\chi_{0} ^ {2}$ is approximately the chi-square distribution with n 1 degrees offreedom.

Example 7

Use the chi-square test with $\alpha = 0 . 05$ to test for whether the data shown next are uniformly distributed. Table 3 contains the essential computations. The test uses $n = 10$ intervals of equal length, namely $[ 0 , 0 . 1 ) , [ 0 . 1 , 0 . 2 ) , \ldots , [ 0 . 9 , 1 . 0 )$. The value of $\chi_{0} ^ {2}$ is 3.4. This is compared with the critical value $\chi_{0 . 05 , 9} ^ {2} = 16 . 9$ from Table A.6. Since $\chi_{0} ^ {2}$ is much smaller than the tabulated value of $\chi_{0 . 05 , 9} ^ {2}$, the null hypothesis of a uniform distribution is not rejected.

Figure 2 Comparison of F(x) and $S_{N} ( x )$| 0.34 | 0.90 | 0.25 | 0.89 | 0.87 | 0.44 | 0.12 | 0.21 | 0.46 | 0.67 |

Different authors have offered considerations concerning the application of the chi-square test. In the application to a data set the size of that in Example 7, the considerations do not apply—that is, if 100 values are in the sample and from 5 to 10 intervals of equal length are used, the test will be acceptable. In general, it is recommended that n and N be chosen so that each $E_{i} \geq 5$ Both the Kolmogorov–Smirnov test and the chi-square test are acceptable for testing the uniformity of a sample of data, provided that the sample size is large. However, the Kolmogorov–Smirnov test is the more powerful of the two and is recommended. Furthermore, the Kolmogorov–Smirnov test can be applied to small sample sizes, whereas the chi-square is valid only for large samples, say $N \geq 50$ Table 3 Computations for Chi-Square Test

Imagine a set of 100 numbers that are being tested for independence, one where the first 10 values are in the range 0.01–0.10, the second 10 values are in the range 0.11–0.20, and so on. This set of numbers would pass the frequency tests with ease, but the ordering of the numbers produced by the generator would not be random. The test in the next section of this chapter is concerned with the independence of random numbers that are generated.

4.2 Tests for Autocorrelation

The tests for autocorrelation are concerned with the dependence between numbers in a sequence. As an example, consider the following sequence of numbers, read from left to right:

From a visual inspection, these numbers appear random, and they would probably pass all the tests presented to this point. However, an examination of the 5th, 10th, 15th (every five numbers beginning with the fifth), and so on, indicates a very large number in that position. Now, 30 numbers is a rather small sample size on which to reject a random number generator, but the notion is that numbers in the sequence might be related. In this particular section, a method for discovering whether such a relationship exists is described. The relationship would not have to be all high numbers. It is possible to have all low numbers in the locations being examined, or the numbers could alternate from very high to very low.

The test to be described shortly requires the computation of the autocorrelation between every  numbers ( is also known as the $l a g ) ,$, starting with the ith number. Thus, the autocorrelation $\rho_{i \ell}$ between the following numbers would be of interest:$R_{i} , R_{i + \ell} , R_{i + 2 \ell} , \ldots , R_{i + ( M + 1 ) \ell}$. The value M is the largest integer such that $i + ( M + 1 ) \ell \leq N$, where N is the total number of values in the sequence. (Thus, a subsequence of length $M + 2$ is being tested.)

Figure 3 Failure to reject hypothesis.

A nonzero autocorrelation implies a lack of independence, so the following two-tailed test is appropriate:

For large values of M, the distribution of the estimator of $\rho_{i \ell}$, denoted $\widehat {\rho} _ {i \ell}$, is approximately normal if the values $R_{i} , R_{i + \ell} , R_{i + 2 \ell} , \ldots , R_{i + ( M + 1 ) \ell}$ are uncorrelated. Then the test statistic can be formed as follows:

which is distributed normally with a mean of zero and a variance of 1, under the assumption of independence, for large M.

The formula for $\widehat {\rho} _ {i \ell}$, in a slightly different form, and the standard deviation of the estimator $\sigma _ {\hat {\rho} _ {i \ell}}$ are given by Schmidt and Taylor [1970] as follows:

and

After computing $Z_{0} ,$ do not reject the null hypothesis of independence $\mathrm{i f} - z_{\alpha / 2} \le Z_{0} \le z_{\alpha / 2}$ where α is the level of significance and $z_{\alpha / 2}$ is obtained from Table A.3. Figure 3 illustrates this test.

If $\rho_{i \ell} > 0 .$, the subsequence is said to exhibit positive autocorrelation. In this case, successive values at lag  have a higher probability than expected of being close in value $( \mathrm{i . e .}$, high random numbers in the subsequence followed by high, and low followed by low). On the other hand, if $\rho_{i \ell} < 0$ the subsequence is exhibiting negative autocorrelation, which means that low random numbers tend

to be followed by high ones, and vice versa. The desired property, independence (which implies zero autocorrelation), means that there is no discernible relationship of the nature discussed here between successive random numbers at lag .

Example 8

Test for whether the 3rd, 8th, 13th, and so on, numbers in the sequence at the beginning of this section are autocorrelated using $\alpha = 0 . 05$. Here, i 3 (beginning with the third number),$\ell = 5$(every five numbers),$N = 30$(30 numbers in the sequence), and $M = 4$(largest integer such tha $3 + ( M + 1 ) 5 \le 30 )$. Then,

and

Then, the test statistic assumes the value

Now, the critical value from Table A.3 is

Therefore, the hypothesis of independence cannot be rejected on the basis of this test.

It can be observed that this test is not very sensitive for small values of M, particularly when the numbers being tested are on the low side. Imagine what would happen if each of the entries in the foregoing computation of $\widehat {\rho} _ {i \ell}$ were equal to zero. Then $\widehat {\rho} _ {i \ell}$ would be equal to 0.25 and the calculated Z would have the value of 1.95, not quite enough to reject the hypothesis of independence.

There are many sequences that can be formed in a set of data, given a large value of N. For example, beginning with the first number in the sequence, possibilities include (a) the sequence of all numbers, (b) the sequence formed from the first, third, fifth, . . . , numbers, (c) the sequence formed from the first, fourth, . . . , numbers, and so on. If $\alpha = 0 . 05$, there is a probability of 0.05 of rejecting a true hypothesis. If 10 independent sequences are examined, the probability of finding no significan autocorrelation, by chance alone, is $( 0 . 95 ) ^ {10}$ or 0.60. Thus, 40% of the time significant autocorrelation would be detected when it does not exist. If α is 0.10 and 10 tests are conducted, there is a 65% chance of finding autocorrelation by chance alone. In conclusion, in “fishing” for autocorrelation by performing numerous tests, autocorrelation might eventually be detected, perhaps by chance alone, even when there is no autocorrelation present.

5 Summary

This chapter described the generation of random numbers and the subsequent testing of the generated numbers for uniformity and independence. Random numbers are used to generate random variates.

Of the many types of random-number generators available, ones based on the linear congruentia method are the most widely used, but they are being replaced by combined linear congruentia generators. Of the many types of statistical tests that are used in testing random-number generators, two different types are described: one testing for uniformity and one testing for independence.

The simulation analyst might never work directly with a random-number generator or with the testing of random numbers from a generator. Most computers and simulation languages have methods that generate a random number, or streams of random numbers, for the asking. But even generators that have been used for years, some of which are still in use, have been found to be inadequate. So this chapter calls the simulation analyst’s attention to such possibilities, with a warning to investigate and confirm that the generator has been tested thoroughly. Some researchers have attained sophisticated expertise in developing methods for generating and testing random numbers and the subsequent application of these methods. This chapter provides only a basic introduction to the subject matter; more depth and breadth are required for the reader to become a specialist in the area. An importan foundational reference is Knuth [1998]; see also the reviews in Bratley, Fox, and Schrage [1996], Law [2007], L’Ecuyer [1998], and Ripley [1987].

One final caution is due. Even if generated numbers pass all the tests (those covered in this chapter and those mentioned in the chapter), some underlying pattern might have gone undetected without the generator’s having been rejected as faulty. However, the generators available in widely used simulation languages have been extensively tested and validated.

REFERENCES

BRATLEY, P., B. L. FOX, AND L. E. SCHRAGE [1996], A Guide to Simulation, 2d ed., Springer-Verlag, New York.

COX, D. R., AND P. A. W. LEWIS [1966], The Statistical Analysis ofSeries ofEvents, Methuen, London.

DURBIN, J. [1967], “Tests of Serial Independence Based on the Cumulated Periodogram,” Bulletin of the International Institute ofStatistics, vol. 42, pp. 1039–1049.

GOOD, I. J. [1953], “The Serial Test for Sampling Numbers and Other Tests of Randomness,” Proceedings of the Cambridge Philosophical Society, Vol. 49, pp. 276–284.

GOOD, I. J. [1967], “The Generalized Serial Test and the Binary Expansion of 4,” Journal ofthe Royal Statistical Society, Ser. A, Vol. 30, No. 1, pp. 102–107.

KNUTH, D. W. [1998], TheArt ofComputer Programming:Vol. 2, Semi-numericalAlgorithms, 2d ed.,Addison– Wesley, Reading, MA.

LAW, A. M. [2007], Simulation Modeling and Analysis, 4th ed., McGraw-Hill, New York.

LEARMONTH, G. P., AND P. A. W. LEWIS [1973], “Statistical Tests of Some Widely Used and Recently Proposed Uniform Random Number Generators,” Proceedings ofthe Conference on Computer Science and Statistics: Seventh Annual Symposium on the Interface, Western Publishing, North Hollywood, CA, pp. 163–171.

L’ECUYER, P. [1988], “Efficient and Portable Combined Random Number Generators,” Communications of the ACM, Vol. 31, pp. 742–749, 774.

L’ECUYER, P. [1996], “Combined Multiple Recursive Random Number Generators,” Operations Research, Vol. 44, pp. 816–822.

L’ECUYER, P. [1998], “Random Number Generation,” Chapter 4 in Handbook of Simulation, J. Banks, ed., pp. 93–137. Wiley, New York.

L’ECUYER, P. [1999], “Good Parameters and Implementations for Combined Multiple Recursive Random Number Generators,” Operations Research, Vol. 47, pp. 159–164.

L’ECUYER, P., R. SIMARD, E. J. CHEN,AND W. D. KELTON [2002], “An Object-Oriented Random-Number Package with Many Long Streams and Substreams,” Operations Research, Vol. 50, pp. 1073–1075.

LEHMER, D. H. [1951], “Mathematical Methods in large-Scale Computing Units,” Proceedings ofthe Second Symposium on Large-Scale Digital Calculating Machinery, Harvard University Press, Cambridge, MA, pp. 141–146.

LEWIS, P. A. W., A. S. GOODMAN, AND J. M. MILLER [1969], “A Pseudo-Random Number Generator for the System/360,” IBM Systems Journal, Vol. 8, pp. 136–145.

RIPLEY, B. D. [1987], Stochastic Simulation, Wiley, New York.

SCHMIDT, J. W., and R. E. TAYLOR [1970], Simulation and Analysis of Industrial Systems, Irwin, Home wood, IL.

EXERCISES

1. Describe a procedure to physically generate random numbers on the interval [0, 1] with 2-digit accuracy. [Hint: Consider drawing something out of a hat.]

2. List applications, other than systems simulation, for pseudo-random numbers—for example, video gambling games.

3. How could random numbers that are uniform on the interval [0, 1] be transformed into random numbers that are uniform on the interval [ 11, 17]?

4. Use the linear congruential method to generate a sequence of three two-digit random integers and corresponding random numbers. Let X 27, a 8, c 47, and m 100.

5. Do we encounter a problem in the previous exercise if $X_{0} = 0 ?$6. Use the multiplicative congruential method to generate a sequence of four three-digit random integers and corresponding random numbers. Let $X_{0} = 117$, a 43, and m 1000.

6. The sequence ofnumbers 0.54, 0.73, 0.98, 0.11, and 0.68 has been generated. Use the Kolmogorov– Smirnov test with $\alpha = 0 . 05$ to learn whether the hypothesis that the numbers are uniformly distributed on the interval [0, 1] can be rejected.

7. Reverse the 100 two-digit random numbers in Example 7 to get a new set of random numbers. Thus, the first random number in the new set will be 0.43. Use the chi-square test, with $\alpha = 0 . 05$ to learn whether the hypothesis that the numbers are uniformly distributed on the interval [0, 1] can be rejected.

8. Figure out whether these linear congruential generators can achieve a maximum period; also, state restrictions on $X_{0}$ to obtain this period:

(a) the mixed congruential method with

(b) the multiplicative congruential generator with

(c) the mixed congruential generator with

(d) the multiplicative congruential generator with

10. Use the mixed congruential method to generate a sequence of three two-digit random integers and corresponding random numbers with $X_{0} = 37 , a = 7 , c = 29$, and $m = 100$11. Use the mixed congruential method to generate a sequence of three two-digit random integers between 0 and 24 and corresponding random numbers with $X_{0} = 13 , a = 9$, and $c = 35$12. Write a computer program that will generate four-digit random numbers, using the multiplicative congruential method. Allow the user to input values of $X_{0} , a , c .$, and $m .$13. If $X_{0} = 3579$ in Exercise $9 ( \mathrm{c} ) .$, generate the first random number in the sequence. Compute the random number to four-place accuracy.

11. Investigate the random-number generator in a spreadsheet program on a computer to which you have access. In many spreadsheets, random numbers are generated by a function called RAND or @RAND.

(a) Check the user’s manual to see whether it describes how the random numbers are generated.

(b) Write macros to conduct each of the tests described in this chapter. Generate 100 sets o random numbers, each set containing 1000 random numbers. Perform each test on each set of random numbers. Draw conclusions.

15. Consider the multiplicative congruential generator under the following circumstances:

(a) $X_{0} = 7 , a = 11 , m = 16$(b) $X_{0} = 8 , a = 11 , m = 16$(c) $X_{0} = 7 , a = 7 , m = 16$

Generate enough values in each case to complete a cycle. What inferences can be drawn? Is maximum period achieved?

16. L’Ecuyer [1988] provides a generator that combines three multiplicative generators, with $a_{1} =$157,$m_{1} = 32 , 363 , a_{2} = 146 , m_{2} = 31 , 727 , a_{3} = 142$, and $m_{3} = 31 , 657$. The period of this generator is approximately $8 \times 10^{12}$. Generate 5 random numbers with the combined generator, using the initial seeds $X_{1 , 0} = 100 , X_{2 , 0} = 300$ and $X_{3 , 0} = 500$, for the individual generators.

17. Apply the tests described in this chapter to the generator given in the previous exercise.

18. Use the principles described in this chapter to develop your own linear congruential randomnumber generator.

19. Use the principles described in this chapter to develop your own combined linear congruentia random-number generator.

20. Test the following sequence of numbers for uniformity and independence, using procedures you learned in this chapter: 0.594, 0.928, 0.515, 0.055, 0.507, 0.351, 0.262, 0.797, 0.788, 0.442, 0.097, 0.798, 0.227, 0.127, 0.474, 0.825, 0.007, 0.182, 0.929, 0.852.

21. In some applications, it is useful to be able to quickly skip ahead in a pseudo-random number sequence without actually generating all of the intermediate values.

(a) For a linear congruential generator with $c = 0 .$, show that $X_{i + n} = ( a^{n} X_{i} )$ mod $m .$.

(b) Next, show that $( a^{n} X_{i} )$ mod $m = ( a^{n}$ mod $m ) X_{i}$ mod m (this result is useful because $a^{n}$ mod m can be precomputed, making it easy to skip ahead n random numbers from any point in the sequence).

(c) In Example 3, use this result to compute $X_{5} ,$, starting with $X_{0} = 63$. Check your answer by computing $X_{5}$ in the usual way.

22. Producing good pseudo-random-number generators is much more difficult than it seems. For instance, mathematical genius John Von Neumann proposed the middle-square method. To generate, say, six-digit pseudo-random numbers by the middle-square method, begin with a six-digit integer as the seed. Square the number and extract the middle six digits as the next number in the sequence. Continue the process by repeatedly squaring the current number and extracting the middle six digits as the next number. To obtain numbers between 0 and 1, a decimal point can be placed in front of the six digits. Try a few examples to show how this method can fail spectacularly. (Von Neumann was aware of these problems.)

List of Materials Available on www.bcnn.net

RandomNumberTools.xls