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 / ignpoi.f < prev    next >
Text File  |  1996-09-28  |  8KB  |  262 lines
  1.       INTEGER FUNCTION ignpoi(mu)
  2. C**********************************************************************
  3. C
  4. C     INTEGER FUNCTION IGNPOI( AV )
  5. C
  6. C                    GENerate POIsson random deviate
  7. C
  8. C
  9. C                              Function
  10. C
  11. C
  12. C     Generates a single random deviate from a Poisson
  13. C     distribution with mean AV.
  14. C
  15. C
  16. C                              Arguments
  17. C
  18. C
  19. C     AV --> The mean of the Poisson distribution from which
  20. C            a random deviate is to be generated.
  21. C                              REAL AV
  22. C
  23. C     GENEXP <-- The random deviate.
  24. C                              REAL GENEXP
  25. C
  26. C
  27. C                              Method
  28. C
  29. C
  30. C     Renames KPOIS from TOMS as slightly modified by BWB to use RANF
  31. C     instead of SUNIF.
  32. C
  33. C     For details see:
  34. C
  35. C               Ahrens, J.H. and Dieter, U.
  36. C               Computer Generation of Poisson Deviates
  37. C               From Modified Normal Distributions.
  38. C               ACM Trans. Math. Software, 8, 2
  39. C               (June 1982),163-179
  40. C
  41. C**********************************************************************
  42. C**********************************************************************C
  43. C**********************************************************************C
  44. C                                                                      C
  45. C                                                                      C
  46. C     P O I S S O N  DISTRIBUTION                                      C
  47. C                                                                      C
  48. C                                                                      C
  49. C**********************************************************************C
  50. C**********************************************************************C
  51. C                                                                      C
  52. C     FOR DETAILS SEE:                                                 C
  53. C                                                                      C
  54. C               AHRENS, J.H. AND DIETER, U.                            C
  55. C               COMPUTER GENERATION OF POISSON DEVIATES                C
  56. C               FROM MODIFIED NORMAL DISTRIBUTIONS.                    C
  57. C               ACM TRANS. MATH. SOFTWARE, 8,2 (JUNE 1982), 163 - 179. C
  58. C                                                                      C
  59. C     (SLIGHTLY MODIFIED VERSION OF THE PROGRAM IN THE ABOVE ARTICLE)  C
  60. C                                                                      C
  61. C**********************************************************************C
  62. C
  63. C      INTEGER FUNCTION IGNPOI(IR,MU)
  64. C
  65. C     INPUT:  IR=CURRENT STATE OF BASIC RANDOM NUMBER GENERATOR
  66. C             MU=MEAN MU OF THE POISSON DISTRIBUTION
  67. C     OUTPUT: IGNPOI=SAMPLE FROM THE POISSON-(MU)-DISTRIBUTION
  68. C
  69. C
  70. C
  71. C     MUPREV=PREVIOUS MU, MUOLD=MU AT LAST EXECUTION OF STEP P OR B.
  72. C     TABLES: COEFFICIENTS A0-A7 FOR STEP F. FACTORIALS FACT
  73. C     COEFFICIENTS A(K) - FOR PX = FK*V*V*SUM(A(K)*V**K)-DEL
  74. C
  75. C
  76. C
  77. C     SEPARATION OF CASES A AND B
  78. C
  79. C     .. Scalar Arguments ..
  80.       REAL mu
  81. C     ..
  82. C     .. Local Scalars ..
  83.       REAL a0,a1,a2,a3,a4,a5,a6,a7,b1,b2,c,c0,c1,c2,c3,d,del,difmuk,e,
  84.      +     fk,fx,fy,g,muold,muprev,omega,p,p0,px,py,q,s,t,u,v,x,xx
  85.       INTEGER j,k,kflag,l,m
  86. C     ..
  87. C     .. Local Arrays ..
  88.       REAL fact(10),pp(35)
  89. C     ..
  90. C     .. External Functions ..
  91.       REAL ranf,sexpo,snorm
  92.       EXTERNAL ranf,sexpo,snorm
  93. C     ..
  94. C     .. Intrinsic Functions ..
  95.       INTRINSIC abs,alog,exp,float,ifix,max0,min0,sign,sqrt
  96. C     ..
  97. C     .. Data statements ..
  98.       DATA muprev,muold/0.,0./
  99.       DATA a0,a1,a2,a3,a4,a5,a6,a7/-.5,.3333333,-.2500068,.2000118,
  100.      +     -.1661269,.1421878,-.1384794,.1250060/
  101.       DATA fact/1.,1.,2.,6.,24.,120.,720.,5040.,40320.,362880./
  102. C     ..
  103. C     .. Executable Statements ..
  104.       IF (mu.EQ.muprev) GO TO 10
  105.       IF (mu.LT.10.0) GO TO 120
  106. C
  107. C     C A S E  A. (RECALCULATION OF S,D,L IF MU HAS CHANGED)
  108. C
  109.       muprev = mu
  110.       s = sqrt(mu)
  111.       d = 6.0*mu*mu
  112. C
  113. C             THE POISSON PROBABILITIES PK EXCEED THE DISCRETE NORMAL
  114. C             PROBABILITIES FK WHENEVER K >= M(MU). L=IFIX(MU-1.1484)
  115. C             IS AN UPPER BOUND TO M(MU) FOR ALL MU >= 10 .
  116. C
  117.       l = ifix(mu-1.1484)
  118. C
  119. C     STEP N. NORMAL SAMPLE - SNORM(IR) FOR STANDARD NORMAL DEVIATE
  120. C
  121.    10 g = mu + s*snorm()
  122.       IF (g.LT.0.0) GO TO 20
  123.       ignpoi = ifix(g)
  124. C
  125. C     STEP I. IMMEDIATE ACCEPTANCE IF IGNPOI IS LARGE ENOUGH
  126. C
  127.       IF (ignpoi.GE.l) RETURN
  128. C
  129. C     STEP S. SQUEEZE ACCEPTANCE - SUNIF(IR) FOR (0,1)-SAMPLE U
  130. C
  131.       fk = float(ignpoi)
  132.       difmuk = mu - fk
  133.       u = ranf()
  134.       IF (d*u.GE.difmuk*difmuk*difmuk) RETURN
  135. C
  136. C     STEP P. PREPARATIONS FOR STEPS Q AND H.
  137. C             (RECALCULATIONS OF PARAMETERS IF NECESSARY)
  138. C             .3989423=(2*PI)**(-.5)  .416667E-1=1./24.  .1428571=1./7.
  139. C             THE QUANTITIES B1, B2, C3, C2, C1, C0 ARE FOR THE HERMITE
  140. C             APPROXIMATIONS TO THE DISCRETE NORMAL PROBABILITIES FK.
  141. C             C=.1069/MU GUARANTEES MAJORIZATION BY THE 'HAT'-FUNCTION.
  142. C
  143.    20 IF (mu.EQ.muold) GO TO 30
  144.       muold = mu
  145.       omega = .3989423/s
  146.       b1 = .4166667E-1/mu
  147.       b2 = .3*b1*b1
  148.       c3 = .1428571*b1*b2
  149.       c2 = b2 - 15.*c3
  150.       c1 = b1 - 6.*b2 + 45.*c3
  151.       c0 = 1. - b1 + 3.*b2 - 15.*c3
  152.       c = .1069/mu
  153.    30 IF (g.LT.0.0) GO TO 50
  154. C
  155. C             'SUBROUTINE' F IS CALLED (KFLAG=0 FOR CORRECT RETURN)
  156. C
  157.       kflag = 0
  158.       GO TO 70
  159. C
  160. C     STEP Q. QUOTIENT ACCEPTANCE (RARE CASE)
  161. C
  162.    40 IF (fy-u*fy.LE.py*exp(px-fx)) RETURN
  163. C
  164. C     STEP E. EXPONENTIAL SAMPLE - SEXPO(IR) FOR STANDARD EXPONENTIAL
  165. C             DEVIATE E AND SAMPLE T FROM THE LAPLACE 'HAT'
  166. C             (IF T <= -.6744 THEN PK < FK FOR ALL MU >= 10.)
  167. C
  168.    50 e = sexpo()
  169.       u = ranf()
  170.       u = u + u - 1.0
  171.       t = 1.8 + sign(e,u)
  172.       IF (t.LE. (-.6744)) GO TO 50
  173.       ignpoi = ifix(mu+s*t)
  174.       fk = float(ignpoi)
  175.       difmuk = mu - fk
  176. C
  177. C             'SUBROUTINE' F IS CALLED (KFLAG=1 FOR CORRECT RETURN)
  178. C
  179.       kflag = 1
  180.       GO TO 70
  181. C
  182. C     STEP H. HAT ACCEPTANCE (E IS REPEATED ON REJECTION)
  183. C
  184.    60 IF (c*abs(u).GT.py*exp(px+e)-fy*exp(fx+e)) GO TO 50
  185.       RETURN
  186. C
  187. C     STEP F. 'SUBROUTINE' F. CALCULATION OF PX,PY,FX,FY.
  188. C             CASE IGNPOI .LT. 10 USES FACTORIALS FROM TABLE FACT
  189. C
  190.    70 IF (ignpoi.GE.10) GO TO 80
  191.       px = -mu
  192.       py = mu**ignpoi/fact(ignpoi+1)
  193.       GO TO 110
  194. C
  195. C             CASE IGNPOI .GE. 10 USES POLYNOMIAL APPROXIMATION
  196. C             A0-A7 FOR ACCURACY WHEN ADVISABLE
  197. C             .8333333E-1=1./12.  .3989423=(2*PI)**(-.5)
  198. C
  199.    80 del = .8333333E-1/fk
  200.       del = del - 4.8*del*del*del
  201.       v = difmuk/fk
  202.       IF (abs(v).LE.0.25) GO TO 90
  203.       px = fk*alog(1.0+v) - difmuk - del
  204.       GO TO 100
  205.  
  206.    90 px = fk*v*v* (((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v+a0) -
  207.      +     del
  208.   100 py = .3989423/sqrt(fk)
  209.   110 x = (0.5-difmuk)/s
  210.       xx = x*x
  211.       fx = -0.5*xx
  212.       fy = omega* (((c3*xx+c2)*xx+c1)*xx+c0)
  213.       IF (kflag) 40,40,60
  214. C
  215. C     C A S E  B. (START NEW TABLE AND CALCULATE P0 IF NECESSARY)
  216. C
  217.   120 muprev = 0.0
  218.       IF (mu.EQ.muold) GO TO 130
  219.       muold = mu
  220.       m = max0(1,ifix(mu))
  221.       l = 0
  222.       p = exp(-mu)
  223.       q = p
  224.       p0 = p
  225. C
  226. C     STEP U. UNIFORM SAMPLE FOR INVERSION METHOD
  227. C
  228.   130 u = ranf()
  229.       ignpoi = 0
  230.       IF (u.LE.p0) RETURN
  231. C
  232. C     STEP T. TABLE COMPARISON UNTIL THE END PP(L) OF THE
  233. C             PP-TABLE OF CUMULATIVE POISSON PROBABILITIES
  234. C             (0.458=PP(9) FOR MU=10)
  235. C
  236.       IF (l.EQ.0) GO TO 150
  237.       j = 1
  238.       IF (u.GT.0.458) j = min0(l,m)
  239.       DO 140 k = j,l
  240.           IF (u.LE.pp(k)) GO TO 180
  241.   140 CONTINUE
  242.       IF (l.EQ.35) GO TO 130
  243. C
  244. C     STEP C. CREATION OF NEW POISSON PROBABILITIES P
  245. C             AND THEIR CUMULATIVES Q=PP(K)
  246. C
  247.   150 l = l + 1
  248.       DO 160 k = l,35
  249.           p = p*mu/float(k)
  250.           q = q + p
  251.           pp(k) = q
  252.           IF (u.LE.q) GO TO 170
  253.   160 CONTINUE
  254.       l = 35
  255.       GO TO 130
  256.  
  257.   170 l = k
  258.   180 ignpoi = k
  259.       RETURN
  260.  
  261.       END
  262.