Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
real, | intent(in) | :: | s | |||
logical, | intent(in) | :: | first |
FUNCTION random_gamma2(s, first) RESULT(fn_val) ! Adapted from Fortran 77 code from the book: ! Dagpunar, J. 'Principles of random variate generation' ! Clarendon Press, Oxford, 1988. ISBN 0-19-852202-9 ! FUNCTION GENERATES A RANDOM VARIATE IN [0,INFINITY) FROM ! A GAMMA DISTRIBUTION WITH DENSITY PROPORTIONAL TO ! GAMMA2**(S-1) * EXP(-GAMMA2), ! USING A SWITCHING METHOD. ! S = SHAPE PARAMETER OF DISTRIBUTION ! (REAL < 1.0) REAL, INTENT(IN) :: s LOGICAL, INTENT(IN) :: first REAL :: fn_val ! Local variables REAL :: r, x, w REAL, SAVE :: a, p, c, uf, vr, d IF (s <= zero .OR. s >= one) THEN WRITE(*, *) 'SHAPE PARAMETER VALUE OUTSIDE PERMITTED RANGE' STOP END IF IF (first) THEN ! Initialization, if necessary a = one - s p = a/(a + s*EXP(-a)) IF (s < vsmall) THEN WRITE(*, *) 'SHAPE PARAMETER VALUE TOO SMALL' STOP END IF c = one/s uf = p*(vsmall/a)**s vr = one - vsmall d = a*LOG(a) END IF DO CALL RANDOM_NUMBER(r) IF (r >= vr) THEN CYCLE ELSE IF (r > p) THEN x = a - LOG((one - r)/(one - p)) w = a*LOG(x)-d ELSE IF (r > uf) THEN x = a*(r/p)**c w = x ELSE fn_val = zero RETURN END IF CALL RANDOM_NUMBER(r) IF (one-r <= w .AND. r > zero) THEN IF (r*(w + one) >= one) CYCLE IF (-LOG(r) <= w) CYCLE END IF EXIT END DO fn_val = x RETURN END FUNCTION random_gamma2