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