Poisson noise

Boost::randomを利用してPoisson乱数を生成する関数

#include <boost/random.hpp>
#include <ctime>

int gen_poisson_rand(double mean)
{
	static boost::mt19937 gen(static_cast<unsigned long>(time(0)));
	return boost::variate_generator
		<boost::mt19937&, boost::poisson_distribution<> >
			(gen, boost::poisson_distribution<>(mean))();
}

ただし、現状では高速な生成アルゴリズムは実装されていないっぽい。

  template<class Engine>
  result_type operator()(Engine& eng)
  {
    // TODO: This is O(_mean), but it should be O(log(_mean)) for large _mean
    RealType product = RealType(1);
    for(result_type m = 0; ; ++m) {
      product *= eng();
      if(product <= _exp_mean)
        return m;
    }
  }

Numerical Recipesには高速化手法が載ってるけどライセンス的に問題があるから使いにくいんだろうなあ。