Anthony's VBA ForumWhere the knowledge is shared

 Who's Online [ Administrator ][ Moderator ]

Author Message
Lim

Joined: 27 Oct 2004
Posts: 7

 Posted: Sun Dec 12, 2004 10:53 pm    Post subject: Beta Distribution Anyone know how to generate beta distribution random number?
Anthony

Joined: 29 Sep 2004
Posts: 36
Location: New York, USA

 Posted: Mon Dec 13, 2004 2:34 pm    Post subject: Do you know this language? Hi lim, I got these codes off the Internet. It seems complex. It is probably written in C or Fortran. As it claims, it should generate the Beta deviate. Anthony ====================================== #include "randlib.h" #include #include #include #define ABS(x) ((x) >= 0 ? (x) : -(x)) #define min(a,b) ((a) <= (b) ? (a) : (b)) #define max(a,b) ((a) >= (b) ? (a) : (b)) void ftnstop2(char*); float genbet(float aa,float bb) /* ********************************************************************** float genbet(float aa,float bb) GeNerate BETa random deviate Function Returns a single random deviate from the beta distribution with parameters A and B. The density of the beta is x^(a-1) * (1-x)^(b-1) / B(a,b) for 0 < x < 1 Arguments aa --> First parameter of the beta distribution bb --> Second parameter of the beta distribution Method R. C. H. Cheng Generating Beta Variatew with Nonintegral Shape Parameters Communications of the ACM, 21:317-322 (1978) (Algorithms BB and BC) ********************************************************************** */ { /* JJV changed expmax (log(1.0E38)==87.49823), and added minlog */ #define expmax 87.49823 #define infnty 1.0E38 #define minlog 1.0E-37 float olda = -1.0E37; float oldb = -1.0E37; float genbet,a,alpha,b,beta,delta,gamma,k1,k2,r,s,t,u1,u2,v,w,y,z; long qsame; qsame = olda == aa && oldb == bb; if(qsame) goto S20; if(!(aa < minlog || bb < minlog)) goto S10; fputs(" AA or BB < 1.0E-37 in GENBET - Abort!\n",stderr); fprintf(stderr," AA: %16.6E BB %16.6E\n",aa,bb); exit(1); S10: olda = aa; oldb = bb; S20: if(!(min(aa,bb) > 1.0)) goto S100; /* Alborithm BB Initialize */ if(qsame) goto S30; a = min(aa,bb); b = max(aa,bb); alpha = a+b; beta = sqrt((alpha-2.0)/(2.0*a*b-alpha)); gamma = a+1.0/beta; S30: S40: u1 = ranf(); /* Step 1 */ u2 = ranf(); v = beta*log(u1/(1.0-u1)); /* JJV altered this */ if(v > expmax) goto S55; /* * JJV added checker to see if a*exp(v) will overflow * JJV S50 _was_ w = a*exp(v); also note here a > 1.0 */ w = exp(v); if(w > infnty/a) goto S55; w *= a; goto S60; S55: w = infnty; S60: z = pow(u1,2.0)*u2; r = gamma*v-1.3862944; s = a+r-w; /* Step 2 */ if(s+2.609438 >= 5.0*z) goto S70; /* Step 3 */ t = log(z); if(s > t) goto S70; /* * Step 4 * * JJV added checker to see if log(alpha/(b+w)) will * JJV overflow. If so, we count the log as -INF, and * JJV consequently evaluate conditional as true, i.e. * JJV the algorithm rejects the trial and starts over * JJV May not need this here since alpha > 2.0 */ if(alpha/(b+w) < minlog) goto S40; if(r+alpha*log(alpha/(b+w)) < t) goto S40; S70: /* Step 5 */ if(!(aa == a)) goto S80; genbet = w/(b+w); goto S90; S80: genbet = b/(b+w); S90: goto S230; S100: /* Algorithm BC Initialize */ if(qsame) goto S110; a = max(aa,bb); b = min(aa,bb); alpha = a+b; beta = 1.0/b; delta = 1.0+a-b; k1 = delta*(1.38889E-2+4.16667E-2*b)/(a*beta-0.777778); k2 = 0.25+(0.5+0.25/delta)*b; S110: S120: u1 = ranf(); /* Step 1 */ u2 = ranf(); if(u1 >= 0.5) goto S130; /* Step 2 */ y = u1*u2; z = u1*y; if(0.25*u2+z-y >= k1) goto S120; goto S170; S130: /* Step 3 */ z = pow(u1,2.0)*u2; if(!(z <= 0.25)) goto S160; v = beta*log(u1/(1.0-u1)); /* * JJV instead of checking v > expmax at top, I will check * JJV if a < 1, then check the appropriate values */ if(a > 1.0) goto S135; /* JJV a < 1 so it can help out if exp(v) would overflow */ if(v > expmax) goto S132; w = a*exp(v); goto S200; S132: w = v + log(a); if(w > expmax) goto S140; w = exp(w); goto S200; S135: /* JJV in this case a > 1 */ if(v > expmax) goto S140; w = exp(v); if(w > infnty/a) goto S140; w *= a; goto S200; S140: w = infnty; goto S200; /* * JJV old code * if(!(v > expmax)) goto S140; * w = infnty; * goto S150; *S140: * w = a*exp(v); *S150: * goto S200; */ S160: if(z >= k2) goto S120; S170: /* Step 4 Step 5 */ v = beta*log(u1/(1.0-u1)); /* JJV same kind of checking as above */ if(a > 1.0) goto S175; /* JJV a < 1 so it can help out if exp(v) would overflow */ if(v > expmax) goto S172; w = a*exp(v); goto S190; S172: w = v + log(a); if(w > expmax) goto S180; w = exp(w); goto S190; S175: /* JJV in this case a > 1.0 */ if(v > expmax) goto S180; w = exp(v); if(w > infnty/a) goto S180; w *= a; goto S190; S180: w = infnty; /* * JJV old code * if(!(v > expmax)) goto S180; * w = infnty; * goto S190; *S180: * w = a*exp(v); */ S190: /* * JJV here we also check to see if log overlows; if so, we treat it * JJV as -INF, which means condition is true, i.e. restart */ if(alpha/(b+w) < minlog) goto S120; if(alpha*(log(alpha/(b+w))+v)-1.3862944 < log(z)) goto S120; S200: /* Step 6 */ if(!(a == aa)) goto S210; genbet = w/(b+w); goto S220; S210: genbet = b/(b+w); S230: S220: return genbet; #undef expmax #undef infnty #undef minlog } float genchi(float df) /*
 Display posts from previous: All Posts1 Day7 Days2 Weeks1 Month3 Months6 Months1 Year Oldest FirstNewest First
 All times are GMT - 4 Hours Page 1 of 1

 Jump to: Select a forum Excel VBA----------------General Excel VBA VBA on Finance and Investment----------------Option Pricing ModelsPortfoilio OptimizationOther Finance Topics VBA on Statistics----------------Probability DistributionsMonte Carlo SimulationRegression Analysis Announcement----------------Updated News
You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot vote in polls in this forum