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 / ignbin.f < prev    next >
Text File  |  1996-09-28  |  9KB  |  304 lines
  1.       INTEGER FUNCTION ignbin(n,pp)
  2. C**********************************************************************
  3. C
  4. C     INTEGER FUNCTION IGNBIN( N, P )
  5. C
  6. C                    GENerate BINomial random deviate
  7. C
  8. C
  9. C                              Function
  10. C
  11. C
  12. C     Generates a single random deviate from a binomial
  13. C     distribution whose number of trials is N and whose
  14. C     probability of an event in each trial is P.
  15. C
  16. C
  17. C                              Arguments
  18. C
  19. C
  20. C     N  --> The number of trials in the binomial distribution
  21. C            from which a random deviate is to be generated.
  22. C                              INTEGER N
  23. C
  24. C     P  --> The probability of an event in each trial of the
  25. C            binomial distribution from which a random deviate
  26. C            is to be generated.
  27. C                              REAL P
  28. C
  29. C     IGNBIN <-- A random deviate yielding the number of events
  30. C                from N independent trials, each of which has
  31. C                a probability of event P.
  32. C                              INTEGER IGNBIN
  33. C
  34. C
  35. C                              Note
  36. C
  37. C
  38. C     Uses RANF so the value of the seeds, ISEED1 and ISEED2 must be set
  39. C     by a call similar to the following
  40. C          DUM = RANSET( ISEED1, ISEED2 )
  41. C
  42. C
  43. C                              Method
  44. C
  45. C
  46. C     This is algorithm BTPE from:
  47. C
  48. C         Kachitvichyanukul, V. and Schmeiser, B. W.
  49. C
  50. C         Binomial Random Variate Generation.
  51. C         Communications of the ACM, 31, 2
  52. C         (February, 1988) 216.
  53. C
  54. C**********************************************************************
  55. C     SUBROUTINE BTPEC(N,PP,ISEED,JX)
  56. C
  57. C     BINOMIAL RANDOM VARIATE GENERATOR
  58. C     MEAN .LT. 30 -- INVERSE CDF
  59. C       MEAN .GE. 30 -- ALGORITHM BTPE:  ACCEPTANCE-REJECTION VIA
  60. C       FOUR REGION COMPOSITION.  THE FOUR REGIONS ARE A TRIANGLE
  61. C       (SYMMETRIC IN THE CENTER), A PAIR OF PARALLELOGRAMS (ABOVE
  62. C       THE TRIANGLE), AND EXPONENTIAL LEFT AND RIGHT TAILS.
  63. C
  64. C     BTPE REFERS TO BINOMIAL-TRIANGLE-PARALLELOGRAM-EXPONENTIAL.
  65. C     BTPEC REFERS TO BTPE AND "COMBINED."  THUS BTPE IS THE
  66. C       RESEARCH AND BTPEC IS THE IMPLEMENTATION OF A COMPLETE
  67. C       USABLE ALGORITHM.
  68. C     REFERENCE:  VORATAS KACHITVICHYANUKUL AND BRUCE SCHMEISER,
  69. C       "BINOMIAL RANDOM VARIATE GENERATION,"
  70. C       COMMUNICATIONS OF THE ACM, FORTHCOMING
  71. C     WRITTEN:  SEPTEMBER 1980.
  72. C       LAST REVISED:  MAY 1985, JULY 1987
  73. C     REQUIRED SUBPROGRAM:  RAND() -- A UNIFORM (0,1) RANDOM NUMBER
  74. C                           GENERATOR
  75. C     ARGUMENTS
  76. C
  77. C       N : NUMBER OF BERNOULLI TRIALS            (INPUT)
  78. C       PP : PROBABILITY OF SUCCESS IN EACH TRIAL (INPUT)
  79. C       ISEED:  RANDOM NUMBER SEED                (INPUT AND OUTPUT)
  80. C       JX:  RANDOMLY GENERATED OBSERVATION       (OUTPUT)
  81. C
  82. C     VARIABLES
  83. C       PSAVE: VALUE OF PP FROM THE LAST CALL TO BTPEC
  84. C       NSAVE: VALUE OF N FROM THE LAST CALL TO BTPEC
  85. C       XNP:  VALUE OF THE MEAN FROM THE LAST CALL TO BTPEC
  86. C
  87. C       P: PROBABILITY USED IN THE GENERATION PHASE OF BTPEC
  88. C       FFM: TEMPORARY VARIABLE EQUAL TO XNP + P
  89. C       M:  INTEGER VALUE OF THE CURRENT MODE
  90. C       FM:  FLOATING POINT VALUE OF THE CURRENT MODE
  91. C       XNPQ: TEMPORARY VARIABLE USED IN SETUP AND SQUEEZING STEPS
  92. C       P1:  AREA OF THE TRIANGLE
  93. C       C:  HEIGHT OF THE PARALLELOGRAMS
  94. C       XM:  CENTER OF THE TRIANGLE
  95. C       XL:  LEFT END OF THE TRIANGLE
  96. C       XR:  RIGHT END OF THE TRIANGLE
  97. C       AL:  TEMPORARY VARIABLE
  98. C       XLL:  RATE FOR THE LEFT EXPONENTIAL TAIL
  99. C       XLR:  RATE FOR THE RIGHT EXPONENTIAL TAIL
  100. C       P2:  AREA OF THE PARALLELOGRAMS
  101. C       P3:  AREA OF THE LEFT EXPONENTIAL TAIL
  102. C       P4:  AREA OF THE RIGHT EXPONENTIAL TAIL
  103. C       U:  A U(0,P4) RANDOM VARIATE USED FIRST TO SELECT ONE OF THE
  104. C           FOUR REGIONS AND THEN CONDITIONALLY TO GENERATE A VALUE
  105. C           FROM THE REGION
  106. C       V:  A U(0,1) RANDOM NUMBER USED TO GENERATE THE RANDOM VALUE
  107. C           (REGION 1) OR TRANSFORMED INTO THE VARIATE TO ACCEPT OR
  108. C           REJECT THE CANDIDATE VALUE
  109. C       IX:  INTEGER CANDIDATE VALUE
  110. C       X:  PRELIMINARY CONTINUOUS CANDIDATE VALUE IN REGION 2 LOGIC
  111. C           AND A FLOATING POINT IX IN THE ACCEPT/REJECT LOGIC
  112. C       K:  ABSOLUTE VALUE OF (IX-M)
  113. C       F:  THE HEIGHT OF THE SCALED DENSITY FUNCTION USED IN THE
  114. C           ACCEPT/REJECT DECISION WHEN BOTH M AND IX ARE SMALL
  115. C           ALSO USED IN THE INVERSE TRANSFORMATION
  116. C       R: THE RATIO P/Q
  117. C       G: CONSTANT USED IN CALCULATION OF PROBABILITY
  118. C       MP:  MODE PLUS ONE, THE LOWER INDEX FOR EXPLICIT CALCULATION
  119. C            OF F WHEN IX IS GREATER THAN M
  120. C       IX1:  CANDIDATE VALUE PLUS ONE, THE LOWER INDEX FOR EXPLICIT
  121. C             CALCULATION OF F WHEN IX IS LESS THAN M
  122. C       I:  INDEX FOR EXPLICIT CALCULATION OF F FOR BTPE
  123. C       AMAXP: MAXIMUM ERROR OF THE LOGARITHM OF NORMAL BOUND
  124. C       YNORM: LOGARITHM OF NORMAL BOUND
  125. C       ALV:  NATURAL LOGARITHM OF THE ACCEPT/REJECT VARIATE V
  126. C
  127. C       X1,F1,Z,W,Z2,X2,F2, AND W2 ARE TEMPORARY VARIABLES TO BE
  128. C       USED IN THE FINAL ACCEPT/REJECT TEST
  129. C
  130. C       QN: PROBABILITY OF NO SUCCESS IN N TRIALS
  131. C
  132. C     REMARK
  133. C       IX AND JX COULD LOGICALLY BE THE SAME VARIABLE, WHICH WOULD
  134. C       SAVE A MEMORY POSITION AND A LINE OF CODE.  HOWEVER, SOME
  135. C       COMPILERS (E.G.,CDC MNF) OPTIMIZE BETTER WHEN THE ARGUMENTS
  136. C       ARE NOT INVOLVED.
  137. C
  138. C     ISEED NEEDS TO BE DOUBLE PRECISION IF THE IMSL ROUTINE
  139. C     GGUBFS IS USED TO GENERATE UNIFORM RANDOM NUMBER, OTHERWISE
  140. C     TYPE OF ISEED SHOULD BE DICTATED BY THE UNIFORM GENERATOR
  141. C
  142. C**********************************************************************
  143.  
  144. C
  145. C
  146. C
  147. C*****DETERMINE APPROPRIATE ALGORITHM AND WHETHER SETUP IS NECESSARY
  148. C
  149. C     ..
  150. C     .. Scalar Arguments ..
  151.       REAL pp
  152.       INTEGER n
  153. C     ..
  154. C     .. Local Scalars ..
  155.       REAL al,alv,amaxp,c,f,f1,f2,ffm,fm,g,p,p1,p2,p3,p4,psave,q,qn,r,u,
  156.      +     v,w,w2,x,x1,x2,xl,xll,xlr,xm,xnp,xnpq,xr,ynorm,z,z2
  157.       INTEGER i,ix,ix1,k,m,mp,nsave
  158. C     ..
  159. C     .. External Functions ..
  160.       REAL ranf
  161.       EXTERNAL ranf
  162. C     ..
  163. C     .. Intrinsic Functions ..
  164.       INTRINSIC abs,alog,amin1,iabs,int,sqrt
  165. C     ..
  166. C     .. Data statements ..
  167.       DATA psave,nsave/-1.,-1/
  168. C     ..
  169. C     .. Executable Statements ..
  170.       IF (pp.NE.psave) GO TO 10
  171.       IF (n.NE.nsave) GO TO 20
  172.       IF (xnp-30.) 150,30,30
  173. C
  174. C*****SETUP, PERFORM ONLY WHEN PARAMETERS CHANGE
  175. C
  176.    10 psave = pp
  177.       p = amin1(psave,1.-psave)
  178.       q = 1. - p
  179.    20 xnp = n*p
  180.       nsave = n
  181.       IF (xnp.LT.30.) GO TO 140
  182.       ffm = xnp + p
  183.       m = ffm
  184.       fm = m
  185.       xnpq = xnp*q
  186.       p1 = int(2.195*sqrt(xnpq)-4.6*q) + 0.5
  187.       xm = fm + 0.5
  188.       xl = xm - p1
  189.       xr = xm + p1
  190.       c = 0.134 + 20.5/ (15.3+fm)
  191.       al = (ffm-xl)/ (ffm-xl*p)
  192.       xll = al* (1.+.5*al)
  193.       al = (xr-ffm)/ (xr*q)
  194.       xlr = al* (1.+.5*al)
  195.       p2 = p1* (1.+c+c)
  196.       p3 = p2 + c/xll
  197.       p4 = p3 + c/xlr
  198. C      WRITE(6,100) N,P,P1,P2,P3,P4,XL,XR,XM,FM
  199. C  100 FORMAT(I15,4F18.7/5F18.7)
  200. C
  201. C*****GENERATE VARIATE
  202. C
  203.    30 u = ranf()*p4
  204.       v = ranf()
  205. C
  206. C     TRIANGULAR REGION
  207. C
  208.       IF (u.GT.p1) GO TO 40
  209.       ix = xm - p1*v + u
  210.       GO TO 170
  211. C
  212. C     PARALLELOGRAM REGION
  213. C
  214.    40 IF (u.GT.p2) GO TO 50
  215.       x = xl + (u-p1)/c
  216.       v = v*c + 1. - abs(xm-x)/p1
  217.       IF (v.GT.1. .OR. v.LE.0.) GO TO 30
  218.       ix = x
  219.       GO TO 70
  220. C
  221. C     LEFT TAIL
  222. C
  223.    50 IF (u.GT.p3) GO TO 60
  224.       ix = xl + alog(v)/xll
  225.       IF (ix.LT.0) GO TO 30
  226.       v = v* (u-p2)*xll
  227.       GO TO 70
  228. C
  229. C     RIGHT TAIL
  230. C
  231.    60 ix = xr - alog(v)/xlr
  232.       IF (ix.GT.n) GO TO 30
  233.       v = v* (u-p3)*xlr
  234. C
  235. C*****DETERMINE APPROPRIATE WAY TO PERFORM ACCEPT/REJECT TEST
  236. C
  237.    70 k = iabs(ix-m)
  238.       IF (k.GT.20 .AND. k.LT.xnpq/2-1) GO TO 130
  239. C
  240. C     EXPLICIT EVALUATION
  241. C
  242.       f = 1.0
  243.       r = p/q
  244.       g = (n+1)*r
  245.       IF (m-ix) 80,120,100
  246.    80 mp = m + 1
  247.       DO 90 i = mp,ix
  248.           f = f* (g/i-r)
  249.    90 CONTINUE
  250.       GO TO 120
  251.  
  252.   100 ix1 = ix + 1
  253.       DO 110 i = ix1,m
  254.           f = f/ (g/i-r)
  255.   110 CONTINUE
  256.   120 IF (v-f) 170,170,30
  257. C
  258. C     SQUEEZING USING UPPER AND LOWER BOUNDS ON ALOG(F(X))
  259. C
  260.   130 amaxp = (k/xnpq)* ((k* (k/3.+.625)+.1666666666666)/xnpq+.5)
  261.       ynorm = -k*k/ (2.*xnpq)
  262.       alv = alog(v)
  263.       IF (alv.LT.ynorm-amaxp) GO TO 170
  264.       IF (alv.GT.ynorm+amaxp) GO TO 30
  265. C
  266. C     STIRLING'S FORMULA TO MACHINE ACCURACY FOR
  267. C     THE FINAL ACCEPTANCE/REJECTION TEST
  268. C
  269.       x1 = ix + 1
  270.       f1 = fm + 1.
  271.       z = n + 1 - fm
  272.       w = n - ix + 1.
  273.       z2 = z*z
  274.       x2 = x1*x1
  275.       f2 = f1*f1
  276.       w2 = w*w
  277.       IF (alv- (xm*alog(f1/x1)+ (n-m+.5)*alog(z/w)+ (ix-
  278.      +    m)*alog(w*p/ (x1*q))+ (13860.- (462.- (132.- (99.-
  279.      +    140./f2)/f2)/f2)/f2)/f1/166320.+ (13860.- (462.- (132.- (99.-
  280.      +    140./z2)/z2)/z2)/z2)/z/166320.+ (13860.- (462.- (132.- (99.-
  281.      +    140./x2)/x2)/x2)/x2)/x1/166320.+ (13860.- (462.- (132.- (99.-
  282.      +    140./w2)/w2)/w2)/w2)/w/166320.)) 170,170,30
  283. C
  284. C     INVERSE CDF LOGIC FOR MEAN LESS THAN 30
  285. C
  286.   140 qn = q**n
  287.       r = p/q
  288.       g = r* (n+1)
  289.   150 ix = 0
  290.       f = qn
  291.       u = ranf()
  292.   160 IF (u.LT.f) GO TO 170
  293.       IF (ix.GT.110) GO TO 150
  294.       u = u - f
  295.       ix = ix + 1
  296.       f = f* (g/ix-r)
  297.       GO TO 160
  298.  
  299.   170 IF (psave.GT.0.5) ix = n - ix
  300.       ignbin = ix
  301.       RETURN
  302.  
  303.       END
  304.