1 DOUBLE PRECISION FUNCTION sgamma(a) 2C**********************************************************************C 3C C 4C C 5C (STANDARD-) G A M M A DISTRIBUTION C 6C C 7C C 8C**********************************************************************C 9C**********************************************************************C 10C C 11C PARAMETER A >= 1.0 ! C 12C C 13C**********************************************************************C 14C C 15C FOR DETAILS SEE: C 16C C 17C AHRENS, J.H. AND DIETER, U. C 18C GENERATING GAMMA VARIATES BY A C 19C MODIFIED REJECTION TECHNIQUE. C 20C COMM. ACM, 25,1 (JAN. 1982), 47 - 54. C 21C C 22C STEP NUMBERS CORRESPOND TO ALGORITHM 'GD' IN THE ABOVE PAPER C 23C (STRAIGHTFORWARD IMPLEMENTATION) C 24C C 25C Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of C 26C SUNIF. The argument IR thus goes away. C 27C C 28C**********************************************************************C 29C C 30C PARAMETER 0.0 < A < 1.0 ! C 31C C 32C**********************************************************************C 33C C 34C FOR DETAILS SEE: C 35C C 36C AHRENS, J.H. AND DIETER, U. C 37C COMPUTER METHODS FOR SAMPLING FROM GAMMA, C 38C BETA, POISSON AND BINOMIAL DISTRIBUTIONS. C 39C COMPUTING, 12 (1974), 223 - 246. C 40C C 41C (ADAPTED IMPLEMENTATION OF ALGORITHM 'GS' IN THE ABOVE PAPER) C 42C C 43C**********************************************************************C 44C 45C 46C INPUT: A =PARAMETER (MEAN) OF THE STANDARD GAMMA DISTRIBUTION 47C OUTPUT: SGAMMA = SAMPLE FROM THE GAMMA-(A)-DISTRIBUTION 48C 49C COEFFICIENTS Q(K) - FOR Q0 = SUM(Q(K)*A**(-K)) 50C COEFFICIENTS A(K) - FOR Q = Q0+(T*T/2)*SUM(A(K)*V**K) 51C COEFFICIENTS E(K) - FOR EXP(Q)-1 = SUM(E(K)*Q**K) 52C 53C .. Scalar Arguments .. 54 DOUBLE PRECISION a 55C .. 56C .. Local Scalars .. (JJV added B0 to fix rare and subtle bug) 57 DOUBLE PRECISION a1,a2,a3,a4,a5,a6,a7,aa,aaa,b,b0,c,d, 58 + e,e1,e2,e3,e4,e5,p,q,q0,q1,q2,q3,q4,q5,q6,q7,r,s,s2,si,sqrt32, 59 + t,u,v,w,x 60C .. 61C .. External Functions .. 62 DOUBLE PRECISION ranf,sexpo,snorm 63 EXTERNAL ranf,sexpo,snorm 64C .. 65C .. Intrinsic Functions .. 66C INTRINSIC abs,alog,exp,sign,sqrt 67C .. 68C .. Save statement .. 69C JJV added Save statement for vars in Data satatements 70 SAVE aa,aaa,s2,s,d,q0,b,si,c,q1,q2,q3,q4,q5,q6,q7,a1,a2,a3,a4,a5, 71 + a6,a7,e1,e2,e3,e4,e5,sqrt32 72C .. 73C .. Data statements .. 74C 75C PREVIOUS A PRE-SET TO ZERO - AA IS A', AAA IS A" 76C SQRT32 IS THE SQUAREROOT OF 32 = 5.656854249492380 77C 78 DATA q1,q2,q3,q4,q5,q6,q7/.04166669,.02083148,.00801191,.00144121, 79 + -.00007388,.00024511,.00024240/ 80 DATA a1,a2,a3,a4,a5,a6,a7/.3333333,-.2500030,.2000062,-.1662921, 81 + .1423657,-.1367177,.1233795/ 82 DATA e1,e2,e3,e4,e5/1.,.4999897,.1668290,.0407753,.0102930/ 83 DATA aa/0.0/,aaa/0.0/,sqrt32/5.656854/ 84C .. 85C .. Executable Statements .. 86C 87 IF (a.EQ.aa) GO TO 10 88 IF (a.LT.1.0) GO TO 130 89C 90C STEP 1: RECALCULATIONS OF S2,S,D IF A HAS CHANGED 91C 92 aa = a 93 s2 = a - 0.5 94 s = sqrt(s2) 95 d = sqrt32 - 12.0*s 96C 97C STEP 2: T=STANDARD NORMAL DEVIATE, 98C X=(S,1/2)-NORMAL DEVIATE. 99C IMMEDIATE ACCEPTANCE (I) 100C 101 10 t = snorm() 102 x = s + 0.5*t 103 sgamma = x*x 104 IF (t.GE.0.0) RETURN 105C 106C STEP 3: U= 0,1 -UNIFORM SAMPLE. SQUEEZE ACCEPTANCE (S) 107C 108 u = ranf() 109 IF (d*u.LE.t*t*t) RETURN 110C 111C STEP 4: RECALCULATIONS OF Q0,B,SI,C IF NECESSARY 112C 113 IF (a.EQ.aaa) GO TO 40 114 aaa = a 115 r = 1.0/a 116 q0 = ((((((q7*r+q6)*r+q5)*r+q4)*r+q3)*r+q2)*r+q1)*r 117C 118C APPROXIMATION DEPENDING ON SIZE OF PARAMETER A 119C THE CONSTANTS IN THE EXPRESSIONS FOR B, SI AND 120C C WERE ESTABLISHED BY NUMERICAL EXPERIMENTS 121C 122 IF (a.LE.3.686) GO TO 30 123 IF (a.LE.13.022) GO TO 20 124C 125C CASE 3: A .GT. 13.022 126C 127 b = 1.77 128 si = .75 129 c = .1515/s 130 GO TO 40 131C 132C CASE 2: 3.686 .LT. A .LE. 13.022 133C 134 20 b = 1.654 + .0076*s2 135 si = 1.68/s + .275 136 c = .062/s + .024 137 GO TO 40 138C 139C CASE 1: A .LE. 3.686 140C 141 30 b = .463 + s + .178*s2 142 si = 1.235 143 c = .195/s - .079 + .16*s 144C 145C STEP 5: NO QUOTIENT TEST IF X NOT POSITIVE 146C 147 40 IF (x.LE.0.0) GO TO 70 148C 149C STEP 6: CALCULATION OF V AND QUOTIENT Q 150C 151 v = t/(s+s) 152 IF (abs(v).LE.0.25) GO TO 50 153 q = q0 - s*t + 0.25*t*t + (s2+s2)*dlog(1.0+v) 154 GO TO 60 155 156 50 q = q0 + 0.5*t*t* ((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v 157C 158C STEP 7: QUOTIENT ACCEPTANCE (Q) 159C 160 60 IF (dlog(1.0-u).LE.q) RETURN 161C 162C STEP 8: E=STANDARD EXPONENTIAL DEVIATE 163C U= 0,1 -UNIFORM DEVIATE 164C T=(B,SI)-DOUBLE EXPONENTIAL (LAPLACE) SAMPLE 165C 166 70 e = sexpo() 167 u = ranf() 168 u = u + u - 1.0 169 t = b + sign(si*e,u) 170C 171C STEP 9: REJECTION IF T .LT. TAU(1) = -.71874483771719 172C 173 80 IF (t.LT. (-.7187449)) GO TO 70 174C 175C STEP 10: CALCULATION OF V AND QUOTIENT Q 176C 177 v = t/ (s+s) 178 IF (abs(v).LE.0.25) GO TO 90 179 q = q0 - s*t + 0.25*t*t + (s2+s2)*dlog(1.0+v) 180 GO TO 100 181 182 90 q = q0 + 0.5*t*t* ((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v 183C 184C STEP 11: HAT ACCEPTANCE (H) (IF Q NOT POSITIVE GO TO STEP 8) 185C 186 100 IF (q.LE.0.0) GO TO 70 187 IF (q.LE.0.5) GO TO 110 188C 189C JJV modified the code through line 125 to handle large Q case 190C 191 IF (q.LT.15.0) GO TO 105 192C 193C JJV Here Q is large enough that Q = log(exp(Q) - 1.0) (for real Q) 194C JJV so reformulate test at 120 in terms of one EXP, if not too big 195C JJV 87.49823 is close to the largest real which can be 196C JJV exponentiated (87.49823 = log(1.0E38)) 197C 198 IF ((q+e-0.5*t*t).GT.87.49823) GO TO 125 199 IF (c*abs(u).GT.exp(q+e-0.5*t*t)) GO TO 70 200 GO TO 125 201 202 105 w = exp(q) - 1.0 203 GO TO 120 204 205 110 w = ((((e5*q+e4)*q+e3)*q+e2)*q+e1)*q 206C 207C IF T IS REJECTED, SAMPLE AGAIN AT STEP 8 208C 209 120 IF (c*abs(u).GT.w*exp(e-0.5*t*t)) GO TO 70 210 125 x = s + 0.5*t 211 sgamma = x*x 212 RETURN 213C 214C ALTERNATE METHOD FOR PARAMETERS A BELOW 1 (.3678794=EXP(-1.)) 215C 216C JJV changed B to B0 (which was added to declarations for this) 217C JJV in 130 to END to fix rare and subtle bug. 218C JJV Line: '130 aa = 0.0' was removed (unnecessary, wasteful). 219C JJV Reasons: the state of AA only serves to tell the A .GE. 1.0 220C JJV case if certain A-dependant constants need to be recalculated. 221C JJV The A .LT. 1.0 case (here) no longer changes any of these, and 222C JJV the recalculation of B (which used to change with an 223C JJV A .LT. 1.0 call) is governed by the state of AAA anyway. 224C 225 130 b0 = 1.0 + .3678794*a 226 140 p = b0*ranf() 227 IF (p.GE.1.0) GO TO 150 228 sgamma = exp(dlog(p)/a) 229 IF (sexpo().LT.sgamma) GO TO 140 230 RETURN 231 232 150 sgamma = -dlog((b0-p)/a) 233 IF (sexpo().LT. (1.0-a)*dlog(sgamma)) GO TO 140 234 RETURN 235 236 END 237