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 / genbet.f < prev    next >
Text File  |  1996-09-28  |  4KB  |  193 lines
  1.       REAL FUNCTION genbet(aa,bb)
  2. C**********************************************************************
  3. C
  4. C     REAL FUNCTION GENBET( A, B )
  5. C               GeNerate BETa random deviate
  6. C
  7. C
  8. C                              Function
  9. C
  10. C
  11. C     Returns a single random deviate from the beta distribution with
  12. C     parameters A and B.  The density of the beta is
  13. C               x^(a-1) * (1-x)^(b-1) / B(a,b) for 0 < x < 1
  14. C
  15. C
  16. C                              Arguments
  17. C
  18. C
  19. C     A --> First parameter of the beta distribution
  20. C                         REAL A
  21. C
  22. C     B --> Second parameter of the beta distribution
  23. C                         REAL B
  24. C
  25. C
  26. C                              Method
  27. C
  28. C
  29. C     R. C. H. Cheng
  30. C     Generating Beta Variatew with Nonintegral Shape Parameters
  31. C     Communications of the ACM, 21:317-322  (1978)
  32. C     (Algorithms BB and BC)
  33. C
  34. C**********************************************************************
  35. C     .. Parameters ..
  36. C     Close to the largest number that can be exponentiated
  37.       REAL expmax
  38.       PARAMETER (expmax=89.0)
  39. C     Close to the largest representable single precision number
  40.       REAL infnty
  41.       PARAMETER (infnty=1.0E38)
  42. C     ..
  43. C     .. Scalar Arguments ..
  44.       REAL aa,bb
  45. C     ..
  46. C     .. Local Scalars ..
  47.       REAL a,alpha,b,beta,delta,gamma,k1,k2,olda,oldb,r,s,t,u1,u2,v,w,y,
  48.      +     z
  49.       LOGICAL qsame
  50. C     ..
  51. C     .. External Functions ..
  52.       REAL ranf
  53.       EXTERNAL ranf
  54. C     ..
  55. C     .. Intrinsic Functions ..
  56.       INTRINSIC exp,log,max,min,sqrt
  57. C     ..
  58. C     .. Save statement ..
  59.       SAVE olda,oldb,alpha,beta,gamma,k1,k2
  60. C     ..
  61. C     .. Data statements ..
  62.       DATA olda,oldb/-1,-1/
  63. C     ..
  64. C     .. Executable Statements ..
  65.       qsame = (olda.EQ.aa) .AND. (oldb.EQ.bb)
  66.       IF (qsame) GO TO 20
  67.       IF (.NOT. (aa.LE.0.0.OR.bb.LE.0.0)) GO TO 10
  68.       WRITE (*,*) ' AA or BB <= 0 in GENBET - Abort!'
  69.       WRITE (*,*) ' AA: ',aa,' BB ',bb
  70.       CALL XSTOPX (' AA or BB <= 0 in GENBET - Abort!')
  71.  
  72.    10 olda = aa
  73.       oldb = bb
  74.    20 IF (.NOT. (min(aa,bb).GT.1.0)) GO TO 100
  75.  
  76.  
  77. C     Alborithm BB
  78.  
  79. C
  80. C     Initialize
  81. C
  82.       IF (qsame) GO TO 30
  83.       a = min(aa,bb)
  84.       b = max(aa,bb)
  85.       alpha = a + b
  86.       beta = sqrt((alpha-2.0)/ (2.0*a*b-alpha))
  87.       gamma = a + 1.0/beta
  88.    30 CONTINUE
  89.    40 u1 = ranf()
  90. C
  91. C     Step 1
  92. C
  93.       u2 = ranf()
  94.       v = beta*log(u1/ (1.0-u1))
  95.       IF (.NOT. (v.GT.expmax)) GO TO 50
  96.       w = infnty
  97.       GO TO 60
  98.  
  99.    50 w = a*exp(v)
  100.    60 z = u1**2*u2
  101.       r = gamma*v - 1.3862944
  102.       s = a + r - w
  103. C
  104. C     Step 2
  105. C
  106.       IF ((s+2.609438).GE. (5.0*z)) GO TO 70
  107. C
  108. C     Step 3
  109. C
  110.       t = log(z)
  111.       IF (s.GT.t) GO TO 70
  112. C
  113. C     Step 4
  114. C
  115.       IF ((r+alpha*log(alpha/ (b+w))).LT.t) GO TO 40
  116. C
  117. C     Step 5
  118. C
  119.    70 IF (.NOT. (aa.EQ.a)) GO TO 80
  120.       genbet = w/ (b+w)
  121.       GO TO 90
  122.  
  123.    80 genbet = b/ (b+w)
  124.    90 GO TO 230
  125.  
  126.  
  127. C     Algorithm BC
  128.  
  129. C
  130. C     Initialize
  131. C
  132.   100 IF (qsame) GO TO 110
  133.       a = max(aa,bb)
  134.       b = min(aa,bb)
  135.       alpha = a + b
  136.       beta = 1.0/b
  137.       delta = 1.0 + a - b
  138.       k1 = delta* (0.0138889+0.0416667*b)/ (a*beta-0.777778)
  139.       k2 = 0.25 + (0.5+0.25/delta)*b
  140.   110 CONTINUE
  141.   120 u1 = ranf()
  142. C
  143. C     Step 1
  144. C
  145.       u2 = ranf()
  146.       IF (u1.GE.0.5) GO TO 130
  147. C
  148. C     Step 2
  149. C
  150.       y = u1*u2
  151.       z = u1*y
  152.       IF ((0.25*u2+z-y).GE.k1) GO TO 120
  153.       GO TO 170
  154. C
  155. C     Step 3
  156. C
  157.   130 z = u1**2*u2
  158.       IF (.NOT. (z.LE.0.25)) GO TO 160
  159.       v = beta*log(u1/ (1.0-u1))
  160.       IF (.NOT. (v.GT.expmax)) GO TO 140
  161.       w = infnty
  162.       GO TO 150
  163.  
  164.   140 w = a*exp(v)
  165.   150 GO TO 200
  166.  
  167.   160 IF (z.GE.k2) GO TO 120
  168. C
  169. C     Step 4
  170. C
  171. C
  172. C     Step 5
  173. C
  174.   170 v = beta*log(u1/ (1.0-u1))
  175.       IF (.NOT. (v.GT.expmax)) GO TO 180
  176.       w = infnty
  177.       GO TO 190
  178.  
  179.   180 w = a*exp(v)
  180.   190 IF ((alpha* (log(alpha/ (b+w))+v)-1.3862944).LT.log(z)) GO TO 120
  181. C
  182. C     Step 6
  183. C
  184.   200 IF (.NOT. (a.EQ.aa)) GO TO 210
  185.       genbet = w/ (b+w)
  186.       GO TO 220
  187.  
  188.   210 genbet = b/ (b+w)
  189.   220 CONTINUE
  190.   230 RETURN
  191.  
  192.       END
  193.