! Last change: DOS 1 Aug 2000 8:11 pm ! *** copyright 2000 *** ! *** filename ghrua.f95 *** John F. Monahan ** ! ********************** PROGRAM PGHRUA ! TEST OF GHRUA -- RATIO OF UNIFORMS METHOD FOR GENERATING RANDOM ! VARIABLES FROM HYPERGEOMETRIC DISTRIBUTION ! IMPLICIT NONE REAL, DIMENSION(51) :: P,EXPECT REAL EXPCTD,U,CHISQ REAL RAN,RBICOF INTEGER, PARAMETER :: NOBS = 2500 ! SAMPLE SIZE INTEGER, DIMENSION(51) :: IP,OBS INTEGER N,NP1,I,L,OB,J,K,KCELL,GHRUA,NN,MM INTEGER, DIMENSION(3) :: NLIST,MLIST ! 20 FORMAT(' HYPERGEOMETRIC DISTRIBUTION'/' DRAWING',I4, & & ' BALLS FROM URN WITH',I4,' RED BALLS OUT OF TOTAL OF', & & I4,' BALLS'/' PROBABILITY OF J RED BALLS'/11F7.4/2X,11F7.4) 21 FORMAT(11I7) 22 FORMAT(' EXPECTED NUMBER FOR EACH CELL'/11F7.1/2X,11F7.1) 23 FORMAT(/' CHI-SQUARE TEST WITH',I4,' CELLS, STATISTIC =',F12.4//) 25 FORMAT(/' CELL COUNTS FOR ',I6,' OBSERVATIONS') ! ! DATA NLIST/ 100, 50, 40 / DATA MLIST/ 40, 30, 16 / ! OPEN( UNIT=6, FILE = 'ghrua.out' ) ! INITIALIZE UNIFORM GENERATOR U = RAN(5151917) ! DO SOME HYPERGEOMETRIC DISTRIBUTIONS ! N = NUMBER OF BALLS DRAWN FROM URN DO N = 5,15,5 DO L = 1,3 ! NN = TOTAL NUMBER OF BALLS IN URN NN = NLIST(L) ! MM = NUMBER OF RED BALLS IN URN MM = MLIST(L) ! FILL UP TABLE OF PROBABILITIES NP1 = N + 1 DO I = 1,NP1 P(I) = RBICOF(MM,I-1)*RBICOF(NN-MM,N-I+1) / RBICOF(NN,N) ! 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,MM,NN,(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 = GHRUA(NN,MM,N) ! 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 L END DO ! LOOP ON N STOP END PROGRAM PGHRUA INTEGER FUNCTION GHRUA(NN,NM,NNS) ! GENERATES HYPERGEOMETRIC RANDOM VARIABLES ! INPUT NN TOTAL NUMBER OF BALLS IN THE URN ! INPUT NM NUMBER OF RED BALLS ! INPUT NNS NUMBER OF BALLS TO BE DRAWN FROM THE URN ! ! ** 2 LE NNS LS N/2 *** 2 LE NM LE NN/2 ** NNS*NM GE NN *** ! ! OUTPUT GHRUA HYPERGEOMETRIC DEVIATE (NUMBER OF RED BALLS DRAWN) ! ! REQUIRED FUNCTION SUBPROGRAMS ! RAN UNIFORM RANDOM NUMBER GENERATOR ! SLGAMK COMPUTES LN K! IN DOUBLE PRECISION ! ! MODIFICATION OF ALGORITHM HRUAT FROM: ! ERNST STADLOBER, "SAMPLING FROM POISSON, BINOMIAL AND HYPERGEOMETRIC ! DISTRIBUTIONS: RATIO OF UNIFORMS AS A SIMPLE AND FAST ALTERNATIVE," ! TECHNISCHE UNIVERSITAT GRAZ, 1989. ! SEE ALSO: ! ERNST STADLOBER (1990) "THE RATIO OF UNIFORMS APPROACH FOR GENERATING ! DISCRETE RANDOM VARIABLES," JOURNAL OF COMPUTATIONAL AND APPLIED ! MATHEMATICS, VOLUME 31, PP. 181-189. ! ! JFM SEPT 1990 ! J F MONAHAN (REVISED JUNE 1997) DEPT OF STATISTICS, NC STATE UNIV ! RECODED FEBRUARY, MARCH 2000 FOR FORTRAN 95 ! IMPLICIT NONE INTEGER, INTENT(IN) :: NN,NM,NNS REAL, PARAMETER :: D1 = 1.71552777 ! D1 = 2*SQRT(2/E) REAL, PARAMETER :: D2 = .898916162 ! D2 = 3 - 2*SQRT(3/E) REAL, PARAMETER :: BB = 7. REAL U,T,X REAL RAN REAL, SAVE :: P,Q,A,C,H REAL(KIND=8), SAVE :: G REAL(KIND=8) :: SLGAMK INTEGER, SAVE :: MS,IB INTEGER, SAVE :: N = 0 INTEGER, SAVE :: M = 0 INTEGER, SAVE :: NS = 0 ! ! CHECK ON INITIALIZATIONS IF( (NN.NE.N) .OR. (NM.NE.M) .OR. (NNS.NE.NS) ) THEN N = NN M = NM NS = NNS P = REAL(M)/REAL(N) Q = 1.-P A = NS*P + .5 C = SQRT( (N-NS)*NS*P*Q/(N-1.) + .5 ) H = D1*C + D2 MS = ( (NS+1)*(M+1) )/ (N+2) G = SLGAMK(MS) + SLGAMK(M-MS) + SLGAMK(NS-MS) + SLGAMK(N-M-NS+MS) IB = MIN( MIN(NS,M), INT(A+BB*C) ) END IF ! CHANGE ! END OF INITIALIZATIONS ! DO ! UNRESTRICTED DO ! RETURN TO HERE FOR A NEW DEVIATE U = RAN(1) X = A + H*( RAN(2) -.5 ) / U IF( X .LE. 0.) CYCLE GHRUA = INT(X) IF( GHRUA .GT. IB ) CYCLE T = REAL( G - SLGAMK(GHRUA) - SLGAMK(M-GHRUA) & & - SLGAMK(NS-GHRUA)-SLGAMK(N-M-NS+GHRUA) ) ! QUICK ACCEPT CHECK IF( T .GE. U*(4.-U) - 3. ) RETURN ! QUICK REJECT CHECK IF( U*(U-T) .GE. 1. ) CYCLE ! MAIN RATIO OF UNIFORMS CHECK IF( 2.*LOG(U) .LE. T ) RETURN END DO ! END UNRESTRICTED DO RETURN END FUNCTION GHRUA 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 FUNCTION RBICOF(I,J) ! COMPUTES BINOMIAL COEFFICIENTS, USUALLY BY TABLE ! IMPLICIT NONE INTEGER, INTENT(IN) :: I,J INTEGER K,L,N INTEGER, DIMENSION(21) :: PUSH REAL, DIMENSION(122) :: BICFTB REAL(KIND=8) SLGAMK ! ! TEST WHETHER IT IS IN TABLE - ! REMEMBER THE RESULT WILL GO HAYWIRE IN INTEGERS IF TOO BIG ! DATA PUSH/1,2,4,6,8,11,14,18,22,27,32,38,44,51,58,66,74,83,92, & & 102,112/ DATA BICFTB/1.,1.,1.,1.,2.,1.,3.,1.,4.,6.,1.,5.,10.,1.,6.,15., & & 20.,1.,7.,21.,35.,1.,8.,28.,56.,70.,1.,9.,36.,84.,126.,1.,10., & & 45.,120.,210.,252.,1.,11.,55.,165.,330.,462.,1.,12.,66.,220., & & 495.,792.,924.,1.,13.,78.,286.,715.,1287.,1716.,1.,14.,91.,364.& &,1001.,2002.,3003.,3432.,1.,15.,105.,455.,1365.,3003.,5005.,6435.& &,1.,16.,120.,560.,1820.,4368.,8008.,11440.,12870.,1.,17.,136., & & 680.,2380.,6188.,12376.,19448.,24310., 1.,18.,153.,816.,3060., & & 8568.,18564.,31824.,43758.,48620., 1.,19.,171.,969.,3876.,11628.& &,27132.,50388.,75582.,92378., 1.,20.,190.,1140.,4845.,15504., & & 38760.,77520.,125970.,167960.,184756. / ! ! TAKE CARE OF SPECIAL CASES RBICOF = 0. IF( ( J .LT. 0 ) .OR. ( J .GT. I ) ) RETURN RBICOF = 1. IF( ( J .EQ. 0 ) .OR. ( J .EQ. I ) ) RETURN ! GENERAL CASE ! L = MIN(J,I-J) IF( I .LE. 20 ) THEN ! LOOK IT UP IN THE TABLE K = PUSH(I+1) ! FIND THE STARTING POINT FOR I+1 RBICOF = BICFTB(K+L) ! THEN GET THE L-TH ELEMENT ELSE IF( L .LE. 20 ) THEN ! FOR BIG VALUES COMPUTE BY MULTIPLICATION RBICOF = 1. DO N = 1,L RBICOF = RBICOF*REAL(I-N+1) / REAL(N) END DO ! LOOP ON N ELSE ! FOR LONG MULTIPLICATIONS, USE STIRLING VIA SLGAMK RBICOF = REAL( SLGAMK(I) - SLGAMK(J) - SLGAMK(I-K) ) RBICOF = EXP( RBICOF ) END IF ! ( L .LE. 20 ) END IF ! ( I .LE. 20 ) RETURN END FUNCTION RBICOF 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 ghrua.f95 *********************