!*********************************************************************************************************************************** ! ! D E V I A T E S ! ! Module: DEVIATES ! ! Programmer: David G. Simpson ! NASA Goddard Space Flight Center ! Greenbelt, Maryland 20771 ! ! Date: March 6, 2008 ! ! Language: Fortran-90 ! ! Version: 1.00a ! ! Description: Generate random deviates for a variety of distributions. ! ! Files: Source files: ! ! deviates.f90 Main module ! ! Notes: Distributions included in this module: ! ! Beta ! Binomial (small t) ! Binomial (large t) ! Chi-squared ! Exponential ! F ! Gamma ! Geometric ! Normal ! Poisson (Ahrens and Dieter algorithm) ! Spherical ! t ! ! References: Knuth, Donald E. "Seminumerical Algorithms" (2nd ed.) (Section 3.4.1). Addison-Wesley, Reading, ! Massachusetts, 1981. ! ! Ahrens, J.H. and Dieter, U., "Computer Generation of Poisson Deviates from Modified Normal Distributions". ! ACM Trans. Math. Software, 8, 2, 163-179 (1982). ! !*********************************************************************************************************************************** MODULE DEVIATES IMPLICIT NONE CONTAINS !*********************************************************************************************************************************** ! BETA_DIST !*********************************************************************************************************************************** FUNCTION BETA_DIST (A, B) RESULT (X) IMPLICIT NONE INTEGER, INTENT(IN) :: A, B DOUBLE PRECISION :: X DOUBLE PRECISION :: X1, X2 X1 = GAMMA_DIST (A) X2 = GAMMA_DIST (B) X = X1/(X1+X2) RETURN END FUNCTION BETA_DIST !*********************************************************************************************************************************** ! BINOMIAL_DIST1 ! ! For small t. !*********************************************************************************************************************************** FUNCTION BINOMIAL_DIST1 (T, P) RESULT (N) IMPLICIT NONE INTEGER, INTENT(IN) :: T DOUBLE PRECISION, INTENT(IN) :: P INTEGER :: N, I DOUBLE PRECISION :: U N = 0 DO I = 1, T CALL RANDOM_NUMBER (U) IF (U .LT. P) N = N + 1 END DO RETURN END FUNCTION BINOMIAL_DIST1 !*********************************************************************************************************************************** ! BINOMIAL_DIST2 ! ! For large t. !*********************************************************************************************************************************** RECURSIVE FUNCTION BINOMIAL_DIST2 (T, P, CALLS, T2) RESULT (N) IMPLICIT NONE INTEGER, INTENT(IN) :: T DOUBLE PRECISION, INTENT(IN) :: P INTEGER, INTENT(IN OUT) :: CALLS, T2 INTEGER :: N INTEGER :: A, B DOUBLE PRECISION :: X INTEGER :: TARG DOUBLE PRECISION, PARAMETER :: LOG2 = 0.693147180559945309417232121458176568075500134360255254120680009493393621969694715606D0 CALLS = CALLS + 1 TARG = FLOOR(LOG(DBLE(T2))/LOG2) A = 1 + T/2 B = T-A+1 X = BETA_DIST (A, B) IF (X .GE. P) THEN IF (CALLS .LE. TARG) THEN N = BINOMIAL_DIST2 (A-1, P/X, CALLS, T2) ELSE N = BINOMIAL_DIST1 (A-1, P/X) END IF ELSE IF (CALLS .LE. TARG) THEN N = A + BINOMIAL_DIST2 (B-1, (P-X)/(1.0D0-X), CALLS, T2) ELSE N = A + BINOMIAL_DIST1 (B-1, (P-X)/(1.0D0-X)) END IF END IF RETURN END FUNCTION BINOMIAL_DIST2 !*********************************************************************************************************************************** ! CHI_SQUARED_DIST !*********************************************************************************************************************************** FUNCTION CHI_SQUARED_DIST (NU) RESULT (X) IMPLICIT NONE INTEGER, INTENT(IN) :: NU DOUBLE PRECISION :: X X = 2.0D0 * GAMMA_DIST (NU/2) RETURN END FUNCTION CHI_SQUARED_DIST !*********************************************************************************************************************************** ! EXPONENTIAL_DIST !*********************************************************************************************************************************** FUNCTION EXPONENTIAL_DIST (MU) RESULT (X) IMPLICIT NONE DOUBLE PRECISION, INTENT(IN) :: MU DOUBLE PRECISION :: X DOUBLE PRECISION :: U 10 CALL RANDOM_NUMBER (U) IF (U .EQ. 0.0D0) GO TO 10 X = -MU*LOG(U) RETURN END FUNCTION EXPONENTIAL_DIST !*********************************************************************************************************************************** ! F_DIST !*********************************************************************************************************************************** FUNCTION F_DIST (NU1, NU2) RESULT (X) IMPLICIT NONE INTEGER, INTENT(IN) :: NU1, NU2 DOUBLE PRECISION :: X DOUBLE PRECISION :: Y1, Y2 Y1 = CHI_SQUARED_DIST (NU1) Y2 = CHI_SQUARED_DIST (NU2) X = (Y1*NU1)/(Y2*NU2) RETURN END FUNCTION F_DIST !*********************************************************************************************************************************** ! GAMMA_DIST !*********************************************************************************************************************************** FUNCTION GAMMA_DIST (A) RESULT (X) IMPLICIT NONE INTEGER, INTENT(IN) :: A DOUBLE PRECISION :: X DOUBLE PRECISION, PARAMETER :: PI = 3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803D0 DOUBLE PRECISION :: U, V, Y, TARG 10 CALL RANDOM_NUMBER (U) Y = TAN(PI*U) X = SQRT(2.0D0*A)*Y+A-1 IF (X .LE. 0) GO TO 10 CALL RANDOM_NUMBER (V) TARG = (1.0D0+Y**2)*EXP((A-1.0D0)*LOG(X/(A-1.0D0))-SQRT(2.0D0*A-1.0D0)*Y) IF (V .GT. TARG) GO TO 10 RETURN END FUNCTION GAMMA_DIST !*********************************************************************************************************************************** ! GEOMETRIC_DIST !*********************************************************************************************************************************** FUNCTION GEOMETRIC_DIST (P) RESULT (N) IMPLICIT NONE DOUBLE PRECISION, INTENT(IN) :: P INTEGER :: N DOUBLE PRECISION :: U CALL RANDOM_NUMBER (U) N = CEILING(LOG(U)/LOG(1.0D0-P)) RETURN END FUNCTION GEOMETRIC_DIST !*********************************************************************************************************************************** ! NORMAL_DIST !*********************************************************************************************************************************** SUBROUTINE NORMAL_DIST (X1, X2) IMPLICIT NONE DOUBLE PRECISION, INTENT(OUT) :: X1, X2 DOUBLE PRECISION :: S, U1, U2, V1, V2 10 CALL RANDOM_NUMBER (U1) CALL RANDOM_NUMBER (U2) V1 = 2.0D0*U1-1.0D0 V2 = 2.0D0*U2-1.0D0 S = V1**2+V2**2 IF (S .GE. 1.0D0) GO TO 10 X1 = V1*SQRT(-2.0D0*LOG(S)/S) X2 = V2*SQRT(-2.0D0*LOG(S)/S) RETURN END SUBROUTINE NORMAL_DIST !*********************************************************************************************************************************** ! POISSON_DIST ! ! Ahrens and Dieter algorithm, from ACM Trans. Math. Software, 8, 2, 163-179 (1982). ! Downloaded from: http://www.netlib.org/toms/599 ! ! Converted to Fortran-90 by David G. Simpson 3/6/2008. ! ! FOR DETAILS SEE: ! ! AHRENS, J.H. AND DIETER, U. ! COMPUTER GENERATION OF POISSON DEVIATES ! FROM MODIFIED NORMAL DISTRIBUTIONS. ! ACM TRANS. MATH. SOFTWARE, 8,2 (JUNE 1982), 163 - 179. ! ! (SLIGHTLY MODIFIED VERSION OF THE PROGRAM IN THE ABOVE ARTICLE) ! ! INPUT: MU = MEAN MU OF THE POISSON DISTRIBUTION ! OUTPUT: KPOISS = SAMPLE FROM THE POISSON-(MU)-DISTRIBUTION ! !*********************************************************************************************************************************** FUNCTION POISSON_DIST (MU) RESULT (KPOISS) IMPLICIT NONE DOUBLE PRECISION, INTENT(IN) :: MU INTEGER :: KPOISS ! ! TABLES: COEFFICIENTS A(0)-A(7) FOR STEP F. FACTORIALS FACT ! COEFFICIENTS A(K) - FOR PX = FK*V*V*SUM(A(K)*V**K)-DEL ! DOUBLE PRECISION, DIMENSION(0:7), PARAMETER :: A = (/ -0.5D0, 0.3333333D0, -0.2500068D0, & 0.2000118D0, -0.1661269D0, 0.1421878D0, -0.1384794D0, 0.1250060D0 /) DOUBLE PRECISION, DIMENSION(10), PARAMETER :: FACT = (/ & 1.0D0, 1.0D0, 2.0D0, 6.0D0, 24.0D0, 120.0D0, 720.0D0, 5040.0D0, 40320.0D0, 362880.0D0 /) DOUBLE PRECISION, DIMENSION(35) :: PP INTEGER :: J, K, KFLAG, L, M DOUBLE PRECISION :: MUPREV, MUOLD, S, D, RN, RN2, G, FK, DIFMUK, U, OMEGA, B1, B2, C0, C1, C2, C3, & C, FX, FY, PX, PY, E, T, DEL, V, X, XX, P, P0, Q !----------------------------------------------------------------------------------------------------------------------------------- ! ! Initialize data. ! MUPREV = 0.0D0 ! PREVIOUS MU MUOLD = 0.0D0 ! MU AT LAST EXECUTION OF STEP P OR B. ! ! SEPARATION OF CASES A AND B ! IF (MU .EQ. MUPREV) GO TO 1 IF (MU .LT. 10.0D0) GO TO 12 ! ! C A S E A. (RECALCULATION OF S,D,L IF MU HAS CHANGED) ! MUPREV = MU S = SQRT(MU) D = 6.0D0*MU*MU ! ! THE POISSON PROBABILITIES PK EXCEED THE DISCRETE NORMAL ! PROBABILITIES FK WHENEVER K >= M(MU). L=INT(MU-1.1484) ! IS AN UPPER BOUND TO M(MU) FOR ALL MU >= 10 . ! L = INT(MU-1.1484D0) ! ! STEP N. NORMAL SAMPLE - SNORM(IR) FOR STANDARD NORMAL DEVIATE ! 1 CALL NORMAL_DIST (RN, RN2) G = MU + S*RN IF (G .LT. 0.0D0) GO TO 2 KPOISS = INT(G) ! ! STEP I. IMMEDIATE ACCEPTANCE IF KPOISS IS LARGE ENOUGH ! IF (KPOISS .GE. L) RETURN ! ! STEP S. SQUEEZE ACCEPTANCE - SUNIF(IR) FOR (0,1)-SAMPLE U ! FK = DBLE(KPOISS) DIFMUK = MU - FK CALL RANDOM_NUMBER (U) IF (D*U .GE. DIFMUK*DIFMUK*DIFMUK) RETURN ! ! STEP P. PREPARATIONS FOR STEPS Q AND H. ! (RECALCULATIONS OF PARAMETERS IF NECESSARY) ! .3989423=(2*PI)**(-.5) .416667E-1=1./24. .1428571=1./7. ! THE QUANTITIES B1, B2, C3, C2, C1, C0 ARE FOR THE HERMITE ! APPROXIMATIONS TO THE DISCRETE NORMAL PROBABILITIES FK. ! C=.1069/MU GUARANTEES MAJORIZATION BY THE 'HAT'-FUNCTION. ! 2 IF (MU .EQ. MUOLD) GO TO 3 MUOLD = MU OMEGA = 0.3989423D0/S B1 = 0.4166667D-1/MU B2 = 0.3D0*B1*B1 C3 = 0.1428571D0*B1*B2 C2 = B2-15.0D0*C3 C1 = B1-6.0D0*B2+45.0D0*C3 C0 = 1.0D0 - B1 + 3.0D0*B2 - 15.0D0*C3 C = 0.1069D0/MU 3 IF (G .LT. 0.0D0) GO TO 5 ! ! 'SUBROUTINE' F IS CALLED (KFLAG=0 FOR CORRECT RETURN) ! KFLAG = 0 GO TO 7 ! ! STEP Q. QUOTIENT ACCEPTANCE (RARE CASE) ! 4 IF (FY-U*FY .LE. PY*EXP(PX-FX)) RETURN ! ! STEP E. EXPONENTIAL SAMPLE - SEXPO(IR) FOR STANDARD EXPONENTIAL ! DEVIATE E AND SAMPLE T FROM THE LAPLACE 'HAT' ! (IF T <= -.6744 THEN PK < FK FOR ALL MU >= 10.) ! 5 E = EXPONENTIAL_DIST (0.0D0) CALL RANDOM_NUMBER (U) U = U + U - 1.0D0 T = 1.8D0 + SIGN(E,U) IF (T .LE. (-0.6744D0)) GO TO 5 KPOISS = INT(MU+S*T) FK = DBLE(KPOISS) DIFMUK = MU-FK ! ! 'SUBROUTINE' F IS CALLED (KFLAG=1 FOR CORRECT RETURN) ! KFLAG = 1 GO TO 7 ! ! STEP H. HAT ACCEPTANCE (E IS REPEATED ON REJECTION) ! 6 IF (C*ABS(U) .GT. PY*EXP(PX+E)-FY*EXP(FX+E)) GO TO 5 RETURN ! ! STEP F. 'SUBROUTINE' F. CALCULATION OF PX,PY,FX,FY. ! CASE KPOISS .LT. 10 USES FACTORIALS FROM TABLE FACT ! 7 IF (KPOISS .GE. 10) GO TO 8 PX = -MU PY = MU**KPOISS/FACT(KPOISS+1) GO TO 11 ! ! CASE KPOISS .GE. 10 USES POLYNOMIAL APPROXIMATION ! A0-A7 FOR ACCURACY WHEN ADVISABLE ! .8333333E-1=1./12. .3989423=(2*PI)**(-.5) ! 8 DEL = 0.8333333D-1/FK DEL = DEL - 4.8D0*DEL*DEL*DEL V = DIFMUK/FK IF (ABS(V) .LE. 0.25D0) GO TO 9 PX = FK*LOG(1.0D0+V) - DIFMUK - DEL GO TO 10 9 PX = FK*V*V*(((((((A(7)*V+A(6))*V+A(5))*V+A(4))*V+A(3))*V+A(2))*V+A(1))*V+A(0))-DEL 10 PY = 0.3989423D0/SQRT(FK) 11 X = (0.5D0-DIFMUK)/S XX = X*X FX = -0.5D0*XX FY = OMEGA*(((C3*XX+C2)*XX+C1)*XX+C0) IF (KFLAG) 4, 4, 6 ! ! C A S E B. (START NEW TABLE AND CALCULATE P0 IF NECESSARY) ! 12 MUPREV = 0.0D0 IF (MU .EQ. MUOLD) GO TO 13 MUOLD = MU M = MAX(1,INT(MU)) L = 0 P = EXP(-MU) Q = P P0 = P ! ! STEP U. UNIFORM SAMPLE FOR INVERSION METHOD ! 13 CALL RANDOM_NUMBER (U) KPOISS = 0 IF (U .LE. P0) RETURN ! ! STEP T. TABLE COMPARISON UNTIL THE END PP(L) OF THE ! PP-TABLE OF CUMULATIVE POISSON PROBABILITIES ! (0.458=PP(9) FOR MU=10) ! IF (L .EQ. 0) GO TO 15 J = 1 IF (U .GT. 0.458D0) J = MIN(L,M) DO K = J, L IF (U .LE. PP(K)) GO TO 18 END DO IF (L .EQ. 35) GO TO 13 ! ! STEP C. CREATION OF NEW POISSON PROBABILITIES P ! AND THEIR CUMULATIVES Q=PP(K) ! 15 L = L+1 DO K = L, 35 P = P*MU/DBLE(K) Q = Q + P PP(K) = Q IF (U .LE. Q) GO TO 17 END DO L = 35 GO TO 13 17 L = K 18 KPOISS = K RETURN END FUNCTION POISSON_DIST !*********************************************************************************************************************************** ! SPHERICAL_DIST !*********************************************************************************************************************************** SUBROUTINE SPHERICAL_DIST (D, XARR) IMPLICIT NONE INTEGER, INTENT(IN) :: D DOUBLE PRECISION, DIMENSION(D), INTENT(OUT) :: XARR INTEGER :: I DOUBLE PRECISION :: R2, X1, X2 R2 = 0.0D0 DO I = 1, D CALL NORMAL_DIST (X1, X2) XARR(I) = X1 R2 = R2 + X1**2 END DO XARR = XARR/SQRT(R2) RETURN END SUBROUTINE SPHERICAL_DIST !*********************************************************************************************************************************** ! T_DIST !*********************************************************************************************************************************** FUNCTION T_DIST (NU) RESULT (X) IMPLICIT NONE INTEGER, INTENT(IN) :: NU DOUBLE PRECISION :: X DOUBLE PRECISION :: X2, Y1, Y2 CALL NORMAL_DIST (Y1, X2) Y2 = CHI_SQUARED_DIST (NU) X = Y1/SQRT(Y2/DBLE(NU)) RETURN END FUNCTION T_DIST END MODULE DEVIATES