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