home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / octave-1.1.1p1-src.tgz / tar.out / fsf / octave / libcruft / ranlib / sgamma.f < prev    next >
Text File  |  1996-09-28  |  7KB  |  192 lines
  1.       REAL FUNCTION sgamma(a)
  2. C**********************************************************************C
  3. C                                                                      C
  4. C                                                                      C
  5. C     (STANDARD-)  G A M M A  DISTRIBUTION                             C
  6. C                                                                      C
  7. C                                                                      C
  8. C**********************************************************************C
  9. C**********************************************************************C
  10. C                                                                      C
  11. C               PARAMETER  A >= 1.0  !                                 C
  12. C                                                                      C
  13. C**********************************************************************C
  14. C                                                                      C
  15. C     FOR DETAILS SEE:                                                 C
  16. C                                                                      C
  17. C               AHRENS, J.H. AND DIETER, U.                            C
  18. C               GENERATING GAMMA VARIATES BY A                         C
  19. C               MODIFIED REJECTION TECHNIQUE.                          C
  20. C               COMM. ACM, 25,1 (JAN. 1982), 47 - 54.                  C
  21. C                                                                      C
  22. C     STEP NUMBERS CORRESPOND TO ALGORITHM 'GD' IN THE ABOVE PAPER     C
  23. C                                 (STRAIGHTFORWARD IMPLEMENTATION)     C
  24. C                                                                      C
  25. C     Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of   C
  26. C     SUNIF.  The argument IR thus goes away.                          C
  27. C                                                                      C
  28. C**********************************************************************C
  29. C                                                                      C
  30. C               PARAMETER  0.0 < A < 1.0  !                            C
  31. C                                                                      C
  32. C**********************************************************************C
  33. C                                                                      C
  34. C     FOR DETAILS SEE:                                                 C
  35. C                                                                      C
  36. C               AHRENS, J.H. AND DIETER, U.                            C
  37. C               COMPUTER METHODS FOR SAMPLING FROM GAMMA,              C
  38. C               BETA, POISSON AND BINOMIAL DISTRIBUTIONS.              C
  39. C               COMPUTING, 12 (1974), 223 - 246.                       C
  40. C                                                                      C
  41. C     (ADAPTED IMPLEMENTATION OF ALGORITHM 'GS' IN THE ABOVE PAPER)    C
  42. C                                                                      C
  43. C**********************************************************************C
  44. C
  45. C
  46. C     INPUT: A =PARAMETER (MEAN) OF THE STANDARD GAMMA DISTRIBUTION
  47. C     OUTPUT: SGAMMA = SAMPLE FROM THE GAMMA-(A)-DISTRIBUTION
  48. C
  49. C     COEFFICIENTS Q(K) - FOR Q0 = SUM(Q(K)*A**(-K))
  50. C     COEFFICIENTS A(K) - FOR Q = Q0+(T*T/2)*SUM(A(K)*V**K)
  51. C     COEFFICIENTS E(K) - FOR EXP(Q)-1 = SUM(E(K)*Q**K)
  52. C
  53.       DATA q1,q2,q3,q4,q5,q6,q7/.04166669,.02083148,.00801191,.00144121,
  54.      +     -.00007388,.00024511,.00024240/
  55.       DATA a1,a2,a3,a4,a5,a6,a7/.3333333,-.2500030,.2000062,-.1662921,
  56.      +     .1423657,-.1367177,.1233795/
  57.       DATA e1,e2,e3,e4,e5/1.,.4999897,.1668290,.0407753,.0102930/
  58. C
  59. C     PREVIOUS A PRE-SET TO ZERO - AA IS A', AAA IS A"
  60. C     SQRT32 IS THE SQUAREROOT OF 32 = 5.656854249492380
  61. C
  62.       DATA aa/0.0/,aaa/0.0/,sqrt32/5.656854/
  63. C
  64.       IF (a.EQ.aa) GO TO 10
  65.       IF (a.LT.1.0) GO TO 120
  66. C
  67. C     STEP  1:  RECALCULATIONS OF S2,S,D IF A HAS CHANGED
  68. C
  69.       aa = a
  70.       s2 = a - 0.5
  71.       s = sqrt(s2)
  72.       d = sqrt32 - 12.0*s
  73. C
  74. C     STEP  2:  T=STANDARD NORMAL DEVIATE,
  75. C               X=(S,1/2)-NORMAL DEVIATE.
  76. C               IMMEDIATE ACCEPTANCE (I)
  77. C
  78.    10 t = snorm()
  79.       x = s + 0.5*t
  80.       sgamma = x*x
  81.       IF (t.GE.0.0) RETURN
  82. C
  83. C     STEP  3:  U= 0,1 -UNIFORM SAMPLE. SQUEEZE ACCEPTANCE (S)
  84. C
  85.       u = ranf()
  86.       IF (d*u.LE.t*t*t) RETURN
  87. C
  88. C     STEP  4:  RECALCULATIONS OF Q0,B,SI,C IF NECESSARY
  89. C
  90.       IF (a.EQ.aaa) GO TO 40
  91.       aaa = a
  92.       r = 1.0/a
  93.       q0 = ((((((q7*r+q6)*r+q5)*r+q4)*r+q3)*r+q2)*r+q1)*r
  94. C
  95. C               APPROXIMATION DEPENDING ON SIZE OF PARAMETER A
  96. C               THE CONSTANTS IN THE EXPRESSIONS FOR B, SI AND
  97. C               C WERE ESTABLISHED BY NUMERICAL EXPERIMENTS
  98. C
  99.       IF (a.LE.3.686) GO TO 30
  100.       IF (a.LE.13.022) GO TO 20
  101. C
  102. C               CASE 3:  A .GT. 13.022
  103. C
  104.       b = 1.77
  105.       si = .75
  106.       c = .1515/s
  107.       GO TO 40
  108. C
  109. C               CASE 2:  3.686 .LT. A .LE. 13.022
  110. C
  111.    20 b = 1.654 + .0076*s2
  112.       si = 1.68/s + .275
  113.       c = .062/s + .024
  114.       GO TO 40
  115. C
  116. C               CASE 1:  A .LE. 3.686
  117. C
  118.    30 b = .463 + s - .178*s2
  119.       si = 1.235
  120.       c = .195/s - .079 + .016*s
  121. C
  122. C     STEP  5:  NO QUOTIENT TEST IF X NOT POSITIVE
  123. C
  124.    40 IF (x.LE.0.0) GO TO 70
  125. C
  126. C     STEP  6:  CALCULATION OF V AND QUOTIENT Q
  127. C
  128.       v = t/ (s+s)
  129.       IF (abs(v).LE.0.25) GO TO 50
  130.       q = q0 - s*t + 0.25*t*t + (s2+s2)*alog(1.0+v)
  131.       GO TO 60
  132.  
  133.    50 q = q0 + 0.5*t*t* ((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v
  134. C
  135. C     STEP  7:  QUOTIENT ACCEPTANCE (Q)
  136. C
  137.    60 IF (alog(1.0-u).LE.q) RETURN
  138. C
  139. C     STEP  8:  E=STANDARD EXPONENTIAL DEVIATE
  140. C               U= 0,1 -UNIFORM DEVIATE
  141. C               T=(B,SI)-DOUBLE EXPONENTIAL (LAPLACE) SAMPLE
  142. C
  143.    70 e = sexpo()
  144.       u = ranf()
  145.       u = u + u - 1.0
  146.       t = b + sign(si*e,u)
  147. C
  148. C     STEP  9:  REJECTION IF T .LT. TAU(1) = -.71874483771719
  149. C
  150.       IF (t.LT. (-.7187449)) GO TO 70
  151. C
  152. C     STEP 10:  CALCULATION OF V AND QUOTIENT Q
  153. C
  154.       v = t/ (s+s)
  155.       IF (abs(v).LE.0.25) GO TO 80
  156.       q = q0 - s*t + 0.25*t*t + (s2+s2)*alog(1.0+v)
  157.       GO TO 90
  158.  
  159.    80 q = q0 + 0.5*t*t* ((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v
  160. C
  161. C     STEP 11:  HAT ACCEPTANCE (H) (IF Q NOT POSITIVE GO TO STEP 8)
  162. C
  163.    90 IF (q.LE.0.0) GO TO 70
  164.       IF (q.LE.0.5) GO TO 100
  165.       w = exp(q) - 1.0
  166.       GO TO 110
  167.  
  168.   100 w = ((((e5*q+e4)*q+e3)*q+e2)*q+e1)*q
  169. C
  170. C               IF T IS REJECTED, SAMPLE AGAIN AT STEP 8
  171. C
  172.   110 IF (c*abs(u).GT.w*exp(e-0.5*t*t)) GO TO 70
  173.       x = s + 0.5*t
  174.       sgamma = x*x
  175.       RETURN
  176. C
  177. C     ALTERNATE METHOD FOR PARAMETERS A BELOW 1  (.3678794=EXP(-1.))
  178. C
  179.   120 aa = 0.0
  180.       b = 1.0 + .3678794*a
  181.   130 p = b*ranf()
  182.       IF (p.GE.1.0) GO TO 140
  183.       sgamma = exp(alog(p)/a)
  184.       IF (sexpo().LT.sgamma) GO TO 130
  185.       RETURN
  186.  
  187.   140 sgamma = -alog((b-p)/a)
  188.       IF (sexpo().LT. (1.0-a)*alog(sgamma)) GO TO 130
  189.       RETURN
  190.  
  191.       END
  192.