00001 REAL FUNCTION sgamma(a) 00002 C**********************************************************************C 00003 C C 00004 C C 00005 C (STANDARD-) G A M M A DISTRIBUTION C 00006 C C 00007 C C 00008 C**********************************************************************C 00009 C**********************************************************************C 00010 C C 00011 C PARAMETER A >= 1.0 ! C 00012 C C 00013 C**********************************************************************C 00014 C C 00015 C FOR DETAILS SEE: C 00016 C C 00017 C AHRENS, J.H. AND DIETER, U. C 00018 C GENERATING GAMMA VARIATES BY A C 00019 C MODIFIED REJECTION TECHNIQUE. C 00020 C COMM. ACM, 25,1 (JAN. 1982), 47 - 54. C 00021 C C 00022 C STEP NUMBERS CORRESPOND TO ALGORITHM 'GD' IN THE ABOVE PAPER C 00023 C (STRAIGHTFORWARD IMPLEMENTATION) C 00024 C C 00025 C Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of C 00026 C SUNIF. The argument IR thus goes away. C 00027 C C 00028 C**********************************************************************C 00029 C C 00030 C PARAMETER 0.0 < A < 1.0 ! C 00031 C C 00032 C**********************************************************************C 00033 C C 00034 C FOR DETAILS SEE: C 00035 C C 00036 C AHRENS, J.H. AND DIETER, U. C 00037 C COMPUTER METHODS FOR SAMPLING FROM GAMMA, C 00038 C BETA, POISSON AND BINOMIAL DISTRIBUTIONS. C 00039 C COMPUTING, 12 (1974), 223 - 246. C 00040 C C 00041 C (ADAPTED IMPLEMENTATION OF ALGORITHM 'GS' IN THE ABOVE PAPER) C 00042 C C 00043 C**********************************************************************C 00044 C 00045 C 00046 C INPUT: A =PARAMETER (MEAN) OF THE STANDARD GAMMA DISTRIBUTION 00047 C OUTPUT: SGAMMA = SAMPLE FROM THE GAMMA-(A)-DISTRIBUTION 00048 C 00049 C COEFFICIENTS Q(K) - FOR Q0 = SUM(Q(K)*A**(-K)) 00050 C COEFFICIENTS A(K) - FOR Q = Q0+(T*T/2)*SUM(A(K)*V**K) 00051 C COEFFICIENTS E(K) - FOR EXP(Q)-1 = SUM(E(K)*Q**K) 00052 C 00053 C .. Scalar Arguments .. 00054 REAL a 00055 C .. 00056 C .. Local Scalars .. (JJV added B0 to fix rare and subtle bug) 00057 REAL a1,a2,a3,a4,a5,a6,a7,aa,aaa,b,b0,c,d,e,e1,e2,e3,e4,e5,p,q,q0, 00058 + q1,q2,q3,q4,q5,q6,q7,r,s,s2,si,sqrt32,t,u,v,w,x 00059 C .. 00060 C .. External Functions .. 00061 REAL ranf,sexpo,snorm 00062 EXTERNAL ranf,sexpo,snorm 00063 C .. 00064 C .. Intrinsic Functions .. 00065 INTRINSIC abs,alog,exp,sign,sqrt 00066 C .. 00067 C .. Save statement .. 00068 C JJV added Save statement for vars in Data satatements 00069 SAVE aa,aaa,s2,s,d,q0,b,si,c,q1,q2,q3,q4,q5,q6,q7,a1,a2,a3,a4,a5, 00070 + a6,a7,e1,e2,e3,e4,e5,sqrt32 00071 C .. 00072 C .. Data statements .. 00073 C 00074 C PREVIOUS A PRE-SET TO ZERO - AA IS A', AAA IS A" 00075 C SQRT32 IS THE SQUAREROOT OF 32 = 5.656854249492380 00076 C 00077 DATA q1,q2,q3,q4,q5,q6,q7/.04166669,.02083148,.00801191,.00144121, 00078 + -.00007388,.00024511,.00024240/ 00079 DATA a1,a2,a3,a4,a5,a6,a7/.3333333,-.2500030,.2000062,-.1662921, 00080 + .1423657,-.1367177,.1233795/ 00081 DATA e1,e2,e3,e4,e5/1.,.4999897,.1668290,.0407753,.0102930/ 00082 DATA aa/0.0/,aaa/0.0/,sqrt32/5.656854/ 00083 C .. 00084 C .. Executable Statements .. 00085 C 00086 IF (a.EQ.aa) GO TO 10 00087 IF (a.LT.1.0) GO TO 130 00088 C 00089 C STEP 1: RECALCULATIONS OF S2,S,D IF A HAS CHANGED 00090 C 00091 aa = a 00092 s2 = a - 0.5 00093 s = sqrt(s2) 00094 d = sqrt32 - 12.0*s 00095 C 00096 C STEP 2: T=STANDARD NORMAL DEVIATE, 00097 C X=(S,1/2)-NORMAL DEVIATE. 00098 C IMMEDIATE ACCEPTANCE (I) 00099 C 00100 10 t = snorm() 00101 x = s + 0.5*t 00102 sgamma = x*x 00103 IF (t.GE.0.0) RETURN 00104 C 00105 C STEP 3: U= 0,1 -UNIFORM SAMPLE. SQUEEZE ACCEPTANCE (S) 00106 C 00107 u = ranf() 00108 IF (d*u.LE.t*t*t) RETURN 00109 C 00110 C STEP 4: RECALCULATIONS OF Q0,B,SI,C IF NECESSARY 00111 C 00112 IF (a.EQ.aaa) GO TO 40 00113 aaa = a 00114 r = 1.0/a 00115 q0 = ((((((q7*r+q6)*r+q5)*r+q4)*r+q3)*r+q2)*r+q1)*r 00116 C 00117 C APPROXIMATION DEPENDING ON SIZE OF PARAMETER A 00118 C THE CONSTANTS IN THE EXPRESSIONS FOR B, SI AND 00119 C C WERE ESTABLISHED BY NUMERICAL EXPERIMENTS 00120 C 00121 IF (a.LE.3.686) GO TO 30 00122 IF (a.LE.13.022) GO TO 20 00123 C 00124 C CASE 3: A .GT. 13.022 00125 C 00126 b = 1.77 00127 si = .75 00128 c = .1515/s 00129 GO TO 40 00130 C 00131 C CASE 2: 3.686 .LT. A .LE. 13.022 00132 C 00133 20 b = 1.654 + .0076*s2 00134 si = 1.68/s + .275 00135 c = .062/s + .024 00136 GO TO 40 00137 C 00138 C CASE 1: A .LE. 3.686 00139 C 00140 30 b = .463 + s + .178*s2 00141 si = 1.235 00142 c = .195/s - .079 + .16*s 00143 C 00144 C STEP 5: NO QUOTIENT TEST IF X NOT POSITIVE 00145 C 00146 40 IF (x.LE.0.0) GO TO 70 00147 C 00148 C STEP 6: CALCULATION OF V AND QUOTIENT Q 00149 C 00150 v = t/ (s+s) 00151 IF (abs(v).LE.0.25) GO TO 50 00152 q = q0 - s*t + 0.25*t*t + (s2+s2)*alog(1.0+v) 00153 GO TO 60 00154 00155 50 q = q0 + 0.5*t*t* ((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v 00156 C 00157 C STEP 7: QUOTIENT ACCEPTANCE (Q) 00158 C 00159 60 IF (alog(1.0-u).LE.q) RETURN 00160 C 00161 C STEP 8: E=STANDARD EXPONENTIAL DEVIATE 00162 C U= 0,1 -UNIFORM DEVIATE 00163 C T=(B,SI)-DOUBLE EXPONENTIAL (LAPLACE) SAMPLE 00164 C 00165 70 e = sexpo() 00166 u = ranf() 00167 u = u + u - 1.0 00168 t = b + sign(si*e,u) 00169 C 00170 C STEP 9: REJECTION IF T .LT. TAU(1) = -.71874483771719 00171 C 00172 80 IF (t.LT. (-.7187449)) GO TO 70 00173 C 00174 C STEP 10: CALCULATION OF V AND QUOTIENT Q 00175 C 00176 v = t/ (s+s) 00177 IF (abs(v).LE.0.25) GO TO 90 00178 q = q0 - s*t + 0.25*t*t + (s2+s2)*alog(1.0+v) 00179 GO TO 100 00180 00181 90 q = q0 + 0.5*t*t* ((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v 00182 C 00183 C STEP 11: HAT ACCEPTANCE (H) (IF Q NOT POSITIVE GO TO STEP 8) 00184 C 00185 100 IF (q.LE.0.0) GO TO 70 00186 IF (q.LE.0.5) GO TO 110 00187 C 00188 C JJV modified the code through line 125 to handle large Q case 00189 C 00190 IF (q.LT.15.0) GO TO 105 00191 C 00192 C JJV Here Q is large enough that Q = log(exp(Q) - 1.0) (for real Q) 00193 C JJV so reformulate test at 120 in terms of one EXP, if not too big 00194 C JJV 87.49823 is close to the largest real which can be 00195 C JJV exponentiated (87.49823 = log(1.0E38)) 00196 C 00197 IF ((q+e-0.5*t*t).GT.87.49823) GO TO 125 00198 IF (c*abs(u).GT.exp(q+e-0.5*t*t)) GO TO 70 00199 GO TO 125 00200 00201 105 w = exp(q) - 1.0 00202 GO TO 120 00203 00204 110 w = ((((e5*q+e4)*q+e3)*q+e2)*q+e1)*q 00205 C 00206 C IF T IS REJECTED, SAMPLE AGAIN AT STEP 8 00207 C 00208 120 IF (c*abs(u).GT.w*exp(e-0.5*t*t)) GO TO 70 00209 125 x = s + 0.5*t 00210 sgamma = x*x 00211 RETURN 00212 C 00213 C ALTERNATE METHOD FOR PARAMETERS A BELOW 1 (.3678794=EXP(-1.)) 00214 C 00215 C JJV changed B to B0 (which was added to declarations for this) 00216 C JJV in 130 to END to fix rare and subtle bug. 00217 C JJV Line: '130 aa = 0.0' was removed (unnecessary, wasteful). 00218 C JJV Reasons: the state of AA only serves to tell the A .GE. 1.0 00219 C JJV case if certain A-dependant constants need to be recalculated. 00220 C JJV The A .LT. 1.0 case (here) no longer changes any of these, and 00221 C JJV the recalculation of B (which used to change with an 00222 C JJV A .LT. 1.0 call) is governed by the state of AAA anyway. 00223 C 00224 130 b0 = 1.0 + .3678794*a 00225 140 p = b0*ranf() 00226 IF (p.GE.1.0) GO TO 150 00227 sgamma = exp(alog(p)/a) 00228 IF (sexpo().LT.sgamma) GO TO 140 00229 RETURN 00230 00231 150 sgamma = -alog((b0-p)/a) 00232 IF (sexpo().LT. (1.0-a)*alog(sgamma)) GO TO 140 00233 RETURN 00234 00235 END
1.7.1
