! Last change: DOS 1 Aug 2000 8:23 pm ! *** copyright 2000 *** ! *** filename gprua.f95 *** John F. Monahan ** ! ********************** PROGRAM PGPRUA ! TEST OF GPRUA -- DISCRETE RATIO OF UNIFORMS ALGORITHM ! FOR GENERATING RANDOM VARIABLES FROM POISSON DISTRIBUTION ! IMPLICIT NONE REAL RMU,RMULOG,P,PLOG,SUMTOJ,PSOFAR,U,CHISQ REAL, DIMENSION(251) :: EXPECT REAL RAN REAL(KIND=8) SLGAMK INTEGER M,N,NP1,I,J,K,LUP INTEGER, PARAMETER :: NOBS = 1000 ! SAMPLE SIZE INTEGER, DIMENSION(251) :: IQ,OBS INTEGER GPRUA,IJFIND ! 20 FORMAT(/' TEST OF GPRUA WITH MEAN=',F9.2,8X,I8,' OBSERVATIONS') 21 FORMAT(13I6) 22 FORMAT(' IN TAIL, OBS',I4,' U=',F6.4,' INDEX',I3,' REAL X',I3,& & ' CELL',I3) 23 FORMAT(/' CHI-SQUARE TEST WITH',I4,' CELLS, STATISTIC =',F12.4//) 24 FORMAT(' EXPECTEDS AND CELL INDICES'/13F6.1/2X,13F6.1) 25 FORMAT(/' CELL COUNTS FOR ',I6,' OBSERVATIONS') ! ! OPEN( UNIT=6, FILE = 'gprua.out' ) ! INITIALIZE UNIFORM GENERATOR U = RAN(5151917) ! DO SOME POISSON DISTRIBUTIONS -- CHANGE MEAN DO M = 1,20 RMU = REAL(M*M) WRITE(6,20) RMU,NOBS ! NOW DO CHI-SQUARE TEST ! ! INITIALIZE CELLS ! ONLY MAKE CELLS BIG ENOUGH TO HAVE 10 EXPECTEDS K = 1 RMULOG = LOG(RMU) SUMTOJ = 0. PSOFAR = 0. ! SET UPPER LIMIT FOR NOW MEAN + 5*STDEV LUP = M*(M+5) DO J = 0,LUP ! GET PR( X = J ) = EXP(-RMU) * ( RMU**J ) / J! P = 0. PLOG = J*RMULOG - RMU - SLGAMK(J) IF( PLOG .GT. -50. ) P = EXP(PLOG) SUMTOJ = SUMTOJ + P EXPECT(K) = ( SUMTOJ - PSOFAR ) * NOBS IF( EXPECT(K) .GE. 10. ) THEN ! BIG ENOUGH TO COUNT -- MAKE A NEW CELL IQ(K) = J OBS(K) = 0 PSOFAR = SUMTOJ K = K + 1 END IF ! ( EXPECT(K) .GT. 10. ) END DO ! LOOP ON J N = K - 1 NP1 = N + 1 EXPECT(NP1) = NOBS*(1.-PSOFAR) OBS(NP1) = 0 IF( EXPECT(NP1) .LT. 10. ) THEN EXPECT(N) = EXPECT(N) + EXPECT(NP1) NP1 = N N = N - 1 END IF ! ( EXPECT(NP1) .LT. 10. ) WRITE(6,24) (EXPECT(J),J=1,NP1) WRITE(6,21) (IQ(J),J=1,N) ! DO NOBS DRAWS DO I = 1,NOBS ! RANDOM VARIABLE FROM GPRUA K = GPRUA(RMU) J = IJFIND(K,IQ,N) OBS(J) = OBS(J) + 1 END DO ! LOOP ON I ! WRITE OUT CELL COUNTS WRITE(6,25) NOBS WRITE(6,21) (OBS(J),J=1,NP1) ! COMPUTE CHI-SQUARE STATISTIC CHISQ = 0. DO J = 1,NP1 CHISQ = CHISQ + (OBS(J)-EXPECT(J)) * (OBS(J)-EXPECT(J))/EXPECT(J) END DO ! LOOP ON J WRITE(6,23) NP1,CHISQ END DO ! LOOP ON M STOP END PROGRAM PGPRUA INTEGER FUNCTION GPRUA(MU) ! GENERATES POISSON RANDOM VARIABLES ! INPUT MU MEAN OF POISSON DISTRIBUTION ! OUTPUT GPRUA POISSON DEVIATE ! ! REQUIRED FUNCTION SUBPROGRAMS ! RAN UNIFORM RANDOM NUMBER GENERATOR ! SLGAMK COMPUTES LN K! IN DOUBLE PRECISION ! ! MODIFICATION BY STADLOBER OF AN ALGORITHM BY AHRENS AND DIETER ! E. STADLOBER, SAMPLING FROM POISSON, BINOMIAL AND HYPERGEOMETRIC ! DISTRIBUTIONS: RATIO OF UNIFORMS AS A SIMPLE AND FAST ALTERNATIVE ! JFM SEPT 1990 ! RECODING FOR CLARITY -- JULY 1996 ! RECODING FOR FORTRAN 95 FEBRUARY, MARCH 2000 IMPLICIT NONE REAL, INTENT(IN) :: MU REAL SRA,U,X REAL RAN REAL, PARAMETER :: C = 1.71552777 ! C = 2*SQRT(2/E) REAL, PARAMETER :: D = .898916162 ! D = 3 - 2*SQRT(3/E) REAL, PARAMETER :: BB = 7. REAL, SAVE :: SMU = 0. REAL, SAVE :: A,H REAL(KIND=8) T,SLGAMK REAL(KIND=8), SAVE :: G,Q INTEGER M,K INTEGER, SAVE :: IBIG ! ! CHECK ON INITIALIZATIONS IF( MU .NE. SMU ) THEN SMU = MU A = SMU + .5 SRA = SQRT(A) H = C*SRA + D G = LOG( REAL(SMU,8) ) !KIND M = INT(SMU) Q = M*G - SLGAMK(M) ! PROB IS ZERO BEYOND IBIG, CANCELLATION ERROR BIGGER IBIG = INT( A + BB*(SRA+1.) ) END IF ! ( MU .NE. SMU ) ! END OF INITIALIZATIONS DO ! UNRESTRICTED DO U = RAN(1) X = A + H*( RAN(2) -.5 ) / U IF( (X .LE. 0.) .OR. (X .GE. IBIG) ) CYCLE ! IN RANGE CHECK K = INT(X) GPRUA = K T = K*G - SLGAMK(K) - Q IF( REAL(T) .GE. U*(4.-U) - 3. ) RETURN ! INNER CHECK IF( U*(U-REAL(T)) .GE. 1. ) CYCLE ! OUTER (REJECT) CHECK IF( 2.*LOG(U) .LT. REAL(T) ) RETURN END DO ! END UNRESTRICTED DO RETURN END FUNCTION GPRUA INTEGER FUNCTION IJFIND(J,ILIST,N) ! FIND I SUCH THAT ILIST(I-1) LT J LE ILIST(I) ! ! ARGUMENTS ! J INTEGER INPUT ITEM WHOSE PLACE IS TO BE FOUND IN LIST ! ILIST INTEGER VECTOR INPUT LIST OF INTEGER INDICES ! ** ILIST MUST BE IN INCREASING ORDER ** ! N INTEGER INPUT LENGTH OF LIST ILIST ! ! IJFIND = 1 IF J LE ILIST(1) ! IJFIND = N+1 IF J GT ILIST(N) ! ! J F MONAHAN JULY 1996 DEPT OF STAT, N C STATE U, RALEIGH, NC 27695 ! ADAPTION TO INTEGER CASE OF IFIND AND JFIND ! RECODED FEBRUARY, 2000 FOR FORTRAN 95 INTEGER, INTENT(IN) :: J,N INTEGER, DIMENSION(N) :: ILIST INTEGER IL,IU ! TAKE CARE OF SPECIAL CASES FIRST IJFIND = 1 IF( J .LE. ILIST(1) ) RETURN IJFIND = N + 1 IF( J .GT. ILIST(N) ) RETURN ! GENERAL CASE IL = 1 IU = N DO WHILE( IU-IL .GT. 1 ) IJFIND = (IU+IL)/2 ! DISCRETE BISECTION SEARCH ON INTEGERS IF( J .EQ. ILIST(IJFIND) ) RETURN IF( J .LT. ILIST(IJFIND) ) THEN ! I TOO BIG, MOVE UPPER LIMIT DOWN IU = IJFIND ELSE ! I TOO SMALL, MOVE LOWER LIMIT UP IL = IJFIND END IF ! ( J .LT. ILIST(IJFIND) ) END DO ! WHILE( IU-IL .GT. 1 ) ! INTERVAL LENGTH DOWN TO 1, FINISH IT OFF IJFIND = IU RETURN END FUNCTION IJFIND 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 gprua.f95 *********************