GNU Octave  4.4.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
sgamma.f
Go to the documentation of this file.
1  REAL FUNCTION sgamma(a)
2 C**********************************************************************C
3 C C
4 C C
5 C (STANDARD-) G A M M A DISTRIBUTION C
6 C C
7 C C
8 C**********************************************************************C
9 C**********************************************************************C
10 C C
11 C PARAMETER A >= 1.0 ! C
12 C C
13 C**********************************************************************C
14 C C
15 C FOR DETAILS SEE: C
16 C C
17 C AHRENS, J.H. AND DIETER, U. C
18 C GENERATING GAMMA VARIATES BY A C
19 C MODIFIED REJECTION TECHNIQUE. C
20 C COMM. ACM, 25,1 (JAN. 1982), 47 - 54. C
21 C C
22 C STEP NUMBERS CORRESPOND TO ALGORITHM 'GD' IN THE ABOVE PAPER C
23 C (STRAIGHTFORWARD IMPLEMENTATION) C
24 C C
25 C Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of C
26 C SUNIF. The argument IR thus goes away. C
27 C C
28 C**********************************************************************C
29 C C
30 C PARAMETER 0.0 < A < 1.0 ! C
31 C C
32 C**********************************************************************C
33 C C
34 C FOR DETAILS SEE: C
35 C C
36 C AHRENS, J.H. AND DIETER, U. C
37 C COMPUTER METHODS FOR SAMPLING FROM GAMMA, C
38 C BETA, POISSON AND BINOMIAL DISTRIBUTIONS. C
39 C COMPUTING, 12 (1974), 223 - 246. C
40 C C
41 C (ADAPTED IMPLEMENTATION OF ALGORITHM 'GS' IN THE ABOVE PAPER) C
42 C C
43 C**********************************************************************C
44 C
45 C
46 C INPUT: A =PARAMETER (MEAN) OF THE STANDARD GAMMA DISTRIBUTION
47 C OUTPUT: SGAMMA = SAMPLE FROM THE GAMMA-(A)-DISTRIBUTION
48 C
49 C COEFFICIENTS Q(K) - FOR Q0 = SUM(Q(K)*A**(-K))
50 C COEFFICIENTS A(K) - FOR Q = Q0+(T*T/2)*SUM(A(K)*V**K)
51 C COEFFICIENTS E(K) - FOR EXP(Q)-1 = SUM(E(K)*Q**K)
52 C
53 C .. Scalar Arguments ..
54  REAL a
55 C ..
56 C .. 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
59 C ..
60 C .. External Functions ..
61  REAL ranf,sexpo,snorm
62  EXTERNAL ranf,sexpo,snorm
63 C ..
64 C .. Intrinsic Functions ..
65  INTRINSIC abs,alog,exp,sign,sqrt
66 C ..
67 C .. Save statement ..
68 C 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
71 C ..
72 C .. Data statements ..
73 C
74 C PREVIOUS A PRE-SET TO ZERO - AA IS A', AAA IS A"
75 C SQRT32 IS THE SQUAREROOT OF 32 = 5.656854249492380
76 C
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/
83 C ..
84 C .. Executable Statements ..
85 C
86  IF (a.EQ.aa) GO TO 10
87  IF (a.LT.1.0) GO TO 130
88 C
89 C STEP 1: RECALCULATIONS OF S2,S,D IF A HAS CHANGED
90 C
91  aa = a
92  s2 = a - 0.5
93  s = sqrt(s2)
94  d = sqrt32 - 12.0*s
95 C
96 C STEP 2: T=STANDARD NORMAL DEVIATE,
97 C X=(S,1/2)-NORMAL DEVIATE.
98 C IMMEDIATE ACCEPTANCE (I)
99 C
100  10 t = snorm()
101  x = s + 0.5*t
102  sgamma = x*x
103  IF (t.GE.0.0) RETURN
104 C
105 C STEP 3: U= 0,1 -UNIFORM SAMPLE. SQUEEZE ACCEPTANCE (S)
106 C
107  u = ranf()
108  IF (d*u.LE.t*t*t) RETURN
109 C
110 C STEP 4: RECALCULATIONS OF Q0,B,SI,C IF NECESSARY
111 C
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
116 C
117 C APPROXIMATION DEPENDING ON SIZE OF PARAMETER A
118 C THE CONSTANTS IN THE EXPRESSIONS FOR B, SI AND
119 C C WERE ESTABLISHED BY NUMERICAL EXPERIMENTS
120 C
121  IF (a.LE.3.686) GO TO 30
122  IF (a.LE.13.022) GO TO 20
123 C
124 C CASE 3: A .GT. 13.022
125 C
126  b = 1.77
127  si = .75
128  c = .1515/s
129  GO TO 40
130 C
131 C CASE 2: 3.686 .LT. A .LE. 13.022
132 C
133  20 b = 1.654 + .0076*s2
134  si = 1.68/s + .275
135  c = .062/s + .024
136  GO TO 40
137 C
138 C CASE 1: A .LE. 3.686
139 C
140  30 b = .463 + s + .178*s2
141  si = 1.235
142  c = .195/s - .079 + .16*s
143 C
144 C STEP 5: NO QUOTIENT TEST IF X NOT POSITIVE
145 C
146  40 IF (x.LE.0.0) GO TO 70
147 C
148 C STEP 6: CALCULATION OF V AND QUOTIENT Q
149 C
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
156 C
157 C STEP 7: QUOTIENT ACCEPTANCE (Q)
158 C
159  60 IF (alog(1.0-u).LE.q) RETURN
160 C
161 C STEP 8: E=STANDARD EXPONENTIAL DEVIATE
162 C U= 0,1 -UNIFORM DEVIATE
163 C T=(B,SI)-DOUBLE EXPONENTIAL (LAPLACE) SAMPLE
164 C
165  70 e = sexpo()
166  u = ranf()
167  u = u + u - 1.0
168  t = b + sign(si*e,u)
169 C
170 C STEP 9: REJECTION IF T .LT. TAU(1) = -.71874483771719
171 C
172  80 IF (t.LT. (-.7187449)) GO TO 70
173 C
174 C STEP 10: CALCULATION OF V AND QUOTIENT Q
175 C
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
182 C
183 C STEP 11: HAT ACCEPTANCE (H) (IF Q NOT POSITIVE GO TO STEP 8)
184 C
185  100 IF (q.LE.0.0) GO TO 70
186  IF (q.LE.0.5) GO TO 110
187 C
188 C JJV modified the code through line 125 to handle large Q case
189 C
190  IF (q.LT.15.0) GO TO 105
191 C
192 C JJV Here Q is large enough that Q = log(exp(Q) - 1.0) (for real Q)
193 C JJV so reformulate test at 120 in terms of one EXP, if not too big
194 C JJV 87.49823 is close to the largest real which can be
195 C JJV exponentiated (87.49823 = log(1.0E38))
196 C
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
205 C
206 C IF T IS REJECTED, SAMPLE AGAIN AT STEP 8
207 C
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
212 C
213 C ALTERNATE METHOD FOR PARAMETERS A BELOW 1 (.3678794=EXP(-1.))
214 C
215 C JJV changed B to B0 (which was added to declarations for this)
216 C JJV in 130 to END to fix rare and subtle bug.
217 C JJV Line: '130 aa = 0.0' was removed (unnecessary, wasteful).
218 C JJV Reasons: the state of AA only serves to tell the A .GE. 1.0
219 C JJV case if certain A-dependant constants need to be recalculated.
220 C JJV The A .LT. 1.0 case (here) no longer changes any of these, and
221 C JJV the recalculation of B (which used to change with an
222 C JJV A .LT. 1.0 call) is governed by the state of AAA anyway.
223 C
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
sexpo
real function sexpo()
Definition: sexpo.f:2
abs
static T abs(T x)
Definition: pr-output.cc:1696
ranf
real function ranf()
Definition: ranf.f:2
snorm
real function snorm()
Definition: snorm.f:2
sgamma
real function sgamma(a)
Definition: sgamma.f:2