Anthony's VBA Forum Index Anthony's VBA Forum
Where the knowledge is shared
Menu
 Anthony's VBA Forum IndexHome Page
 Anthony's VBA Forum IndexForum Index
FAQFAQ
MemberlistMemberlist
UsergroupsUsergroups
RegisterRegister
ProfileProfile
Log in to check your private messagesMessages
Log inLogin/Out

Quick Search

Advanced Search

Links
Consulting
Products

Who's Online
[ Administrator ]
[ Moderator ]


Beta Distribution

 
Post new topic   Reply to topic     Anthony's VBA Forum Index -> Probability Distributions
View previous topic :: View next topic  
Author Message
Lim



Joined: 27 Oct 2004
Posts: 7

PostPosted: Sun Dec 12, 2004 10:53 pm    Post subject: Beta Distribution Reply with quote

Anyone know how to generate beta distribution random number? Question Neutral
Back to top
View user's profile Send private message Send e-mail
Anthony
Site Admin
Site Admin


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

PostPosted: Mon Dec 13, 2004 2:34 pm    Post subject: Do you know this language? Reply with quote

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 <stdio.h>
#include <math.h>
#include <stdlib.h>
#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)
/*
Back to top
View user's profile Send private message
Display posts from previous:   
Post new topic   Reply to topic     Anthony's VBA Forum Index -> Probability Distributions All times are GMT - 4 Hours
Page 1 of 1

 
Jump to:  
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


Powered by phpBB © 2001, 2002 phpBB Group