! Last change: DOS 1 Aug 2000 8:05 pm ! *** copyright 2000 *** ! *** filename gdinv.f95 *** John F. Monahan ** ! ********************** PROGRAM GDINV ! TEST OF GDINV -- DISCRETE INVERSION FOR GENERATING RANDOM VARIABLES ! FROM FINITE DISCRETE DISTRIBUTION ! IMPLICIT NONE REAL, DIMENSION(51) :: P,Q,EXPECT REAL GN,ALGR,ALGONP,ALP,EXPCTD,U,CHISQ REAL, PARAMETER :: PP = .40 ! PROBABILITY FOR TEST REAL RAN REAL(KIND=8) SLGAMK INTEGER, PARAMETER :: NOBS = 1024 ! NUMBER OF OBS TO TEST INTEGER, DIMENSION(51) :: IP,OBS INTEGER N,NP1,I,OB,J,K,KCELL INTEGER IFIND ! 20 FORMAT(' N=',I4,' PROBS AND INDEX VECTOR'/11F7.4/2X,11F7.4/9F7.4) 21 FORMAT(11I7) 22 FORMAT(' EXPECTED NUMBER FOR EACH CELL'/11F7.1/2X,11F7.1/9F7.1) 23 FORMAT(/' CHI-SQUARE TEST WITH',I4,' CELLS, STATISTIC =',F12.4//) 24 FORMAT(' CUMULATIVE PROBS & CELL INDICES'/11F7.4/2X,11F7.4/9F7.4) 25 FORMAT(/' CELL COUNTS FOR ',I6,' OBSERVATIONS') ! ! OPEN( UNIT=6, FILE = 'gdinv.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) ! SORT THE PROBABILITIES CALL HKSORT(P,IP,NP1) WRITE(6,20) N,(P(J),J=1,NP1) WRITE(6,21) (IP(J),J=1,NP1) ! GET THE CUMULATIVES Q(1) = P(1) DO I = 2,NP1 Q(I) = Q(I-1) + P(I) END DO ! LOOP ON I ! ! DO NOBS DRAWS FROM THIS DISTRIBUTION DO I = 1,NOBS U = RAN(I) ! RANDOM VARIABLE FROM INVERSION J = IFIND(U,Q,NP1) + 1 ! IP(J) IS REAL VALUE K = IP(J) ! 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 ! 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 GDINV INTEGER FUNCTION IFIND(X,XK,N) ! FIND I SUCH THAT XK(I) LE X LT XK(I+1) ! IFIND=0 IF X LT XK(1) ! IFIND=N IF X GT XK(N) ! J F MONAHAN JAN 1982 DEPT OF STAT, N C S U, RALEIGH, NC 27650 ! RECODED FEBRUARY, 2000 FOR FORTRAN 95 ! IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL, INTENT(IN) :: X REAL, DIMENSION(N), INTENT(IN) :: XK INTEGER IL,IU ! SPECIAL CASES IFIND = 0 IF( X .LT. XK(1) ) RETURN IFIND = N IF( X .GE. XK(N) ) RETURN IL = 1 IU = N DO WHILE( IU-IL .GT. 1 ) IFIND = ( IU + IL ) / 2 IF( X .EQ. XK(IFIND) ) RETURN IF( X .LT. XK(IFIND) ) THEN IU = IFIND ELSE IL = IFIND END IF ! ( X .LT. XK(IFIND) ) END DO ! WHILE ( IU-IL .GT. 1 ) IFIND = IL RETURN END FUNCTION IFIND 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 gdinv.f95 *********************