Rejection loop in C
- From: Darius <my_vision2@xxxxxxxxx>
- Date: Thu, 31 Aug 2006 10:30:22 EDT
I understand what the below codes doin except that m_MaxDensityVal, DistributionDensity, are they pre-defined functions? is there any formula for them cause I need to define them in Pascal. Any help would be apprecited.
-----------
The sampling range of the needed distribution is [a,b] and the probability
density function (pdf) is denoted by f(x).
x1 = a + r1 * (b -a );
if( r2 * Maximum_of_PDF <= f( x1 ) )
then accept x1
else
throw away r1 & r2 and start again.
In principle a random number is accepted more often if the PDF for this value is bigger. Thus the desired distribution is sampled.
This is the so called rejection loop
Here is some code to illustrate the algorithm:
const double PI = 4.0 * atan(1.0);
/////////////////////////////////////////////////////////
double CRejectionMethod::DistributionDensity( double x )
// The implementation of the distribution density function.
// In this case its a normal distribution!
/////////////////////////////////////////////////////////{
double RetVal;
RetVal = 1 / ( sqrt( 2 * PI ) ) * exp( - (x * x) / 2. );
return RetVal;
}
/////////////////////////////////////////////////////////
double CRejectionMethod::GetRandomNr( double& a, double& b )
// This is the main routine, that implements the specific
// algorithm to deliver the random numbers between [a,b]
/////////////////////////////////////////////////////////{
double R1, R2;
double x;
int bReject = true;
m_MaxDensityVal = FindMaxDensityVal(); // for a normal distribution
this would be m_MaxDensityVal = DistributionDensity( 0.0 );
// main algorithm, the rejection loop
while( bReject == true ) {
R1 = GetlRandomNr();
x = a + (b - a) * R1;
R2 = GetlRandomNr();
if( (R2 * m_MaxDensityVal) <= (DistributionDensity(x)) )
return x;
}
return false;
}
.
- Prev by Date: Re: A statistical list problem
- Next by Date: Re: A statistical list problem
- Previous by thread: A statistical list problem
- Next by thread: Multinomial approximation to Poisson ??
- Index(es):
Relevant Pages
|
|
