00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036 #ifndef BZ_RANDOM_BETA
00037 #define BZ_RANDOM_BETA
00038
00039 #ifndef BZ_RANDOM_UNIFORM
00040 #include <random/uniform.h>
00041 #endif
00042
00043 #ifndef BZ_NUMINQUIRE_H
00044 #include <blitz/numinquire.h>
00045 #endif
00046
00047 BZ_NAMESPACE(ranlib)
00048
00049 template<typename T = double, typename IRNG = defaultIRNG,
00050 typename stateTag = defaultState>
00051 class Beta : public UniformOpen<T,IRNG,stateTag>
00052 {
00053 public:
00054 typedef T T_numtype;
00055
00056 Beta(T a, T b)
00057 {
00058 setParameters(a, b);
00059 }
00060
00061 Beta(T a, T b, unsigned int i) : UniformOpen<T, IRNG, stateTag>(i)
00062 {
00063 setParameters(a, b);
00064 }
00065
00066 T random();
00067
00068 void setParameters(T a, T b)
00069 {
00070 aa = a;
00071 bb = b;
00072 infnty = 0.3 * blitz::huge(T());
00073 minlog = 0.085 * blitz::tiny(T());
00074 expmax = log(infnty);
00075 }
00076
00077 protected:
00078 T ranf()
00079 {
00080 return UniformOpen<T,IRNG,stateTag>::random();
00081 }
00082
00083 T aa, bb;
00084 T infnty, minlog, expmax;
00085 };
00086
00087 template<typename T, typename IRNG, typename stateTag>
00088 T Beta<T,IRNG,stateTag>::random()
00089 {
00090
00091
00092
00093
00094
00095
00096
00097
00098 static T olda = minlog;
00099 static T oldb = minlog;
00100 static T genbet,a,alpha,b,beta,delta,gamma,k1,k2,r,s,t,u1,u2,v,w,y,z;
00101 static long qsame;
00102
00103
00104
00105
00106 qsame = (olda == aa) && (oldb == bb);
00107
00108 if (!qsame)
00109 {
00110 BZPRECHECK((aa > minlog) && (bb > minlog),
00111 "Beta RNG: parameters must be >= " << minlog << endl
00112 << "Parameters are aa=" << aa << " and bb=" << bb << endl);
00113 olda = aa;
00114 oldb = bb;
00115 }
00116
00117 if (!(min(aa,bb) > 1.0))
00118 goto S100;
00119
00120
00121
00122
00123 if(qsame) goto S30;
00124 a = min(aa,bb);
00125 b = max(aa,bb);
00126 alpha = a+b;
00127 beta = sqrt((alpha-2.0)/(2.0*a*b-alpha));
00128 gamma = a+1.0/beta;
00129 S30:
00130 S40:
00131 u1 = ranf();
00132
00133
00134
00135 u2 = ranf();
00136 v = beta*log(u1/(1.0-u1));
00137
00138 if(v > expmax) goto S55;
00139
00140
00141
00142
00143 w = exp(v);
00144 if(w > infnty/a) goto S55;
00145 w *= a;
00146 goto S60;
00147 S55:
00148 w = infnty;
00149 S60:
00150 z = pow(u1,2.0)*u2;
00151 r = gamma*v-1.3862944;
00152 s = a+r-w;
00153
00154
00155
00156 if(s+2.609438 >= 5.0*z) goto S70;
00157
00158
00159
00160 t = log(z);
00161 if(s > t) goto S70;
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171 if(alpha/(b+w) < minlog) goto S40;
00172 if(r+alpha*log(alpha/(b+w)) < t) goto S40;
00173 S70:
00174
00175
00176
00177 if(!(aa == a)) goto S80;
00178 genbet = w/(b+w);
00179 goto S90;
00180 S80:
00181 genbet = b/(b+w);
00182 S90:
00183 goto S230;
00184 S100:
00185
00186
00187
00188
00189 if(qsame) goto S110;
00190 a = max(aa,bb);
00191 b = min(aa,bb);
00192 alpha = a+b;
00193 beta = 1.0/b;
00194 delta = 1.0+a-b;
00195 k1 = delta*(1.38889E-2+4.16667E-2*b)/(a*beta-0.777778);
00196 k2 = 0.25+(0.5+0.25/delta)*b;
00197 S110:
00198 S120:
00199 u1 = ranf();
00200
00201
00202
00203 u2 = ranf();
00204 if(u1 >= 0.5) goto S130;
00205
00206
00207
00208 y = u1*u2;
00209 z = u1*y;
00210 if(0.25*u2+z-y >= k1) goto S120;
00211 goto S170;
00212 S130:
00213
00214
00215
00216 z = pow(u1,2.0)*u2;
00217 if(!(z <= 0.25)) goto S160;
00218 v = beta*log(u1/(1.0-u1));
00219
00220
00221
00222
00223 if(a > 1.0) goto S135;
00224
00225 if(v > expmax) goto S132;
00226 w = a*exp(v);
00227 goto S200;
00228 S132:
00229 w = v + log(a);
00230 if(w > expmax) goto S140;
00231 w = exp(w);
00232 goto S200;
00233 S135:
00234
00235 if(v > expmax) goto S140;
00236 w = exp(v);
00237 if(w > infnty/a) goto S140;
00238 w *= a;
00239 goto S200;
00240 S140:
00241 w = infnty;
00242 goto S200;
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253 S160:
00254 if(z >= k2) goto S120;
00255 S170:
00256
00257
00258
00259
00260 v = beta*log(u1/(1.0-u1));
00261
00262 if(a > 1.0) goto S175;
00263
00264 if(v > expmax) goto S172;
00265 w = a*exp(v);
00266 goto S190;
00267 S172:
00268 w = v + log(a);
00269 if(w > expmax) goto S180;
00270 w = exp(w);
00271 goto S190;
00272 S175:
00273
00274 if(v > expmax) goto S180;
00275 w = exp(v);
00276 if(w > infnty/a) goto S180;
00277 w *= a;
00278 goto S190;
00279 S180:
00280 w = infnty;
00281
00282
00283
00284
00285
00286
00287
00288
00289 S190:
00290
00291
00292
00293
00294 if(alpha/(b+w) < minlog) goto S120;
00295 if(alpha*(log(alpha/(b+w))+v)-1.3862944 < log(z)) goto S120;
00296 S200:
00297
00298
00299
00300 if(!(a == aa)) goto S210;
00301 genbet = w/(b+w);
00302 goto S220;
00303 S210:
00304 genbet = b/(b+w);
00305 S230:
00306 S220:
00307 return genbet;
00308 }
00309
00310 BZ_NAMESPACE_END
00311
00312 #endif // BZ_RANDOM_BETA