GNU Octave  4.0.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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
base_det::exp
int exp(void) const
Definition: DET.h:64
abs
T abs(T x)
Definition: pr-output.cc:3062
octave_value::sqrt
octave_value sqrt(void) const
Definition: ov.h:1200
Generated on Thu Jun 4 2015 23:30:24 for GNU Octave by   doxygen 1.8.8