! Last change: DOS 1 Aug 2000 8:29 pm ! *** copyright 2000 *** ! *** filename walker.f95 *** John F. Monahan ** ! ********************** PROGRAM WALKER ! TEST OF WALKER'S METHOD FOR GENERATING RANDOM VARIABLES ! FROM FINITE DISCRETE DISTRIBUTION ! REAL, DIMENSION(51) :: P,Q,EXPECT REAL, PARAMETER :: PP = .40 ! BINOMIAL PROBABILITY REAL GN,ALGR,ALGONP,ALP,EXPCTD,U,CHISQ REAL(KIND=8) SLGAMK INTEGER, DIMENSION(51) :: IP,OBS,ALIAS INTEGER, PARAMETER :: NOBS = 1024 ! SAMPLE SIZE TO TEST INTEGER N,NP1,I,OB,J,K,KCELL ! 20 FORMAT(' N =',I4,' PROBABILITIES'/11F7.4/2X,11F7.4/4F7.4) 21 FORMAT(13I6) 22 FORMAT(' EXPECTED NUMBER FOR EACH CELL'/13F6.1/2X,13F6.1) 23 FORMAT(/' CHI-SQUARE TEST WITH',I4,' CELLS, STATISTIC =',F12.4//) 25 FORMAT(/' CELL COUNTS FOR ',I6,' OBSERVATIONS') 26 FORMAT(' ALIAS PROBABILITIES'/10F7.4/2X,10F7.4/4X,10F7.4) 27 FORMAT(9X,'ALIAS LIST FOR N=',I4/25I3) 28 FORMAT(' CHECK OF DISTRIBUTION: MAX DIFF OF PROBS',F8.5) ! ! OPEN( UNIT=6, FILE = 'walker.out' ) ! INITIALIZE UNIFORM GENERATOR U = RAN(5151917) ! DO SOME BINOMIAL DISTRIBUTIONS -- SAME P, CHANGE N DO N = 5,25,4 ! FILL UP TABLE OF PROBABILITIES NP1 = N + 1 GN = SLGAMK(N) ALGR = LOG(PP/(1.-PP)) ALGONP = LOG( 1. - PP ) DO I = 1,NP1 ALP = GN - SLGAMK(I-1) - SLGAMK(N+1-I) + (I-1)*ALGR + N*ALGONP P(I) = 0. IF( ALP .GT. -50. ) P(I) = EXP(ALP) ! CARRY AUXILIARY VECTOR -- INDEXED FROM ZERO IP(I) = I-1 ! INITIALIZE CELL COUNT TO ZERO OBS(I) = 0 ! GET EXPECTED VALUES EXPECT(I) = REAL(NOBS)*P(I) END DO ! LOOP ON I ! WRITE OUT PROBABILITIES WRITE(6,20) N,(P(J),J=1,NP1) ! CALL SET-UP ALGORITHM CALL WALKST(NP1,P,Q,ALIAS) ! WRITE OUT TARGET IF( N .LT. 25 ) WRITE(6,26) (Q(I),I=1,NP1) IF( N .LT. 25 ) WRITE(6,27) NP1,(ALIAS(I),I=1,NP1) ! CHECK THE DISTRIBUTION DO I = 1,NP1 P(I) = P(I) - Q(I)/REAL(NP1) J = ALIAS(I) IF( J .NE. 0 ) P(J) = P(J) - (1.-Q(I))/REAL(NP1) END DO ! IF WALKST IS RIGHT THEN P'S SHOULD NOW BE ZERO CHISQ = 0. DO I = 1,NP1 CHISQ = MAX( CHISQ, ABS(P(I)) ) END DO WRITE(6,28) CHISQ ! ! DO NOBS DRAWS FROM THIS DISTRIBUTION DO I = 1,NOBS U = RAN(I) K = INT( REAL(NP1)*RAN(I+1) ) + 1 IF( U .GT. Q(K) ) K = ALIAS(K) ! UPDATE CELL COUNT OBS(K) = OBS(K) + 1 END DO ! LOOP ON I ! WRITE OUT CELL COUNTS AND EXPECTEDS WRITE(6,25) NOBS WRITE(6,21) (OBS(J),J=1,NP1) WRITE(6,22) (EXPECT(J),J=1,NP1) ! ! NOW DO CHI-SQUARE TEST ! ! SORT THE PROBABILITIES CALL HKSORT(P,IP,NP1) ! CHISQ = 0. EXPCTD = 0. OB = 0 KCELL = 0 DO J = 1,NP1 K = IP(J) + 1 ! GO THROUGH CELLS FROM SMALLEST AND ! MAKE SURE CELLS BIG ENOUGH TO HAVE 10 EXPECTEDS EXPCTD = EXPCTD + EXPECT(K) OB = OB + OBS(K) IF( EXPCTD .GE. 10. ) THEN ! BIG ENOUGH TO COUNT CHISQ = CHISQ + (OB-EXPCTD) * (OB-EXPCTD)/EXPCTD KCELL = KCELL + 1 EXPCTD = 0. OB = 0 END IF ! ( EXPCTD .GE. 10. ) END DO ! LOOP ON J WRITE(6,23) KCELL,CHISQ ! END DO ! LOOP ON N STOP END PROGRAM WALKER SUBROUTINE WALKST(N,P,F,ALIAS) ! ROUTINE TO COMPUTE NEW PROBABILITIES AND ALIASES FOR WALKER'S ALIAS ! METHOD FOR GENERATING RANDOM VARIABLES FROM FINITE, DISCRETE ! DISTRIBUTION WITH PROBABILITY VECTOR P ! ! ARGUMENTS ! N IN INTEGER NUMBER OF ATOMS OF DISCRETE DISTRIBUTION ! P IN REAL VECTOR OF PROBABILITIES ! F OUT REAL VECTOR OF PROBABILITIES FOR ALIAS METHOD ! ALIAS OUT INTEGER VECTOR OF ALIASES ! ! REFERENCES ! ALASTAIR J. WALKER (1977) "AN EFFICIENT METHOD FOR GENERATING RANDOM ! VARIABLES WITH GENERAL DISTRIBUTIONS," ACM TRANSACTIONS ON ! MATHEMATICAL SOFTWARE, VOLUME 3, PP. 253-256. ! ! RICHARD A. KRONMAL AND ARTHUR V. PETERSON (1979) "ON THE ALIAS METHOD ! FOR GENERATING RANDOM VARIABLES FROM A DISCRETE DISTRIBUTION," THE ! AMERICAN STATISTICIAN, VOLUME 33, PP. 214-218. ! ! J F MONAHAN (JULY, 1996) DEPT OF STATISTICS, N C STATE U, 27695-8203 ! RECODED FEBRUARY, 2000 FOR FORTRAN 95 ! INTEGER, INTENT(IN) :: N INTEGER, DIMENSION(N), INTENT(OUT) :: ALIAS REAL, DIMENSION(N), INTENT(IN) :: P REAL, DIMENSION(N), INTENT(OUT) :: F INTEGER PXS,PND,I,J,K REAL FN ! INITIALIZATIONS FN = REAL(N) PXS = 0 PND = 0 ! PXS POINTS TO EXCESS CELL, PND TO NEEDY CELL DO I = 1,N F(I) = P(I) * FN ALIAS(I) = 0 ! NEEDY, EXCESS OR OK ? IF( F(I) .LT. 1. ) CALL ADDXS(I) IF( F(I) .GT. 1. ) CALL ADDND(I) DO WHILE( ( PND .NE. 0 ) .AND. (PXS .NE. 0 ) ) ! HAVE NEEDY AND EXCESS -- GET READY TO PROCESS J = PXS K = PND CALL PROCESS ! IS K STILL NEEDY ? IF( F(K) .LT. 1. ) THEN ! K NOW EXCESS, REMOVE FROM NEEDY LIST PND = ALIAS(K) CALL ADDXS(K) END IF ! ( F(K) .LT. 1. ) END DO ! WHILE ( HAVE BOTH NEEDY AND EXCESS ) ! K SATISFIED, DONE WITH IT IF( (PXS .EQ. 0 ) .AND. (PND .EQ. 0 ) ) ALIAS(K) = 0 END DO ! LOOP ON I RETURN CONTAINS SUBROUTINE ADDXS(I) INTEGER I ! ADD TO EXCESS LIST ALIAS(I) = PXS PXS = I J = I RETURN END SUBROUTINE ADDXS SUBROUTINE ADDND(I) INTEGER I ! ADD TO THE NEEDY LIST ALIAS(I) = PND PND = I K = I RETURN END SUBROUTINE ADDND SUBROUTINE PROCESS ! PROCESS A NEEDY AND AN EXCESS ! FIRST REMOVE CONTRIBUTOR PXS = ALIAS(J) ALIAS(J) = K ! USE THE CONTRIBUTOR'S PROBABILITY F(K) = F(K) - ( 1. - F(J) ) RETURN END SUBROUTINE PROCESS END SUBROUTINE WALKST SUBROUTINE HKSORT(K,M,N) ! HEAPSORT ALGORITHM FOR SORTING ON VECTOR OF KEYS K OF LENGTH N ! ! ARGUMENTS ! K REAL VECTOR OF KEYS TO BE SORTED ! M INTEGER VECTOR TO BE MOVED IN PARALLEL TO K ! N NUMBER OF ITEMS TO BE SORTED ! ! TO SORT A PARALLEL VECTOR OF RECORDS, USE THIS ONE ! ! J F MONAHAN (DEC, 1999) FORTRAN 95 IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL, DIMENSION(N), INTENT(IN OUT) :: K INTEGER, DIMENSION(N), INTENT(IN OUT) :: M REAL KK INTEGER I,L,NCUR,MM ! ! DO NOTHING IF THERE'S NOTHING TO DO IF( N .LE. 1 ) RETURN ! INITIALIZE TO BUILDHEAP PART (LOOP ON L) L = N/2 + 1 NCUR = N DO I = L,1,-1 CALL HEAPIFY(I) END DO ! LOOP ON I DO I = 2,N ! SWITCH CURRENT LARGEST WITH BOTTOM KK = K(1) K(1) = K(NCUR) K(NCUR) = KK ! ALSO SWITCH PARALLEL VECTOR MM = M(1) M(1) = M(NCUR) M(NCUR) = MM ! REHEAP WITH ONE SHORTER NCUR = NCUR - 1 CALL HEAPIFY(1) END DO ! LOOP ON NCUR RETURN CONTAINS SUBROUTINE HEAPIFY(II) INTEGER, INTENT(IN) :: II INTEGER I,J I = II DO J = 2*I ! IS IT A LEAF OR ARE THERE SONS? IF( J > NCUR ) EXIT ! A LEAF IF( J < NCUR ) THEN ! ANOTHER SON OF I IF( K(J+1) > K(J) ) J = J+1 ! LARGER SON IS K END IF ! A LEAF IF( K(J) > K(I) ) THEN ! EXCHANGE KK = K(J) K(J) = K(I) K(I) = KK ! ALSO SWITCH PARALLEL VECTOR MM = M(J) M(J) = M(I) M(I) = MM I = J ELSE ! EXIT -- HEAP PROPERTY EXIT END IF END DO ! WHILE NOT A LEAF END SUBROUTINE HEAPIFY END SUBROUTINE HKSORT REAL FUNCTION RAN(IDUM) ! PORTABLE IMPLEMENTATION OF UNIFORM PSEUDORANDOM NUMBER GENERATOR ! LEWIS, GOODMAN, & MILLER MULTIPLICATIVE GENERATOR ! X(N+1) = MOD( 16807*X(N), 2**31-1 ) ! ! P. BRANTLEY, B.L. FOX, L. SCHRAGE (1983) A GUIDE TO SIMULATION ! SPRINGER-VERLAG, NEW YORK. PP. 200-202. ! UPDATED VERSION OF ! LINUS SCHRAGE (1979)'A MORE PORTABLE FORTRAN RANDOM NUMBER GENERATOR' ! ACM TRANSACTIONS ON MATHEMATICAL SOFTWARE, VOLUME 5, PP. 132-138. ! ! ARGUMENT ! IDUM INTEGER FIRST CALL SETS SEED, IGNORED IN SUBSEQUENT CALLS ! IMPLICIT NONE INTEGER, INTENT(IN) :: IDUM REAL, PARAMETER :: RP = 4.656612875E-10 ! 1/P INTEGER, PARAMETER :: A = 16807 ! MULTIPLIER INTEGER, PARAMETER :: B = 127773 ! B = P / A INTEGER, PARAMETER :: C = 2836 ! C = P MOD A INTEGER, PARAMETER :: P = 2147483647 ! MODULUS 2**31 - 1 INTEGER, SAVE :: IX = 0 INTEGER K1 ! ! IF NOT FIRST CALL, THEN SKIP SETTING SEED IF( IX .EQ. 0) IX = IDUM ! WRITE NUMBER AS ALPHA*2**16 + BETA K1 = IX / B IX = A*( IX - K1*B) - K1*C ! ABOVE DOES A*IX MOD B -K1*C IF( IX .LT. 0 ) IX = IX + P ! RP BELOW IS RECIPROCAL OF P RAN = REAL(IX)*RP RETURN END FUNCTION RAN REAL(KIND=8) FUNCTION SLGAMK(K) ! PRODUCES NATURAL LOGARITHM OF K! ! FOR SMALL VALUES OF K -- USES TABLE ! FOR LARGER VALUES -- USES STIRLING'S APPROXIMATION ! ! J F MONAHAN (SEPT, 1990) DEPT OF STATISTICS, NCSU, RALEIGH, NC USA ! RECODED OCTOBER 1999, APRIL 2000 FOR FORTRAN 95 IMPLICIT NONE INTEGER, INTENT(IN) :: K REAL(KIND=8), DIMENSION(36) :: TAB REAL(KIND=8) F,DLF REAL(KIND=8), PARAMETER :: CHL2PI = .918938533205D0 ! log(2pi)/2 ! DATA TAB/ 0.D0, 0.6931471806D0, 1.7917594692D0, 3.1780538303D0& &, 4.7874917428D0, 6.5792512120D0, 8.5251613611D0,10.6046029027D0& &,12.8018274801D0,15.1044125731D0,17.5023078459D0,19.9872144957D0& &,22.5521638531D0,25.1912211827D0,27.8992713838D0,30.6718601061D0& &,33.5050734501D0,36.3954452080D0,39.3398841872D0,42.3356164608D0& &,45.3801388985D0,48.4711813518D0,51.6066755678D0,54.7847293981D0& &,58.0036052230D0,61.2617017610D0,64.5575386270D0,67.8897431372D0& &,71.2570389672D0,74.6582363488D0,78.0922235533D0,81.5579594561D0& &,85.0544670176D0,88.5808275422D0,92.1361756037D0,95.7196945421D0/ ! SLGAMK = 0.D0 IF( (K.EQ.0) .OR. (K.EQ.1) ) RETURN IF( K .LE. 36 ) THEN SLGAMK = TAB(K) RETURN ! USE STIRLING'S FORMULA ELSE F = DBLE(K) DLF = LOG(F) SLGAMK = ( (F+.5D0)*DLF - F ) & & + ( CHL2PI + ( 1.D0 - 1.D0/(30.D0*F*F) )/(12.D0*F) ) END IF ! ( K .LE. 36 ) ! RETURN END FUNCTION SLGAMK ! *** end of filename walker.f95 *********************