I implemented this function to generate a poisson random variable
typedef long unsigned int luint;
luint poisson(luint lambda) {
double L = exp(-double(lambda));
luint k = 0;
double p = 1;
do {
k++;
p *= mrand.rand();
} while( p > L);
return (k-1);
}
where mrand is the MersenneTwister random number generator. I find that, as I increase lambda, the expected distribution is going to be wrong, with a mean that saturates at around 750. Is it due to numerical approximations or did I make any mistakes?
If you go the "existing library" route, your compiler may already support the C++11 std::random package. Here is how you use it:
#include <random>
#include <ctime>
#include <iostream>
std::mt19937 mrand(std::time(0)); // seed however you want
typedef long unsigned int luint;
luint poisson(luint lambda)
{
std::poisson_distribution<luint> d(lambda);
return d(mrand);
}
int main()
{
std::cout << poisson(750) << '\n';
std::poisson_distribution<luint> d(750);
std::cout << d(mrand) << '\n';
std::cout << d(mrand) << '\n';
}
I've used it two ways above:
I tried to imitate your existing interface.
If you create a std::poisson_distribution with a mean, it is more efficient to use that distribution over and over for the same mean (as done in main()).
Here is sample output for me:
751
730
779