00001 INTEGER FUNCTION ignpoi(mu) 00002 C********************************************************************** 00003 C 00004 C INTEGER FUNCTION IGNPOI( MU ) 00005 C 00006 C GENerate POIsson random deviate 00007 C 00008 C 00009 C Function 00010 C 00011 C 00012 C Generates a single random deviate from a Poisson 00013 C distribution with mean MU. 00014 C 00015 C 00016 C Arguments 00017 C 00018 C 00019 C MU --> The mean of the Poisson distribution from which 00020 C a random deviate is to be generated. 00021 C REAL MU 00022 C JJV (MU >= 0.0) 00023 C 00024 C IGNPOI <-- The random deviate. 00025 C INTEGER IGNPOI (non-negative) 00026 C 00027 C 00028 C Method 00029 C 00030 C 00031 C Renames KPOIS from TOMS as slightly modified by BWB to use RANF 00032 C instead of SUNIF. 00033 C 00034 C For details see: 00035 C 00036 C Ahrens, J.H. and Dieter, U. 00037 C Computer Generation of Poisson Deviates 00038 C From Modified Normal Distributions. 00039 C ACM Trans. Math. Software, 8, 2 00040 C (June 1982),163-179 00041 C 00042 C********************************************************************** 00043 C**********************************************************************C 00044 C**********************************************************************C 00045 C C 00046 C C 00047 C P O I S S O N DISTRIBUTION C 00048 C C 00049 C C 00050 C**********************************************************************C 00051 C**********************************************************************C 00052 C C 00053 C FOR DETAILS SEE: C 00054 C C 00055 C AHRENS, J.H. AND DIETER, U. C 00056 C COMPUTER GENERATION OF POISSON DEVIATES C 00057 C FROM MODIFIED NORMAL DISTRIBUTIONS. C 00058 C ACM TRANS. MATH. SOFTWARE, 8,2 (JUNE 1982), 163 - 179. C 00059 C C 00060 C (SLIGHTLY MODIFIED VERSION OF THE PROGRAM IN THE ABOVE ARTICLE) C 00061 C C 00062 C**********************************************************************C 00063 C 00064 C INTEGER FUNCTION IGNPOI(IR,MU) 00065 C 00066 C INPUT: IR=CURRENT STATE OF BASIC RANDOM NUMBER GENERATOR 00067 C MU=MEAN MU OF THE POISSON DISTRIBUTION 00068 C OUTPUT: IGNPOI=SAMPLE FROM THE POISSON-(MU)-DISTRIBUTION 00069 C 00070 C 00071 C 00072 C MUPREV=PREVIOUS MU, MUOLD=MU AT LAST EXECUTION OF STEP P OR CASE B 00073 C TABLES: COEFFICIENTS A0-A7 FOR STEP F. FACTORIALS FACT 00074 C COEFFICIENTS A(K) - FOR PX = FK*V*V*SUM(A(K)*V**K)-DEL 00075 C 00076 C 00077 C 00078 C SEPARATION OF CASES A AND B 00079 C 00080 C .. Scalar Arguments .. 00081 REAL mu 00082 C .. 00083 C .. Local Scalars .. 00084 REAL a0,a1,a2,a3,a4,a5,a6,a7,b1,b2,c,c0,c1,c2,c3,d,del,difmuk,e, 00085 + fk,fx,fy,g,muold,muprev,omega,p,p0,px,py,q,s,t,u,v,x,xx 00086 C JJV I added a variable 'll' here - it is the 'l' for CASE A 00087 INTEGER j,k,kflag,l,ll,m 00088 C .. 00089 C .. Local Arrays .. 00090 REAL fact(10),pp(35) 00091 C .. 00092 C .. External Functions .. 00093 REAL ranf,sexpo,snorm 00094 EXTERNAL ranf,sexpo,snorm 00095 C .. 00096 C .. Intrinsic Functions .. 00097 INTRINSIC abs,alog,exp,float,ifix,max0,min0,sign,sqrt 00098 C .. 00099 C JJV added this for case: mu unchanged 00100 C .. Save statement .. 00101 SAVE s, d, l, ll, omega, c3, c2, c1, c0, c, m, p, q, p0, 00102 + a0, a1, a2, a3, a4, a5, a6, a7, fact, muprev, muold 00103 C .. 00104 C JJV end addition - I am including vars in Data statements 00105 C .. Data statements .. 00106 C JJV changed initial values of MUPREV and MUOLD to -1.0E37 00107 C JJV if no one calls IGNPOI with MU = -1.0E37 the first time, 00108 C JJV the code shouldn't break 00109 DATA muprev,muold/-1.0E37,-1.0E37/ 00110 DATA a0,a1,a2,a3,a4,a5,a6,a7/-.5,.3333333,-.2500068,.2000118, 00111 + -.1661269,.1421878,-.1384794,.1250060/ 00112 DATA fact/1.,1.,2.,6.,24.,120.,720.,5040.,40320.,362880./ 00113 C .. 00114 C .. Executable Statements .. 00115 00116 IF (mu.EQ.muprev) GO TO 10 00117 IF (mu.LT.10.0) GO TO 120 00118 C 00119 C C A S E A. (RECALCULATION OF S,D,LL IF MU HAS CHANGED) 00120 C 00121 C JJV This is the case where I changed 'l' to 'll' 00122 C JJV Here 'll' is set once and used in a comparison once 00123 00124 muprev = mu 00125 s = sqrt(mu) 00126 d = 6.0*mu*mu 00127 C 00128 C THE POISSON PROBABILITIES PK EXCEED THE DISCRETE NORMAL 00129 C PROBABILITIES FK WHENEVER K >= M(MU). LL=IFIX(MU-1.1484) 00130 C IS AN UPPER BOUND TO M(MU) FOR ALL MU >= 10 . 00131 C 00132 ll = ifix(mu-1.1484) 00133 C 00134 C STEP N. NORMAL SAMPLE - SNORM(IR) FOR STANDARD NORMAL DEVIATE 00135 C 00136 10 g = mu + s*snorm() 00137 IF (g.LT.0.0) GO TO 20 00138 ignpoi = ifix(g) 00139 C 00140 C STEP I. IMMEDIATE ACCEPTANCE IF IGNPOI IS LARGE ENOUGH 00141 C 00142 IF (ignpoi.GE.ll) RETURN 00143 C 00144 C STEP S. SQUEEZE ACCEPTANCE - SUNIF(IR) FOR (0,1)-SAMPLE U 00145 C 00146 fk = float(ignpoi) 00147 difmuk = mu - fk 00148 u = ranf() 00149 IF (d*u.GE.difmuk*difmuk*difmuk) RETURN 00150 C 00151 C STEP P. PREPARATIONS FOR STEPS Q AND H. 00152 C (RECALCULATIONS OF PARAMETERS IF NECESSARY) 00153 C .3989423=(2*PI)**(-.5) .416667E-1=1./24. .1428571=1./7. 00154 C THE QUANTITIES B1, B2, C3, C2, C1, C0 ARE FOR THE HERMITE 00155 C APPROXIMATIONS TO THE DISCRETE NORMAL PROBABILITIES FK. 00156 C C=.1069/MU GUARANTEES MAJORIZATION BY THE 'HAT'-FUNCTION. 00157 C 00158 20 IF (mu.EQ.muold) GO TO 30 00159 muold = mu 00160 omega = .3989423/s 00161 b1 = .4166667E-1/mu 00162 b2 = .3*b1*b1 00163 c3 = .1428571*b1*b2 00164 c2 = b2 - 15.*c3 00165 c1 = b1 - 6.*b2 + 45.*c3 00166 c0 = 1. - b1 + 3.*b2 - 15.*c3 00167 c = .1069/mu 00168 30 IF (g.LT.0.0) GO TO 50 00169 C 00170 C 'SUBROUTINE' F IS CALLED (KFLAG=0 FOR CORRECT RETURN) 00171 C 00172 kflag = 0 00173 GO TO 70 00174 C 00175 C STEP Q. QUOTIENT ACCEPTANCE (RARE CASE) 00176 C 00177 40 IF (fy-u*fy.LE.py*exp(px-fx)) RETURN 00178 C 00179 C STEP E. EXPONENTIAL SAMPLE - SEXPO(IR) FOR STANDARD EXPONENTIAL 00180 C DEVIATE E AND SAMPLE T FROM THE LAPLACE 'HAT' 00181 C (IF T <= -.6744 THEN PK < FK FOR ALL MU >= 10.) 00182 C 00183 50 e = sexpo() 00184 u = ranf() 00185 u = u + u - 1.0 00186 t = 1.8 + sign(e,u) 00187 IF (t.LE. (-.6744)) GO TO 50 00188 ignpoi = ifix(mu+s*t) 00189 fk = float(ignpoi) 00190 difmuk = mu - fk 00191 C 00192 C 'SUBROUTINE' F IS CALLED (KFLAG=1 FOR CORRECT RETURN) 00193 C 00194 kflag = 1 00195 GO TO 70 00196 C 00197 C STEP H. HAT ACCEPTANCE (E IS REPEATED ON REJECTION) 00198 C 00199 60 IF (c*abs(u).GT.py*exp(px+e)-fy*exp(fx+e)) GO TO 50 00200 RETURN 00201 C 00202 C STEP F. 'SUBROUTINE' F. CALCULATION OF PX,PY,FX,FY. 00203 C CASE IGNPOI .LT. 10 USES FACTORIALS FROM TABLE FACT 00204 C 00205 70 IF (ignpoi.GE.10) GO TO 80 00206 px = -mu 00207 py = mu**ignpoi/fact(ignpoi+1) 00208 GO TO 110 00209 C 00210 C CASE IGNPOI .GE. 10 USES POLYNOMIAL APPROXIMATION 00211 C A0-A7 FOR ACCURACY WHEN ADVISABLE 00212 C .8333333E-1=1./12. .3989423=(2*PI)**(-.5) 00213 C 00214 80 del = .8333333E-1/fk 00215 del = del - 4.8*del*del*del 00216 v = difmuk/fk 00217 IF (abs(v).LE.0.25) GO TO 90 00218 px = fk*alog(1.0+v) - difmuk - del 00219 GO TO 100 00220 00221 90 px = fk*v*v* (((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v+a0) - 00222 + del 00223 100 py = .3989423/sqrt(fk) 00224 110 x = (0.5-difmuk)/s 00225 xx = x*x 00226 fx = -0.5*xx 00227 fy = omega* (((c3*xx+c2)*xx+c1)*xx+c0) 00228 IF (kflag.LE.0) GO TO 40 00229 GO TO 60 00230 C 00231 C C A S E B. (START NEW TABLE AND CALCULATE P0 IF NECESSARY) 00232 C 00233 C JJV changed MUPREV assignment from 0.0 to initial value 00234 120 muprev = -1.0E37 00235 IF (mu.EQ.muold) GO TO 130 00236 C JJV added argument checker here 00237 IF (mu.GE.0.0) GO TO 125 00238 WRITE (*,*) 'MU < 0 in IGNPOI - ABORT' 00239 WRITE (*,*) 'Value of MU: ',mu 00240 CALL XSTOPX ('MU < 0 in IGNPOI - ABORT') 00241 C JJV added line label here 00242 125 muold = mu 00243 m = max0(1,ifix(mu)) 00244 l = 0 00245 p = exp(-mu) 00246 q = p 00247 p0 = p 00248 C 00249 C STEP U. UNIFORM SAMPLE FOR INVERSION METHOD 00250 C 00251 130 u = ranf() 00252 ignpoi = 0 00253 IF (u.LE.p0) RETURN 00254 C 00255 C STEP T. TABLE COMPARISON UNTIL THE END PP(L) OF THE 00256 C PP-TABLE OF CUMULATIVE POISSON PROBABILITIES 00257 C (0.458=PP(9) FOR MU=10) 00258 C 00259 IF (l.EQ.0) GO TO 150 00260 j = 1 00261 IF (u.GT.0.458) j = min0(l,m) 00262 DO 140 k = j,l 00263 IF (u.LE.pp(k)) GO TO 180 00264 140 CONTINUE 00265 IF (l.EQ.35) GO TO 130 00266 C 00267 C STEP C. CREATION OF NEW POISSON PROBABILITIES P 00268 C AND THEIR CUMULATIVES Q=PP(K) 00269 C 00270 150 l = l + 1 00271 DO 160 k = l,35 00272 p = p*mu/float(k) 00273 q = q + p 00274 pp(k) = q 00275 IF (u.LE.q) GO TO 170 00276 160 CONTINUE 00277 l = 35 00278 GO TO 130 00279 00280 170 l = k 00281 180 ignpoi = k 00282 RETURN 00283 00284 END
1.7.1
