adonis wrote:
hello.
i need some help with a c or c++ program.
a program for generating random numbers with exponential
Take a variable X uniformly distributed in [0,1) and return
- a * log( 1.0 - X )
and poisson distribution.
That is much more tricky. The following is an implementation of the
poission_distribution from C++0X. You will not be able to use the code as
it is, but it should give you an idea as to how one can go about this.
// poisson_distribution.cc, Kai-Uwe Bux [2007-2008]
// ================================================
// Copyright notice
// ================
/*
This code is placed into the public domain by the author.
*/
// DISCLAIMER
// ==========
/*
This code may have undefined behavior and/or
fry your computer. Use at your own risk. You
have been warned.
*/
// Description
// ===========
/*
This implements poisson_distribution<from std.
*/
// Method
// ======
/*
J.H. Ahrens, U. Dieter: "Computer Generation of Poisson
Deviates from Modified Normal Distributions", ACM Transactions
on Mathematical Software 8 (1982) 163-179
BEWARE: there is a typo in Procedure F in the paper, it is
corrected below.
*/
#ifndef KUBUX_poisson_distribution
#define KUBUX_poisson_distribution
#include "contract.cc"
#include "constants.cc"
#include "polynomial.cc"
#include "uniform_real_distribution.cc"
#include "exponential_distribution.cc"
#include "normal_distribution.cc"
#include <cmath>
#include <math.h>
#include <vector>
#include <algorithm>
#include <limits>
namespace kubux {
using std::vector;
using std::upper_bound;
using std::abs;
using std::exp;
using std::log;
using std::pow;
using std::sqrt;
using std::numeric_limits;
template < typename IntType = int >
class poisson_distribution {
private:
typedef double float_type;
typedef uniform_real_distribution< float_type uniform_dist;
typedef exponential_distribution< float_type exponential_dist;
typedef normal_distribution< float_type normal_dist;
typedef typename uniform_dist::param_type uniform_param;
typedef typename normal_dist::param_type normal_param;
public:
typedef IntType result_type;
class param_type {
friend class poisson_distribution;
fhuge m;
// Case A:
fhuge s;
fhuge d;
fhuge L;
fhuge omega;
fhuge b1;
fhuge b2;
fhuge c3;
fhuge c2;
fhuge c1;
fhuge c0;
fhuge c;
mutable fhuge px;
mutable fhuge py;
mutable fhuge fx;
mutable fhuge fy;
void procedure_F ( result_type k ) const {
static fhuge a [10] = {
-0.5000000002L,
0.3333333343L,
-0.2499998565L,
0.1999997049L,
-0.1666848753L,
0.1428833286L,
-0.1241963125L,
0.1101687109L,
-0.1142650302L,
0.1055093006L };
step_1 :
if ( k < 10 ) {
px = -m;
py = pow( m, fhuge(k) ) / tgamma(fhuge(k+1));
} else {
fhuge delta = ( 1.0L / 12.0L ) / fhuge(k);
delta -= 4.8L * delta * delta * delta;
fhuge v = ( m - fhuge(k) ) / fhuge(k);
px = fhuge(k) * log( 1.0L + v ) - ( m - fhuge(k) ) - delta;
if ( -0.25L <= v && v <= 0.25L ) {
fhuge poly = eval_polynomial( a, a+10, v );
px = fhuge(k) * v*v * poly - delta;
} else {
px = fhuge(k) * log( 1.0L + v ) - ( m - fhuge(k) ) - delta;
}
py = ( 1.0L / sqrt2pi ) / sqrt(fhuge(k));
}
step_2 :
fhuge x = ( fhuge(k) - m + 0.5L ) / s;
x *= x;
// NOTE: [typo in the paper]
/*
The paper states the algorithm with
f_x <-- 0.5 X^2
However, this contradicts the (correct) formula (13).
The sign got lost in translation. We put it back in.
*/
fx = -0.5L * x;
fy = omega *
( ( ( c3 * x + c2 ) * x + c1 ) * x + c0 );
}
// Case B:
mutable fhuge p_last;
mutable vector< fhuge cumm;
result_type grow_table ( float_type uniform_variate ) const {
result_type k =
upper_bound( cumm.begin(), cumm.end(), uniform_variate )
-
cumm.begin();
if ( k == cumm.size() ) {
while ( cumm.size() <= 35 && cumm.back() < uniform_variate ) {
p_last *= m/float_type( cumm.size() );
cumm.push_back ( cumm.back() + p_last );
}
return ( cumm.size() - 1 );
}
return ( k );
}
template < typename UniformRNG >
result_type draw ( UniformRNG & urng,
uniform_dist & uni_d,
normal_dist & nor_d,
exponential_dist & exp_d ) const {
fhuge u;
fhuge g;
result_type k;
fhuge e;
fhuge t;
fhuge dummy;
if ( this->m < 10.0 ) {
return( this->grow_table( uni_d( urng ) ) );
} else {
step_N :
g = nor_d( urng, normal_param( this->m, this->s ) );
k = floor(g);
if ( g < 0.0 ) {
goto step_P;
}
step_I :
if ( fhuge(k) >= this->L ) {
return ( k );
}
step_S:
u = uni_d( urng );
dummy = this->m - fhuge(k);
if ( this->d * u >= dummy * dummy * dummy ) {
return ( k );
}
step_P :
if ( g >= 0.0 ) {
this->procedure_F( k );
} else {
goto step_E;
}
step_Q :
if ( this->fy * ( 1.0L - u )
<=
this->py * exp( this->px - this->fx ) ) {
return ( k );
}
step_E :
e = exp_d( urng );
u = uni_d( urng, uniform_param( -1.0L, 1.0L ) );
t = 1.8L + ( u 0 ? e : -e );
if ( t <= -0.6744L ) {
goto step_E;
}
k = floor( this->m + this->s * t );
this->procedure_F( k );
step_H :
if ( this->c * abs(u)
>
this->py * exp( this->px + e )
-
this->fy * exp( this->fx + e ) ) {
goto step_E;
}
return ( k );
}
}
public:
typedef poisson_distribution distribution_type;
param_type ( float_type m_par )
: m ( m_par )
, s ( sqrt( m ) )
, d ( 6.0L * m*m )
, L ( floor( m - 1.1484L ) )
, omega ( ( 1.0L / sqrt2pi ) / s )
, b1 ( ( 1.0L / 24.0L ) / m )
, b2 ( 0.3L * b1*b1 )
, c3 ( ( 1.0L / 7.0L ) * b1 * b2 )
, c2 ( b2 - 15.0L * c3 )
, c1 ( b1 - 6.0L * b2 + 45.0L * c3 )
, c0 ( 1.0L - b1 + 3.0L * b2 - 15.0L * c3 )
, c ( 0.1069L / m )
, p_last ( exp( -m ) )
, cumm ()
{
cumm.push_back( p_last );
KUBUX_ENFORCE( float_type(0) < m );
}
result_type mean ( void ) const {
return ( m );
}
};
private:
param_type the_parm;
uniform_dist uni_d;
normal_dist nor_d;
exponential_dist exp_d;
public:
explicit
poisson_distribution ( float_type mean = 1.0 )
: the_parm( mean )
, uni_d ( 0, 1 )
, nor_d ( 0, 1 )
, exp_d ( 1 )
{}
explicit
poisson_distribution ( param_type const & parm )
: the_parm ( parm )
, uni_d ( 0, 1 )
, nor_d ( 0, 1 )
, exp_d ( 1 )
{}
param_type param ( void ) const {
return ( the_parm );
}
void param ( param_type const & parm ) {
the_parm = parm;
}
float_type mean ( void ) const {
return ( the_parm.mean() );
}
result_type min ( void ) const {
return ( 0 );
}
result_type max ( void ) const {
return ( numeric_limits<result_type>::max() );
}
template < typename UniformRNG >
result_type operator() ( UniformRNG & urng,
param_type const & parm ) {
return( parm.draw( urng, uni_d, nor_d, exp_d ) );
}
template < typename UniformRNG >
result_type operator() ( UniformRNG & urng ) {
return ( this->operator()( urng, the_parm ) );
}
void reset ( void ) {
uni_d.reset();
nor_d.reset();
exp_d.reset();
}
};
} // namespace kubux
#endif // KUBUX_poisson_distribution
// end of file
Best
Kai-Uwe Bux
