Personal computer corner

Pe ,x

Random Numbers Revisited by John Bongiovanni Telos Corporation, Santa Monica, California, USA During the course of my recent column on random number generators (RNG's), I mentioned performing reliability tests on the MTH$RANDOM linear congruential generator (LCG) found on DIGITAL computer systems (Random N umber Generation on Microcomputers, Environmental Software, March 1987). In that column I glossed over the reliability tests for two reasons: 1) Article size limitations precluded anything more than a "passing glance" at RNG testing, and 2) the tests showed no problems with the MTH$RANDOM LCG. More recently, however, I have personally found that testing a "canned" system RNG to be helpful in debugging a MonteCarlo simulation program on another computer system. Since there may be readers of this column who have encountered the same RNG problems I have, I have decided to "share the wealth" and dedicate this column to the subject of RNG testing. TEST ALL OR ONE? An absolute (but entirely true) statement is that any desired random number probability distribution can be generated provided one has a guaranteed source of independent, uniformly distributed random numbers, i.e. numbers with no mutual correlation falling with uniform probability anywhere in the interval from 0.0 to 1.0. Uniform distributions for other intervals can be obtained simply by scaling, while other distribution functions are possible by combining uniform numbers using a proper algorithm. Listing 1 is an example FORTRAN subroutine for generating normal (or Gaussian) random numbers using the system uniform RNG as a base. By and large, system routines on various computers for generating non-uniform distributions

of random numbers are based on the computer's uniform RNG, which is most commonly an LCG. If the system LCG doesn't produce truly independent, uniformly distributed random numbers, all other RNG's in the system mathematics library are in trouble. Conversely, if you test the system uniform RNG for independence and uniformity, you are basically testing all the other system RNG's as well, assuming the distribution algorithms are properly implemented on the other system routines. (Actually, LCG's have a built-in limitation in their independence that is not obvious. This fact will be brought forward later.) CHI-SQUARE MADE SIMPLE All the tests mentioned in this. column (except one) work by generating a statistic based on the RNG's output. Since it is based on the "sums-of-squaresof-the-difference-betweenobserved-and-expected", it's customarily called a Chisquare" statistic. In the independence or"ru ns" series of tests for uniform RNG's, the chi-square statistic is calculated using a special formula developed by Knuth. For the uniformity or "serial" series of tests, the calculation is more straightforward. In both series of tests, a sequence of random numbers is generated, calculations are performed, and a single number, the chi-square statistic for that test sequence, is the result. Since it is the sums of squares, chi-square is non-negative. Also, the farther your observed (RNG output) values are from your expected (true uniform RNG), the larger the numbers you're summing, so the larger your chisquare value will be. Logically, the larger the chi-square, the less like a true uniform RNG your generator is. Simple, right? However, when it comes to interpreting the chi-square statistic, many books on statistical testing get bogged down in nomenclature like "null hypothesis", "rejection criterion", and the

double-negative seeming "failed to reject". My own personal view is that with any statistical testing, the best you can do is say how certain you are that something is wrong. My recommendation is that armed with the results of the tests presented here, you should say, for example, "1 am 95% certain this generator is not random", or "1 am only 30% certain this generator is not random. I might as well use it!". In the end, it's entirely up to you to decide how certain you have to be to reject a RNG, although most literature seems to use 90% as the cutoff point. But how to get from a chisquare value to a percentage? The theory is beyond the scope of this article, but in practice, there exists an often tabulated function called, predictably, the chi-square function. It returns a chi-square value given a certainty percentage and another number called the "Degrees-of-freedom", which is most often one less that the number of items you added to get the chi-square value. Tables of chi-square are found in appendices of most books on statistics, but Abramowitz and Stegun is the most comprehensive source, as it also provides chi-square formulae for when the tables are insufficient. In using the tables, you usually go down the left side looking for your degrees-of-freedom, then go across to find your chi-square, then up to the top to get the associated certainty. The FORTRAN test procedures in this article return both chi o square and degrees-of-freedom for the purpose of looking up the certainty in a chi-square table. TESTING FOR INDEPENDENCE AND UNIFORMITY The FORTRAN subroutine RUNS (Listing 2) provides a test of the independence criterion of uniform RNG's. The algorithm (from Knuth via Law and Kelton) entails counting the number of increasing or decreasing sequences of random numbers of lengths up to six from a long

ENVIRONMENTAL SOFTWARE, 1987, Vol. 2, No. 3.

159

list of random numbers. Once these "runs-up" counts (or "runsdown", depending on the choice of increasing or decreasing) have been made, a matrixinvolved calculation is performed. The numerical result of that calculation, according to Knuth, can be treated as a chisquare statistic with six degrees of freedom. The associated certainty percentage can b e obtained from a chi-square table, using only the line on the table for six degrees of freedom, and used to determine how independent the numbers in the random number sequence are. It's a good idea to do both the runs-up and runs-down tests. Subroutine SERIAL (Listing 3) tests.for uniformity of the distribution of a random number sequence. The algorithm involves dividing a region into equal sub-regions, using random numbers to fill the sub-regions, then comparing the "expected" sub-region population of (number of random numbers) / (number of sub-regions) to the "observed" sub-region population, squaring the differences and summing over the number of sub-regions. The real trick of the serial test is that the region could be multi-dimensional, with the random numbers organized into ordered n-tuples to fill the subregions. It is best to use the serial test for several dimensions rather than just one, as will be seen. Both RUNS and SERIAL assume that the uniform RNG function is called RAN D, and that it requires a four byte integer seed value to be passed to it on each call. If either assumption is false, it is a simple task to define a function RAND that calls the real system RNG, orto modify RUNS and SERIAL to work better with your system. Both subroutines also need to be told how long a random number sequence to test. In addition, SERIAL needs to be passed several work arrays that are dimensioned in the main program which hold sub-region counts and other multidimensional information, as well as the number of dimensions to test in. ACTUAL RNG TESTS A sample main FORTRAN program using RUNS and SERIAL is in Listing 4. In addition 160

to providing an informative front end, it includes an implementation of the MTH$RANDOM LCG, for testing purposes. Figure 1 outlines the expected chisquare and degrees-of-freedom results of using listing 4 unmodified on an IBM PC using Microsoft FORTRAN 4.0, plus the corresponding percentages. As you can see, you can be at most 85% sure that the first 65536 values of MTH$RANDOM starting with a seed value of one do not behave like truly independent, uniform random numbers. Also in Figure 1 are similar results of testing the URAN uniform RNG that is an intrinsic function in HP 1000 FORTRAN, the one I had trouble with. As you can see, URAN effectively fails three of the five tests. It is important to note that URAN "passed" the one-dimensional serial test, but failed the two- and three-dimensional tests. Also, it passed the runs-down, but failed the runs-up independence test. "THE RNG IN SPACE FALLS MAINLY ON THE PLANES" The last test applies only to LCG's, which in review are of the form Z (i+1) = ( A * Z(i) + C ) mod M where all operations involve integer arithmetic and Z(i)'s are scaled to be uniform random numbers. Consider a sequence of N random numbers from an LCG: Z(1),Z(2)..... Z(N). These can be organised into N-D-1 Dtuples: (Z(1), Z(2),...,Z(D)), (Z(2), Z(3),...,Z(D+ 1))..... (Z(N-D+ 1), Z(N-D+2),...,Z(N)). If you plot these N-D-1 points in a D-space, because of the module nature of the LCG, the points will fall on at most M**(1/D) (i.e. the Dth root of M) equally spaced D-1 dimensional hyperplanes filling up the D-space. For example, if the output of an LCG based on module 2**32 were organized into ordered pairs (Z(1),Z(2)), (Z(2),Z(3))..... (Z(N-1),Z(N)) and plotted using X-Y coordinates, the points would end up on at most 2"'16 lines filling up the unit square. This would be more than enough for most purposes, even if it is a limitation in independence. However, it is important to note the "at most". A bad choice of A, C, and M could cause the

ENVIRONMENTALSOFTWARE, 1987, Vol. 2, No. 3.

number of planes to be far fewer than M**(1/D). For a good example of this fact, run the BASIC program in listing 5 on an IBM PC with CGA compatible graphics capability using version 2.1 of the IBM BASICA interpreter, if you can dig it up. The banding effect is quite striking, indicating that the number of planes is quite small for BASIC's LCG. This problem is not evident in any of the version 3.x BASICA interpreters or in the Microsoft QuickBasic compiler, indicating that they have changed the LCG's on those systems. If you have access to Pixel graphics on your computer, it may be worth your while to adapt this test to your system. Although programs do exist that can count the number of planes that show up without using graphics, your eye and brain are by far the best pattern recognistion devices for evaluating this LCG independence test. Maybe when neural networks can do pattern recognition .... FURTHER READING If you want to know more about the theory behind RNG testing, beyond the black box approach rve outlined here, you'll find the following books available. Abramowitz, M., and Stegun, I.A., Handbook of Mathematical Functions, National Bureau of Standards, 1964-1970 Knuth, D.E., Seminumerical Algorithms, Addison-Wesley, 1981 Law, A.M., and Kelton, W.D., Simulation, Modeling and Analysis, McGraw-Hill, 1982 Press, W.H., et al, Numerical Recipes, Cambridge University Press, 1986 LISTING 1 real*4 function gauss(ix) Returns gaussian distributed random numbers with zero mean and standard deviation. Other means and standard deviations are obtained simply by scaling: GAUSS2 =SDEV * GAUSS + MEAN From "Numerical Recipes", Press, et al. Actual algorithm is to

C C C C C C C C C

generate a pair of gaussian RN from a pair of uniform RN. One is returned, and the other is stored locally for the next call.

30 The matrix a(i,j) is from Knuth via Law andKeltOn c > > > > > > > > > >

IX is the seed for the Uniform RNG, RAND. integer*4 ix IFLG is used to see if a gaussian HOLD has already been stored integer iflg real*4 hold data iflg/O/ save iflg,hold

c c c c c

real*4 vl ,v2,r,fac c if (iflg.eq.O) then vl = 2.*rand(ix) - 1. v2 = 2.*rand(ix) - 1. r = v l * v l + v2*v2 if (r.ge.1 .) goto 1 fac = sqrt(-2.*alog(r)/r) gauss = v l * f a c hold = v2*fac iflg = 1 else gauss=hold iflg = 0 endif return end

LISTING 2 subroutine runs (ix,idir,n, chisq,dof) C perform the runs test on a C sequence of r a n d o m C numbers from the system C uniform RNG, RAND. IX is C the seed for RAND. IDIR C determines if the test is C runs-up (IDIR }= O, for C C increasing sequences), or runs-down (IDIR { O, for C C decreasing sequences). C The sequence tested is N numbers long. The CHISQ C returned is the calculated C runs chi-square for the test, C while the degrees-ofC freedom returned DOF is C always six (6). C C integer*4 ix,n,dof integer idir real chisq c real a(6,6),bn(6),bd (6),b(6),r(6),rb real ba,bb,rand integer*4 i,j

data a/4529.4,9044.9,13568., 18091.,22615.,27892.,9044.9, 18097.,27139.,36187.,45234., 55789.,13568.,27139.,40721., 54281.,67852.,83685.,18091., 36187.,54281.,72414.,90470., 111580.,22615.,45234. 67852.,90470.,113262., 139476.,27892.,55789., 83685.,111580.,139476., 172860./ The arrays bn(i) and bd(i) are also from Knuth via Law and Kelton data b n / 1 .,5.,11 .,19.,29.,1 ./ data bd/6.,24.,120.,720., 5040.,840./

c

1

continue chisq--rb/n dof=6

c c c c

10 c c c c c c c c c c

20

clear runs counts, calculate b's do 10 i=1,6 r(i)--O. b(i)=bn(i)/bd(i) continue generate N random numbers ba is last number, bb is current if run doesn't change direction, increment count of run, else record count and start again ba = rand(ix) j=l do 20 i=2,n bb=rand(ix) if (((ba-bb)*isign(1 ,idir)). ge.O) then if (j.gt.6)j=6 r(j)=r(j)+l. j=l else j=j+l endif ba--bb continue if (j.gt.6)j=6 r(j)=r(j)+l. Now that runs are counted, calculate chi-square Formula from Knuth via Law and Kelton

40

rb=O. do 30 i=1,6 do 40 j=1,6 rb=rb+a(i.j)*(r(i)-n*b(i))* (r(j)-n*b(j)) continue

return end LISTING 3 subroutine serial (ix,id,n,f,nf, nn,ll,chisq,dof) c perform the serial test on a c sequence of random c numbers from the system c uniform RNG, RAND. IX is c the seed for RAND. ID is c the number of dimensions c for which the test is c performed. The sequence c c tested is N numbers long. F is an interger array to c c hold the compartment c counts; it is of length NF. c NN and LL are arrays that c hold other incidental c calculations. They should be c as large as the maximum number of dimensions that c c will be tested. c The CHISQ returned is the c calculated serial chic square for the test, with the c degrees-of-freedom c returned in DOF. c integer*4 ix,n,f(1 ),nf,nn(1 ), dof integer id real chisq integer*4 i,j,k,jj,kk c c c 5 c c c c c c c

10 c c c c c c c

first, zero out f do 5 i = l , n f f(i)-O continue next, calculate largest integer id-th root of nf, jj that is the number of compartments in each dimension jj=nint(exp(alog(float(nf))/ float(id))) kk=l do 10 i=1 ,id kk=kk*jj continue if (kk.gt.nf) jj=jj-1 nn contains the dimension counts kk is the total number of compartments 11 holds base (jj) powers, to calculate offsets within f

ENVIRONMENTAL SOFTWARE, 1987, Vol. 2, No. 3.

161

20 c c c c c c c c

35 30 c c c c c c c c

40

kk=l do 20 i=l,id nn(i)=jj kk=kk*jj II(i)=kk/jj continue

ix=l n=65536 nf=4096 c c c c

Now, generate N ordered id-tuples based on RAND calculate offset into f that corresponds to compartment increment compartment count

first, the runs up and down tests call runs(ix,1 ,n,chisq,dof) write (*,*) Runs-up write(*,*) CHI-SQUARE = ,chisq,, DEG. OF FREEDOM = ,dof call runs (ix,-1 ,n,chisq,dof) write(*,*) Runs-down write(*,*) CHI-SQUARE = , chisq,, DEG. OF FREEDOM = ,dof

dp 30 il,n k=l do 35 j = l , i d k=k+ll(j)*int(nn(j)* rand(ix)) continue f(k)---f(k)+ 1 continue

next, the serial tests for 1, 2, and 3 dimensions call serial(ix,1 ,n,f,nf,nn,ll, chisq,dof) write(*,*) 1-D serial write(*,*) CHI-SQUARE = , chisq,, DEG. OF FREEDOM = ,dof call serial (ix,2,n,f,nf, nn,ll,chisq,dof) write(*,*) 2-D serial write(*,*) CHI-SQUARE = , chisq,, DEG. OF FREEDOM = ,dof call serial (ix,3,n,f,nf, nn,ll,chisq,dof) write(*,*) 3-D serial write(*,*) CHI-SQUARE = , chisq,, DEG. OF FREEDOM = ,dof

Now calculate chisq loop over all the compartments "expected" count is n/kk degrees-of-freedom is kk-1 chisq=O.O aexpect=float(n)/float(kk) do 40 i = l , k k chisq=chisq+(float(f(i)) -aexpect)**2 continue chisq=chisq/aexpect dof=kk-1

end return end

real*4 function rand(ix) c

LISTING 4 c This is an example main c program for using the c RUNS and SERIAL tests for c random number generators c c integer*4 ix,n,dof integer*4 11(3),nn(3) real*4 chisq c nf is the size of f, a work c c array for SERIAL C integer*4 nf integer*4 f(4096) c c some good numbers c

162

c c c c c c

This is the Microsoft FORTRAN 4.0 version of the MTH$RANDOM LCG. It uses the new MSF 4.0 function ISHFT. integer*4 ix,alo,xlo data alo/3533/ xlo=iand(ix,65535) ix=ishft(alo*ishft(ix,- 16)+

xlo,16)-I-alo*xlo=l rand=float(ishft(ix,-8))*5. 960464477539062e-08 return end

ENVIRONMENTALSOFTWARE, 1987, Vol. 2, No. 3.

LISTING 5 5 'Plot 10000 random points, organized to show LCG planes 10 RANDOMIZE TIMER 20 SCREEN 2:CLS 30 X=INT (640*RND) 40 FOR I=1 TO 10000 50 Y R N D = R N D : Y - I N T (200*YRND) 60 PSET (X,Y),I 70 X=INT (640*YRND) 80 N EXT I 90 END

FIGURE 1 (Figures presented as X/Y/Z are Chi-sq uare/Deg.-of-Freedom/ Percentage) MTH$RANDOM URAN(HPIOOO FORTRAN) Runs-up 9.65 / 6 / 85% 16.18 / 6 / 99.996% Runs down4.85 / 6 / 43% 44.21 / 6 / 35% Serial 1-D 4110.38/4095/57o64131.62/4095/65% Serial 2-D 4154,62/4095/75°64264.7514095/97% Serial 3-d 3987.00/4095/11%6880.62/4095/99.999+%

A N N U A L R E V I E W OF FLUID MECHANICS Milton Van Dyke, Co-Editor J. V. Wehausen, Co-Editor John L. Lumley, Associate Editor Annual Reviews Inc., Palo Alto, California Vol 18, 1986, 522 pp. This new volume of this famous series contains 16 technical papers. The articles directly related to environmental problems are: R.W. Griffiths: Gravity Currents in Rotating Systems Greg Holloway: Eddies, Waves, Circulation, and Mixing: statistical geofluid mechanics Rodney J. Sobey: Wind-Wave Prediction William D. Grant: The Continental -Shelf Bottom Boundary Layer Peter B. Rhines: Vorticity Dynamics of the Oceanic General Circulation. The volume contains, as usual, a subject index and two cumulative (Vol 1-18) indexes of authors and chapter titles.

PZ