blitz
Version 1.0.2
Loading...
Searching...
No Matches
beta.h
Go to the documentation of this file.
1
// -*- C++ -*-
2
// $Id$
3
4
/*
5
* Generate Beta random deviate
6
*
7
* Returns a single random deviate from the beta distribution with
8
* parameters A and B. The density of the beta is
9
* x^(a-1) * (1-x)^(b-1) / B(a,b) for 0 < x < 1
10
*
11
* The mean is a/(a+b).
12
* The variance is ab/((a+b)^2(a+b+1))
13
* The rth moment is (a+r-1)^(r)/(a+b+r-1)^(r)
14
* where a^(r) means a * (a-1) * (a-2) * ... * (a-r+1)
15
*
16
* Method
17
* R. C. H. Cheng
18
* Generating Beta Variates with Nonintegral Shape Parameters
19
* Communications of the ACM, 21:317-322 (1978)
20
* (Algorithms BB and BC)
21
* http://www.acm.org/pubs/citations/journals/toms/1991-17-1/p98-l_ecuyer/
22
*
23
* This class has been adapted from RANDLIB.C 1.3, by
24
* Barry W. Brown, James Lovato, Kathy Russell, and John Venier.
25
*
26
* Adapter's note (TV): This code has gone from Pascal to Fortran to C.
27
* As a result it is a bit of a mess. Note also that the algorithms were
28
* originally designed for 32-bit float, and so some of the constants
29
* below have inadequate precision. This will not cause problems for
30
* casual use, but if you are generating millions of beta variates and
31
* rely on some convergence property, you may have want to worry
32
* about this.
33
*
34
* NEEDS_WORK: dig out the original paper and determine these constants
35
* to precision adequate for 128-bit float.
36
* NEEDS_WORK: turn this into structured code.
37
*/
38
39
#ifndef BZ_RANDOM_BETA
40
#define BZ_RANDOM_BETA
41
42
#ifndef BZ_RANDOM_UNIFORM
43
#include <
random/uniform.h
>
44
#endif
45
46
#ifndef BZ_NUMINQUIRE_H
47
#include <blitz/numinquire.h>
48
#endif
49
50
namespace
ranlib
{
51
52
template
<
typename
T = double,
typename
IRNG =
defaultIRNG
,
53
typename
stateTag =
defaultState
>
54
class
Beta
:
public
UniformOpen
<T,IRNG,stateTag>
55
{
56
public
:
57
typedef
T
T_numtype
;
58
59
Beta
(T a, T b)
60
{
61
setParameters
(a, b);
62
}
63
64
Beta
(T a, T b,
unsigned
int
i) :
UniformOpen
<T, IRNG, stateTag>(i)
65
{
66
setParameters
(a, b);
67
}
68
69
T
random
();
70
71
void
setParameters
(T a, T b)
72
{
73
aa
= a;
74
bb
= b;
75
infnty
= 0.3 * blitz::huge(T());
76
minlog
= 0.085 * blitz::tiny(T());
77
expmax
= log(
infnty
);
78
}
79
80
protected
:
81
T
ranf
()
82
{
83
return
UniformOpen<T,IRNG,stateTag>::random
();
84
}
85
86
T
aa
,
bb
;
87
T
infnty
,
minlog
,
expmax
;
88
};
89
90
template
<
typename
T,
typename
IRNG,
typename
stateTag>
91
T
Beta<T,IRNG,stateTag>::random
()
92
{
93
/* JJV changed expmax (log(1.0E38)==87.49823), and added minlog */
94
95
// TV: The original code had infnty = 1.0E38, minlog = 1.0E-37.
96
// I have changed these to 0.3 * huge and 8.5 * tiny. For IEEE
97
// float this works out to 1.020847E38 and 0.9991702E-37.
98
// I changed expmax from 87.49823 to log(infnty), which works out
99
// to 87.518866 for float.
100
101
static
T olda =
minlog
;
102
static
T oldb =
minlog
;
103
static
T genbet,a,alpha,b,beta,delta,gamma,k1,k2,r,s,t,u1,u2,v,w,y,z;
104
static
long
qsame;
105
106
// This code determines if the aa and bb parameters are unchanged.
107
// If so, some code can be skipped.
108
109
qsame = (olda ==
aa
) && (oldb ==
bb
);
110
111
if
(!qsame)
112
{
113
BZPRECHECK((
aa
>
minlog
) && (
bb
>
minlog
),
114
"Beta RNG: parameters must be >= "
<<
minlog
<< endl
115
<<
"Parameters are aa="
<<
aa
<<
" and bb="
<<
bb
<< endl);
116
olda =
aa
;
117
oldb =
bb
;
118
}
119
120
if
(!(min(
aa
,
bb
) > 1.0))
121
goto
S100;
122
/*
123
Alborithm BB
124
Initialize
125
*/
126
if
(qsame)
goto
S30;
127
a = min(
aa
,
bb
);
128
b = max(
aa
,
bb
);
129
alpha = a+b;
130
beta = sqrt((alpha-2.0)/(2.0*a*b-alpha));
131
gamma = a+1.0/beta;
132
S30:
133
S40:
134
u1 =
ranf
();
135
/*
136
Step 1
137
*/
138
u2 =
ranf
();
139
v = beta*log(u1/(1.0-u1));
140
/* JJV altered this */
141
if
(v >
expmax
)
goto
S55;
142
/*
143
* JJV added checker to see if a*exp(v) will overflow
144
* JJV S50 _was_ w = a*exp(v); also note here a > 1.0
145
*/
146
w = exp(v);
147
if
(w >
infnty
/a)
goto
S55;
148
w *= a;
149
goto
S60;
150
S55:
151
w =
infnty
;
152
S60:
153
z = pow(u1,2.0)*u2;
154
r = gamma*v-1.3862944;
155
s = a+r-w;
156
/*
157
Step 2
158
*/
159
if
(s+2.609438 >= 5.0*z)
goto
S70;
160
/*
161
Step 3
162
*/
163
t = log(z);
164
if
(s > t)
goto
S70;
165
/*
166
* Step 4
167
*
168
* JJV added checker to see if log(alpha/(b+w)) will
169
* JJV overflow. If so, we count the log as -INF, and
170
* JJV consequently evaluate conditional as true, i.e.
171
* JJV the algorithm rejects the trial and starts over
172
* JJV May not need this here since alpha > 2.0
173
*/
174
if
(alpha/(b+w) <
minlog
)
goto
S40;
175
if
(r+alpha*log(alpha/(b+w)) < t)
goto
S40;
176
S70:
177
/*
178
Step 5
179
*/
180
if
(!(
aa
== a))
goto
S80;
181
genbet = w/(b+w);
182
goto
S90;
183
S80:
184
genbet = b/(b+w);
185
S90:
186
goto
S230;
187
S100:
188
/*
189
Algorithm BC
190
Initialize
191
*/
192
if
(qsame)
goto
S110;
193
a = max(
aa
,
bb
);
194
b = min(
aa
,
bb
);
195
alpha = a+b;
196
beta = 1.0/b;
197
delta = 1.0+a-b;
198
k1 = delta*(1.38889E-2+4.16667E-2*b)/(a*beta-0.777778);
199
k2 = 0.25+(0.5+0.25/delta)*b;
200
S110:
201
S120:
202
u1 =
ranf
();
203
/*
204
Step 1
205
*/
206
u2 =
ranf
();
207
if
(u1 >= 0.5)
goto
S130;
208
/*
209
Step 2
210
*/
211
y = u1*u2;
212
z = u1*y;
213
if
(0.25*u2+z-y >= k1)
goto
S120;
214
goto
S170;
215
S130:
216
/*
217
Step 3
218
*/
219
z = pow(u1,2.0)*u2;
220
if
(!(z <= 0.25))
goto
S160;
221
v = beta*log(u1/(1.0-u1));
222
/*
223
* JJV instead of checking v > expmax at top, I will check
224
* JJV if a < 1, then check the appropriate values
225
*/
226
if
(a > 1.0)
goto
S135;
227
/* JJV a < 1 so it can help out if exp(v) would overflow */
228
if
(v >
expmax
)
goto
S132;
229
w = a*exp(v);
230
goto
S200;
231
S132:
232
w = v + log(a);
233
if
(w >
expmax
)
goto
S140;
234
w = exp(w);
235
goto
S200;
236
S135:
237
/* JJV in this case a > 1 */
238
if
(v >
expmax
)
goto
S140;
239
w = exp(v);
240
if
(w >
infnty
/a)
goto
S140;
241
w *= a;
242
goto
S200;
243
S140:
244
w =
infnty
;
245
goto
S200;
246
/*
247
* JJV old code
248
* if(!(v > expmax)) goto S140;
249
* w = infnty;
250
* goto S150;
251
*S140:
252
* w = a*exp(v);
253
*S150:
254
* goto S200;
255
*/
256
S160:
257
if
(z >= k2)
goto
S120;
258
S170:
259
/*
260
Step 4
261
Step 5
262
*/
263
v = beta*log(u1/(1.0-u1));
264
/* JJV same kind of checking as above */
265
if
(a > 1.0)
goto
S175;
266
/* JJV a < 1 so it can help out if exp(v) would overflow */
267
if
(v >
expmax
)
goto
S172;
268
w = a*exp(v);
269
goto
S190;
270
S172:
271
w = v + log(a);
272
if
(w >
expmax
)
goto
S180;
273
w = exp(w);
274
goto
S190;
275
S175:
276
/* JJV in this case a > 1.0 */
277
if
(v >
expmax
)
goto
S180;
278
w = exp(v);
279
if
(w >
infnty
/a)
goto
S180;
280
w *= a;
281
goto
S190;
282
S180:
283
w =
infnty
;
284
/*
285
* JJV old code
286
* if(!(v > expmax)) goto S180;
287
* w = infnty;
288
* goto S190;
289
*S180:
290
* w = a*exp(v);
291
*/
292
S190:
293
/*
294
* JJV here we also check to see if log overlows; if so, we treat it
295
* JJV as -INF, which means condition is true, i.e. restart
296
*/
297
if
(alpha/(b+w) <
minlog
)
goto
S120;
298
if
(alpha*(log(alpha/(b+w))+v)-1.3862944 < log(z))
goto
S120;
299
S200:
300
/*
301
Step 6
302
*/
303
if
(!(a ==
aa
))
goto
S210;
304
genbet = w/(b+w);
305
goto
S220;
306
S210:
307
genbet = b/(b+w);
308
S230:
309
S220:
310
return
genbet;
311
}
312
313
}
314
315
#endif
// BZ_RANDOM_BETA
ranlib::Beta::random
T random()
Definition
beta.h:91
ranlib::Beta::minlog
T minlog
Definition
beta.h:87
ranlib::Beta::setParameters
void setParameters(T a, T b)
Definition
beta.h:71
ranlib::Beta::Beta
Beta(T a, T b)
Definition
beta.h:59
ranlib::Beta::aa
T aa
Definition
beta.h:86
ranlib::Beta::Beta
Beta(T a, T b, unsigned int i)
Definition
beta.h:64
ranlib::Beta::T_numtype
T T_numtype
Definition
beta.h:57
ranlib::Beta::infnty
T infnty
Definition
beta.h:87
ranlib::Beta::ranf
T ranf()
Definition
beta.h:81
ranlib::Beta::expmax
T expmax
Definition
beta.h:87
ranlib::Beta::bb
T bb
Definition
beta.h:86
ranlib::UniformOpen< double, defaultIRNG, defaultState >::UniformOpen
UniformOpen()
Definition
uniform.h:381
ranlib::UniformOpen::random
T random()
Definition
uniform.h:385
ranlib
Definition
beta.h:50
ranlib::defaultState
sharedState defaultState
Definition
default.h:55
ranlib::defaultIRNG
MersenneTwister defaultIRNG
Definition
default.h:120
uniform.h
random
beta.h
Generated by
1.13.1