! Last change: DOS 1 Aug 2000 8:03 pm ! *** copyright 2000 *** ! *** filename gbrus.f95 *** John F. Monahan ** ! ********************** PROGRAM PGBRUS ! TEST OF GBRUS -- RATIO OF UNIFORMS METHOD FOR GENERATING RANDOM ! VARIABLES FROM BINOMIAL DISTRIBUTION ! IMPLICIT NONE REAL, PARAMETER :: PP = .40 ! BINOMIAL PROBABILITY REAL, DIMENSION(51) :: P,EXPECT REAL GN,ALGR,ALGONP,ALP,EXPCTD,U,CHISQ REAL RAN REAL(KIND=8) SLGAMK INTEGER, PARAMETER :: NOBS = 2500 ! SAMPLE SIZE INTEGER, DIMENSION(51) :: IP,OBS INTEGER N,NP1,I,OB,J,K,KCELL INTEGER GBRUS ! 20 FORMAT(' N =',I4,' PROBABILITIES'/11F7.4/2X,11F7.4/4X,4F7.4) 21 FORMAT(11I7) 22 FORMAT(' EXPECTED NUMBER FOR EACH CELL'/11F7.1/2X,11F7.1/4F7.1) 23 FORMAT(/' CHI-SQUARE TEST WITH',I4,' CELLS, STATISTIC =',F12.4//) 25 FORMAT(/' CELL COUNTS FOR ',I6,' OBSERVATIONS') ! ! OPEN( UNIT=6, FILE = 'gbrus.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) WRITE(6,21) (IP(J),J=1,NP1) ! ! DO NOBS DRAWS FROM THIS DISTRIBUTION DO I = 1,NOBS ! RANDOM VARIABLE FROM GBRUS K = GBRUS(N,PP) ! UPDATE CELL COUNT OBS(K+1) = OBS(K+1) + 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 PGBRUS INTEGER FUNCTION GBRUS(N,P) REAL, INTENT(IN) :: P INTEGER, INTENT(IN) :: N ! RATIO OF UNIFORMS ALGORITHM FOR GENERATING FROM BINOMIAL DISTRIBUTION ! WITH N TRIALS WITH PROBABILITY P ! INPUT ! N INTEGER NUMBER OF TRIALS ! P REAL PROBABAILITY OF SUCCESS IN EACH TRIAL !*** *** P .LE. 1/2 *** *** NOTE *** ! ! ADAPTED FROM ALGORITHM BRUS OF ! ERNST STADLOBER (1991) "BINOMIAL RANDOM VARIATE GENERATION: A METHOD ! BASED ON RATIO OF UNIFORMS," IN 'THE FRONTIERS OF STATISTICAL ! COMPUTATION, SIMULATION & MODELING,' ED BY PETER R. NELSON, ! AMERICAN SCIENCES PRESS. ! ! J F MONAHAN (JUNE 1997) DEPT OF STATISTICS, NC STATE UNIV ! RECODED FEBRUARY, MARCH 2000 FOR FORTRAN 95 REAL, PARAMETER :: C1 = .8577639 ! C1 = SQRT(2/E) REAL, PARAMETER :: C2 = .4494581 ! C2 = 3/2 - SQRT(3/E) REAL U,X,T REAL, SAVE :: PS = 0. REAL, SAVE :: A,SNP REAL(KIND=8), SAVE :: R,G REAL(KIND=8) SLGAMK INTEGER, SAVE :: NS = 0 INTEGER, SAVE :: M,NB ! DO WE HAVE TO DO THE SET-UP? IF( (N .NE. NS) .OR. (P .NE. PS) ) THEN ! ! SET-UP CALCULATIONS NS = N PS = P A = REAL(N)*P + .5 U = SQRT( REAL(N)*P*(1.-P) + .5 ) SNP = C1*U + C2 M = INT( REAL(N+1)*P ) R = LOG( REAL(P,8)/(1.D0-REAL(P,8)) ) !KIND G = SLGAMK(M) + SLGAMK(N-M) ! SET UPPER LIMIT OF DISTRIBUTION ON N AND ROUNDOFF NB = MIN( N, INT( A + 7.*U ) ) END IF ! (DO SET-UP) ! END OF SET-UP DO ! UNRESTRICTED DO U = RAN(1) X = A + SNP*(2.*RAN(2)-1.)/U IF( X .LE. 0. ) CYCLE ! REJECT -- OUT OF RANGE GBRUS = INT(X) IF( GBRUS .GT. NB ) CYCLE ! REJECT -- OUT OF RANGE T = REAL( (GBRUS-M)*R + G - SLGAMK(GBRUS) - SLGAMK(NS-GBRUS) ) ! QUICK ACCEPT IF( U*(4.-U) - 3. .LE. T ) RETURN ! QUICK REJECT IF( U*(U-T) .GE. 1. ) CYCLE ! REJECT ! MAIN RATIO OF UNIFORMS CHECK IF( 2.*LOG(U) .LE. T ) RETURN END DO ! UNRESTRICTED DO END FUNCTION GBRUS 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 gbrus.f95 *********************