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