Serial Correlation in the Generation of Pseudo-Random Numbers

• Authors: R. R. Coveyou
• Venue: Journal of the ACM (JACM)
• Year: 1960

Many practiced and proposed methods for the generation of pseudo-random numbers for use in Monte Carlo calculation can be expressed in the following way: One chooses an integer P, the base; an integer &lgr;, the multiplier, prime to P; and an integer &mgr;, the increment, less than P (&mgr; is frequently, but not always, zero). One then defines recursively a sequence x0, x1, x2, … (1) of integers by x0 = a, (2) xn+1 ≡ &lgr;xn + &mgr; (mod P), (3) 0 ≦ xn P. (4) It is clear that such a sequence is periodic with period P or less. Much work has been done on the determination of the period for various choices of P, &lgr; and &mgr; [1-7]. From this work one concludes that it is relatively easy to assure an adequately long period and that other considerations should determine the choice of these parameters.As for P, machine convenience dictates the choices P = 2q or P = 10q for a q-place binary or decimal machine, respectively. No other choice of P appears to have any practical advantage whatever. As for &lgr; and &mgr; it is the thesis of this note that they should be chosen to reduce serial correlation in the sequence (1), and we proceed to show how this may be done. We assume that &lgr; and &mgr; are such that the sequence (1) has an adequately long period. Then, clearly, one may assume, with error of order 1/P or less, that the sequence (1) is continuously and uniformly distributed on the interval (0, P). Then, if (&lpargt;Z⦔) denotes the mean value of Z, &lpargt;xn⦔ = 1/P ∫P0 x dx = P/2 (5) &lpargt;xn2⦔ = 1/P ∫P0 x2 dx = P2/3 (6) and Var (xn) = P2/3 - P2/4 = P2/12. (7) Further, one may write xn+1 = P{(&lgr;xn + &mgr;)P-1}; (8) where {···} denotes “fractional part of.” (9) Also &lpargt;xnxn+1⦔ = 1/P ∫P0 xP {(&lgr;x + &mgr;)P-1} dx. (10) An elementary but tedious integration yields (see Appendix) &lpargt;xnxn+1⦔ = P2/4 + P2 - 6&mgr;(P - &mgr;)/12&lgr;. (10′) Hence cov (xn, xn+1) = P2 - 6&mgr;(P - &mgr;)/12&lgr;, (11) and, if &rgr; (x, y) = the correlation coefficient of x and y, (12) then &rgr; (xn, xn+1) = 1 - 6 &mgr;/P (1 - &mgr;/P)/&lgr;. (13)One must be cautious of the conclusions one draws from (13) because of the approximate nature of its derivation. But one can say that a method which uses very small values of &lgr; and &mgr; is faulty. Unfortunately, such methods have been repeatedly suggested, and in one case within the author's knowledge serious difficulties were encountered in a Monte Carlo calculation because of the use of one such method. One should not lose sight of the fact that the correlation between a random number and its immediate successor is by no means the only one of importance. It is very difficult to give a precise rule, but it seems plausible that the correlation between the current random number and, at least, its next 10 or 20 successors should be controlled. Fortunately, for any such method one may write xn+p ≡ &lgr;pxn + &mgr;p (mod P), (14) with &lgr;p ≡ &lgr;p (mod P) (15) and &mgr;p ≡ &lgr;p - 1/&lgr; - 1 &mgr; (mod P), (16) and again one may estimate the correlation with the use of equation (13).