Newsgroups: sci.stat.math,sci.math,sci.math.num-analysis Path: bga.com!news.sprintlink.net!news.onramp.net!convex!cs.utexas.edu!swrinde!ihnp4.ucsd.edu!library.ucla.edu!csulb.edu!csus.edu!netcom.com!deleyd From: deleyd@netcom.com Subject: Random Numbers - CHIKSN.FOR Message-ID: <deleydCsIBAp.KuF@netcom.com> Organization: NETCOM On-line Communication Services (408 261-4700 guest) Date: Wed, 6 Jul 1994 06:56:49 GMT Lines: 526 Xref: bga.com sci.stat.math:1316 sci.math:15354 sci.math.num-analysis:3355 {Approx. 520 lines} C CHIKSN.FOR C C This is the program which impliments the chi-square test used to test C random number generators. Presented here if you would wish to play C with it yourself, maybe do some testing of your own. See the paper C "Computer Generated Random Numbers" sections 4 and 6 for an explanation C of what this program does. C C This is not polished code you put on a shelf and admire, this is code C you dig your hands into and work with to make it do what you want it to. C C The main routine is the one to tinker with. This code is meant to be C modified to suit your needs. The goal is to fill up the bins with C balls. Here's a brief outline: C C 1. ASK USER INPUT: C a. Number of dimensions? NDIM C b. Number of bins per dimension? NBINSPD C c. Total number of balls to throw at bins? NBALLS C d. Number of tests to run? NCHITESTS C e. Random number generator to use (if more than one defined) C f. SEED value to initialize generator with C C 2. CREATE ARRAY BINS and ZERO ARRAY C C 3. THROW THE BALLS AT THE BINS and CALCULATE PROBABILITY: C C LOOP (do CHITEST=1 to NCHITESTS) C zero BIN array C LOOP (do J=1 to NBALLS) C get random numbers r1,r2,...,rn) C increment BIN(r1,r2,...,rn) C ENDLOOP C CALL CHSONE to calculate chi-square probability C ENDLOOP C CALL KSONE to calculate Kolmogorov-Smirnov probability C C The main routine here is where you put a call to your random number C generator, or for speed you can attempt a direct implimentation of your C random number generator to save the overhead of a call. (It can make a C difference when you call the random number generator 100 million times C in one test). C C The main routine defines a large one-dimensional array called BINS, the C maximum size of which would depend on your account quotas and machine C specific limitations. The array BINS keeps track of how many balls have C fallen into each bin. The size of array BINS determines the maximum C number of bins a user may select for a test. (10,000,000 bins is a C typical number you may want to use if possible.) C C So steps are: C C 1. Check definition of one-dimensional array NBINS C that it's not too large for your account quotas C or system limitations. C C 2. Place your random number generator to be tested where it C bluntly says "PLACE YOUR RANDOM NUMBER GENERATOR HERE". C The output is an integer IRANDOM between 0 and NBINS-1 C (NBINS is the number of bins chosen by the user). C C Currently the program is set up to use subroutine RAN1, a portable C random number generator from the book "NUMERICAL RECIPES: The Art C of Scientific Computing". I've had trouble on our UNIX system C not making array R(97) static even though the code says to. C Compiling with the -static qualifier works. C C The random number generator being tested is used to "randomly" select a C bin for the ball to fall in, and the counter for that bin is C incremented. Note for a multi-dimensional test we calculate the C appropriate index into the linear array BINS by hand. After all the C balls are thrown we call the subroutines to do the heavy math. C C All the subroutines should be fairly standard FORTRAN-77 modified C versions of routines from the book "NUMERICAL RECIPES: The Art of C Scientific Computing" by William H. Press, Brian P. Flannery, Saul A. C Teukolsky, and William T. Vetterling, and you should look there for C further reference as to what the routines are doing. (Note: the book C comes in several programming language forms including C, PASCAL, BASIC, C as well as FORTRAN, so you can take your pick and rewrite this code in C any language you please.) C C Note: C CHIKSN.FOR is currently set up to be run by a process with a very C large page file quota (pgflquo). If you get a 'exceed quota' error C attempting to run this then all you need to do is change the line C which reads: C C INTEGER*2 BINS(20 000 000) !The bins. C C to something smaller like: C C INTEGER*2 BINS(1 000 000) !The bins. C C C To compile: C $ FORTRAN CHIKSN C $ LINK CHIKSN C or C % f77 chiksn.f !(Some UNIX F77 compilers require -save option) C C Sample run: C Test MTH$RANDOM in 3-D with 10 bins per dimension and 10 balls per bin: C C $ RUN CHIKSN C Input number of dimensions NDIM: 3 C Input number of bins per dimension NBINSPD: 10 C Total number of bins = NBINSPD**NDIM = 1000 C Minimum number of balls = 5*NBINS = 5000 C Input total number of balls NBALLS: 10000 C Input number of Chi-Square tests NCHITESTS (min=2) : 2 C Choose random number generator to test C (1) MTH$RANDOM, C (2) RANDU, C (3) ANSI C, C (4) Microsoft C C (5) Turbo Pascal C (6) DES C : 1 C Input starting SEED value: 1 C C BALLS= 10000 CHISQ= 993.0002441 PROB= 0.4524292 C BALLS= 10000 CHISQ= 974.0001831 PROB= 0.2915459 C KS D= 0.5475708 PROB= 0.5863269 C C----------------------------------------------------------------------- PROGRAM CHIKSN C Perform a CHI-SQUARE test on a sequence of sets of N random numbers C NDIM = number of dimensions C NBINSPD = number of bins per dimension C NBINS = total number of bins. NBINS = NBINSPD**NDIM C NBALLS = total number of balls. Should be at least 5*(NBINS**NDIM) C NCHITESTS = Number of chi-square tests to do. Must be 2 or more. C EBINS = Expected value for each bin. EBINS = NBALLS/NBINS C SEED = Initial seed value for random number generator C C Note 1: The maximum size of array NBINS may be determined by the users C page file quota (pgflquo in AUTHORIZE). Also, it is recommended C the user have a very large working set quota (wsquo,wsextent) C to reduce page faulting. This can greatly improve speed. C We use INTEGER*2 array here to save space. C C Note 2: The maximum number of Chi-Square tests that can be saved is C arbitrary (array SAVEPROB). The user may choose any value. C C Note 3: The user may choose any starting seed value. Some restrictions C may apply depending upon the particular random number generator C being used. For example, RANDU should always be started with C an odd value of SEED. MTH$RANDOM may be started with any value C of SEED. C C Note 4: The MTH$RANDOM generator is used by the VAX FORTRAN intrinsic C function RAN and the VAX BASIC function RND. It is defined as: C C SEED = 69069*SEED + 1 mod 2**32 C X = SEED/2**32 C C Note 5: The RANDU generator is obsolite due to very strong correlation C in 3d space. ( Prove to yourself using 65539 = 2**16 + 3 that C SEED[i+2] = 6*SEED[i+1] - 9*SEED[i] ). It is defined as: C C SEED = 65539*SEED mod 2**31 C X = SEED/2**31 C C The RANDU generator should be started with an odd value of SEED. C C Note 6: The C standard library function rand() is defined as: C C SEED = 1103515245*SEED + 12345 mod 2**32 C IX = SEED mod 2**31 C C This standard random number generator is defined in the book: C The C Programming Language C Brian W. Kernighan and Dennis M. Ritchie C Prentice Hall, 1978 C C The same generator is defined in the ANSI C version by the same C authors above, and the same generator is used in VAX C. C C Note 7: The Microsoft C version 4.0 library function rand() impliments C the following: C C SEED = 214013*SEED + 2531011 mod 2**32 C IX = bits 16-31 of SEED C C Note 8: The Turbo Pascal version 6.0 function impliments the following: C C SEED = 134775813*SEED + 1 mod 2**32 C IX = bits 16-32 of SEED C IMPLICIT NONE INTEGER NDIM !Number of dimensions INTEGER NBINS !Number of bins INTEGER NBINSPD !Number of bins per dimension INTEGER NBALLS !Number of random numbers per !chi-square test. INTEGER NCHITESTS !Number of chi-square tests to do C INTEGER*2 BINS(20 000 000) !The bins. (see note 1) INTEGER*2 BINS(200 000) !Less bins. (see note 1) REAL EBINS !Expected number of balls per bin REAL SAVEPROB(100) !Array to save results of !chi-square tests (see note 2) INTEGER*4 SEED(2) !Only SEED(1) result is ever used. INTEGER*2 W(4) !Seeds for RANDU EQUIVALENCE(SEED,W) !for RANDU COMMON / SEEDSTORE / SEED INTEGER I,J,K,MRANDO,NBYTES,CLEAR,CHITEST,INDEX,IRANDOM,NCLEAR REAL*4 FRANDOM, RRANDOM CHARACTER*8 TIMEBUF EQUIVALENCE (IRANDOM,FRANDOM) REAL FOR$IRAN !The RANDU random number generator REAL RAND !UNIX rand() INTEGER*4 xrand !BSD random() REAL RAN1 !test generator supplied REAL D INTEGER JISHFT,IRANDOM2,COUNT REAL RANDES !DES FUNCTION (not supplied) INTEGER KEY(2) REAL CHSQ,PROB !Chi-square value, !chi-square probability REAL*4 FNBINSPD !float(NBINSPD) REAL*4 TWO31F REAL*4 TWO16F REAL*4 TWO15F TWO31F = 2.0**31.0 TWO16F = 2.0**16.0 TWO15F = 2.0**15.0 C*DES KEY(1) = 12345 !Choose any number you want C*DES KEY(2) = 678901 !to initialize DES with C*DES CALL DES_INIT(KEY) !DES code not included. 104 FORMAT(' Input number of dimensions NDIM: ',$) 100 FORMAT(' Input number of bins per dimension NBINSPD: ',$) 105 FORMAT(' Total number of bins = NBINSPD**NDIM = ',I) 106 FORMAT(' Minimum number of balls = 5*NBINS = ',I) 101 FORMAT(' Input total number of balls NBALLS: ',$) 103 FORMAT(' Input number of Chi-Square tests NCHITESTS (min=2) : ',$) 102 FORMAT(' Choose random number generator to test'/, 1 ' /*(1)*/ xrand(),'/ 1 ' /*(2)*/ UNIX rand(),'/ 1 ' /*(3)*/ MTH$RANDOM,'/ 2 ' /*(4)*/ RANDU,'/ 3 ' /*(5)*/ ANSI C,'/ 4 ' /*(6)*/ Microsoft C'/ 5 ' /*(7)*/ Turbo Pascal'/ 7 ' /*(8)*/ DES'/ 8 ' (9) another random number generator (choose this one)'/ 6 ' : ',$) 107 FORMAT(' Input starting SEED value: ',$) 200 FORMAT(BN,I) C ***GET USER INPUT*** 10 WRITE(6,104) !Input number of dimensions READ(5,200) NDIM WRITE(6,100) !Input number of bins per dimension READ(5,200) NBINSPD FNBINSPD = FLOAT(NBINSPD) NBINS = NBINSPD**NDIM !Calculate total number of bins WRITE(6,105) NBINS !Total number of bins is... WRITE(6,106) 5*NBINS !Minimum number of balls is... WRITE(6,101) !Input total number of balls READ(5,200) NBALLS WRITE(6,103) !Input number of chi-square tests to do READ(5,200) NCHITESTS WRITE(6,102) !Choose random number generator to test READ(5,200) MRANDO WRITE(6,107) !Starting SEED value READ(5,200) SEED(1) SEED(2) = 1 !Used only if random number generator !uses bigger than 32 bits C INITIALIZE GENERATOR IF NEEDED C*XRAND CALL sxrand(SEED(1)) !Initialize xrand() CALL RAN1(-SEED(1)) !Initialize RAN1 generator C Calculate expected average number of balls for each bin EBINS = FLOAT(NBALLS)/FLOAT(NBINS) C CALL TIME(TIMEBUF) C WRITE(6,201) TIMEBUF C201 FORMAT(1X,A8) DO CHITEST=1,NCHITESTS C *** ZERO BIN ARRAY *** DO I=1,NBINS BINS(I) = 0 ENDDO C*VMS !Quickly set BINS(k) = 0, k=1,...NBINS C*VMS !Does the equivalent of above C*VMS !but a lot faster. C*VMS K = 1 C*VMS NBYTES = NBINS*2 !total number of bytes to zero C*VMS DO WHILE (NBYTES .GT. 0) C*VMS IF (NBYTES .LE. 65534) THEN !maximum number of bytes we can clear C*VMS NCLEAR = NBYTES !in one call to LIB$MOVC5 is 65535 C*VMS ELSE C*VMS NCLEAR = 65534 !max that LIB$MOVC3 can do in one call C*VMS ENDIF !(make nclear an even number so we can divide evenly by 2) C*VMS CALL LIB$MOVC5(0,0,0,NCLEAR,BINS(K)) !Clear a block of memory C*VMS NBYTES = NBYTES - NCLEAR !Number of bytes still left to clear C*VMS K = K + NCLEAR/2 !Number of bytes cleared so far + 1 C*VMS ENDDO C Main Loop DO J=1,NBALLS INDEX = 1 DO I=0,NDIM-1 C ***PLACE YOUR RANDOM NUMBER GENERATOR HERE*** C Set IRANDOM using whatever random number generator you choose C IRANDOM = integer between 0 and NBINS-1 c IF (MRANDO .EQ. 1) THEN c IRANDOM = INT( ( float( xrand() ) /TWO31F ) *FNBINSPD) c ELSEIF (MRANDO .EQ. 2) THEN c IRANDOM = INT( RAND(SEED(1)) *FNBINSPD) !UNIX rand() c ELSEIF (MRANDO .EQ. 3) THEN c IRANDOM = INT( RAN(SEED(1)) *FNBINSPD) !VMS mth$random c ELSEIF (MRANDO .EQ. 4) THEN c IRANDOM = INT( FOR$IRAN(W(2),W(1)) *FNBINSPD) !Infamous randu c ELSEIF (MRANDO .EQ. 5) THEN c CALL LIB$EMUL(1103515245,SEED,12345,SEED) !ANSI C c IRANDOM = SEED(1) .AND. '7FFFFFFF'X c IRANDOM = INT( FLOAT(IRANDOM)/(TWO31F) *FNBINSPD) c ELSEIF (MRANDO .EQ. 6) THEN c CALL LIB$EMUL(214013,SEED,2531011,SEED) !Microsoft C 4.0 c IRANDOM = W(2) .AND. '7FFF'X c IRANDOM = INT( FLOAT(IRANDOM)/(TWO15F) *FNBINSPD) c ELSEIF (MRANDO .EQ. 7) THEN c CALL LIB$EMUL(134775813,SEED,1,SEED) !Turbo Pascal 6.0 c IRANDOM = SEED(1) .AND. 'FFFF0000'X c IRANDOM = JISHFT(IRANDOM,-16) c IRANDOM = INT( FLOAT(IRANDOM)/(TWO16F) * FNBINSPD) c ELSEIF (MRANDO .EQ. 8) THEN c IRANDOM = INT( RANDES() * FNBINSPD ) !DES (not supplied) c ELSEIF (MRANDO .EQ. 9) THEN IRANDOM = INT( RAN1(SEED(1)) * FNBINSPD ) c ENDIF C Calculate index by hand. INDEX = INDEX + IRANDOM*(NBINSPD**I) ENDDO BINS(INDEX) = BINS(INDEX) + 1 !ball fell in this bin C IF ( MOD(J, 1 000 000) .EQ. 0 ) THEN C CALL TIME(TIMEBUF) C WRITE(6,302) J, TIMEBUF 302 FORMAT(1X,'AT BALL:',I,3X,'TIME=',A8) C WRITE(6,303) SEED(2), SEED(1) 303 FORMAT(1X,'HEX: SEED(2)= ',Z,' SEED(1)= ',Z) C WRITE(6,304) SEED(2), SEED(1) 304 FORMAT(1X,'DEC: SEED(2)= ',I,' SEED(1)= ',I) C ENDIF ENDDO 400 CALL CHSONE(BINS,EBINS,NBINS,CHSQ,PROB) SAVEPROB(CHITEST) = PROB WRITE(6,1) NBALLS,CHSQ,PROB 1 FORMAT(' BALLS=',I,' CHISQ=',F,' PROB=',F) ENDDO C Now see if all the chi-square values are chi-square distributed: IF (NCHITESTS .GT. 1) THEN CALL KSONE(SAVEPROB,NCHITESTS,D,PROB) WRITE(6,2) D,PROB 2 FORMAT(1X,'KS D=',F,' PROB=',F) ENDIF END C============================================================================ C From book NUMERICAL RECIPES: The Art of Scientific Computing C Here for demonstration purposes C Replace this with whatever random number generator you want to test C Initialize with negative number FUNCTION RAN1(IDUM) REAL R(97) SAVE R !(Some UNIX F77 compilers require -save option on compile) PARAMETER (M1=259200,IA1=7141,IC1=54773,RM1=3.8580247E-6) PARAMETER (M2=134456,IA2=8121,IC2=28411,RM2=7.4373773E-6) PARAMETER (M3=243000,IA3=4561,IC3=51349) DATA IFF /0/ IF (IDUM.LT.0.OR.IFF.EQ.0) THEN IFF=1 IX1=MOD(IC1-IDUM,M1) IX1=MOD(IA1*IX1+IC1,M1) IX2=MOD(IX1,M2) IX1=MOD(IA1*IX1+IC1,M1) IX3=MOD(IX1,M3) DO 11 J=1,97 IX1=MOD(IA1*IX1+IC1,M1) IX2=MOD(IA2*IX2+IC2,M2) R(J)=(FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1 11 CONTINUE IDUM=1 ENDIF IX1=MOD(IA1*IX1+IC1,M1) IX2=MOD(IA2*IX2+IC2,M2) IX3=MOD(IA3*IX3+IC3,M3) J=1+(97*IX3)/M3 IF(J.GT.97.OR.J.LT.1)PAUSE write(1,100) R write(1,102) R(J) 100 format(f) 102 format(1x,'RAN1 = ', F) RAN1=R(J) R(J)=(FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1 RETURN END C---------------------------------------------------------------------------- C CALCULATE THE CHI-SQUARE PROBABILITY. SINCE NBINS IS LARGE, IT IS JUST C THE CUMULATIVE GAUSSIAN DISTRIBUTION AFTER WE NORMALIZE THE VARIABLES. C OR ERROR FUNCTION. FUNCTION CHIPROB(NBINS,CHISQ) C Formula is the inverse of one given in Knuth for going the other way. INTEGER NBINS,DF REAL*4 CHISQ,Z DF = NBINS-1 Z = ( SQRT(24.0*CHISQ - 6.0*DF + 16.0) - 3*SQRT(2.0*DF) ) / 4.0 CHIPROB = ERF(Z) RETURN END FUNCTION ERF(X) C Return approximation to the complimentary error function erfc(X). C Return is not normalized err function. See book for details. C Adapted from book NUMERICAL RECIPES: The Art of Scientific Computing C Modified to return normalized error function erf(X) C (It's a polynomial approximation) REAL ERFCC,Z,T Z=ABS(X/1.414213) !Normalize T=1./(1.+0.5*Z) ERFCC=T*EXP(-Z*Z-1.26551223+T*(1.00002368+T*(.37409196+ * T*(.09678418+T*(-.18628806+T*(.27886807+T*(-1.13520398+ * T*(1.48851587+T*(-.82215223+T*.17087277))))))))) IF (X.LT.0.) ERFCC=2.-ERFCC ERF = 1.0 - ERFCC/2.0 !Normalize and compliment RETURN END C---------------------------------------------------------------------------- C THE FOLLOWING SUBROUTINES CALCULATE THE CHI-SQUARE VALUE: SUBROUTINE CHSONE(BINS,EBINS,NBINS,CHSQ,PROB) C Adapted from book NUMERICAL RECIPES: The Art of Scientific Computing INTEGER NBINS INTEGER*2 BINS(NBINS) REAL EBINS,CHSQ,PROB CHSQ=0. IF(EBINS.LE.0.) PAUSE 'CHSONE: EBINS must be > 0' DO 11 J=1,NBINS CHSQ=CHSQ+(BINS(J)-EBINS)**2/EBINS 11 CONTINUE PROB=CHIPROB(NBINS,CHSQ) RETURN END C============================================================================ C THE FOLLOWING SUBROUTINES CALCULATE THE KOLMOGOROV-SMIRNOV PROBABILITY SUBROUTINE KSONE(DATA,N,D,PROB) C Adapted from book NUMERICAL RECIPES: The Art of Scientific Computing C DF - degrees of freedom. Passsed to FUNC INTEGER N REAL DATA(N) REAL D,PROB CALL PIKSRT(N,DATA) EN=N D=0. FO=0. DO 11 J=1,N FN=J/EN FF=DATA(J) DT=AMAX1(ABS(FO-FF),ABS(FN-FF)) IF(DT.GT.D)D=DT FO=FN 11 CONTINUE PROB=PROBKS(SQRT(EN)*D) RETURN END C---------------------------------------------------------------------------- FUNCTION PROBKS(ALAM) C Adapted from book NUMERICAL RECIPES: The Art of Scientific Computing C Note the routine in the Numerical Recipes book erronously returns C 1 instead of 0 for large values of ALAM. PARAMETER (EPS1=0.001, EPS2=1.E-8) A2=-2.*ALAM**2 FAC=2. PROBKS=0. TERMBF=0. DO 11 J=1,100 TERM=FAC*EXP(A2*J**2) PROBKS=PROBKS+TERM C Error in Numerical Recipes book. Terminate if TERM underflows. C** IF(ABS(TERM).LT.EPS1*TERMBF.OR.ABS(TERM).LT.EPS2*PROBKS)RETURN IF(ABS(TERM).LE.EPS1*TERMBF.OR.ABS(TERM).LE.EPS2*PROBKS)RETURN FAC=-FAC TERMBF=ABS(TERM) 11 CONTINUE PROBKS=1.0 RETURN END C---------------------------------------------------------------------------- SUBROUTINE PIKSRT(N,ARR) C Adapted from book NUMERICAL RECIPES: The Art of Scientific Computing C See book for details. INTEGER N REAL ARR(N) DO 12 J=2,N A=ARR(J) DO 11 I=J-1,1,-1 IF(ARR(I).LE.A)GO TO 10 ARR(I+1)=ARR(I) 11 CONTINUE I=0 10 ARR(I+1)=A 12 CONTINUE RETURN END