/* boost random/poisson_distribution.hpp header file * * Copyright Jens Maurer 2002 * Copyright Steven Watanabe 2010 * Distributed under the Boost Software License, Version 1.0. (See * accompanying file LICENSE_1_0.txt or copy at * http://www.boost.org/LICENSE_1_0.txt) * * See http://www.boost.org for most recent version including documentation. * * $Id$ * */ #ifndef BOOST_RANDOM_POISSON_DISTRIBUTION_HPP #define BOOST_RANDOM_POISSON_DISTRIBUTION_HPP #include <boost/config/no_tr1/cmath.hpp> #include <cstdlib> #include <iosfwd> #include <boost/assert.hpp> #include <boost/limits.hpp> #include <boost/random/uniform_01.hpp> #include <boost/random/detail/config.hpp> #include <boost/random/detail/disable_warnings.hpp> namespace boost { namespace random { namespace detail { template<class RealType> struct poisson_table { static RealType value[10]; }; template<class RealType> RealType poisson_table<RealType>::value[10] = { 0.0, 0.0, 0.69314718055994529, 1.7917594692280550, 3.1780538303479458, 4.7874917427820458, 6.5792512120101012, 8.5251613610654147, 10.604602902745251, 12.801827480081469 }; } /** * An instantiation of the class template @c poisson_distribution is a * model of \random_distribution. The poisson distribution has * \f$p(i) = \frac{e^{-\lambda}\lambda^i}{i!}\f$ * * This implementation is based on the PTRD algorithm described * * @blockquote * "The transformed rejection method for generating Poisson random variables", * Wolfgang Hormann, Insurance: Mathematics and Economics * Volume 12, Issue 1, February 1993, Pages 39-45 * @endblockquote */ template<class IntType = int, class RealType = double> class poisson_distribution { public: typedef IntType result_type; typedef RealType input_type; class param_type { public: typedef poisson_distribution distribution_type; /** * Construct a param_type object with the parameter "mean" * * Requires: mean > 0 */ explicit param_type(RealType mean_arg = RealType(1)) : _mean(mean_arg) { BOOST_ASSERT(_mean > 0); } /* Returns the "mean" parameter of the distribution. */ RealType mean() const { return _mean; } #ifndef BOOST_RANDOM_NO_STREAM_OPERATORS /** Writes the parameters of the distribution to a @c std::ostream. */ template<class CharT, class Traits> friend std::basic_ostream<CharT, Traits>& operator<<(std::basic_ostream<CharT, Traits>& os, const param_type& parm) { os << parm._mean; return os; } /** Reads the parameters of the distribution from a @c std::istream. */ template<class CharT, class Traits> friend std::basic_istream<CharT, Traits>& operator>>(std::basic_istream<CharT, Traits>& is, param_type& parm) { is >> parm._mean; return is; } #endif /** Returns true if the parameters have the same values. */ friend bool operator==(const param_type& lhs, const param_type& rhs) { return lhs._mean == rhs._mean; } /** Returns true if the parameters have different values. */ friend bool operator!=(const param_type& lhs, const param_type& rhs) { return !(lhs == rhs); } private: RealType _mean; }; /** * Constructs a @c poisson_distribution with the parameter @c mean. * * Requires: mean > 0 */ explicit poisson_distribution(RealType mean_arg = RealType(1)) : _mean(mean_arg) { BOOST_ASSERT(_mean > 0); init(); } /** * Construct an @c poisson_distribution object from the * parameters. */ explicit poisson_distribution(const param_type& parm) : _mean(parm.mean()) { init(); } /** * Returns a random variate distributed according to the * poisson distribution. */ template<class URNG> IntType operator()(URNG& urng) const { if(use_inversion()) { return invert(urng); } else { return generate(urng); } } /** * Returns a random variate distributed according to the * poisson distribution with parameters specified by param. */ template<class URNG> IntType operator()(URNG& urng, const param_type& parm) const { return poisson_distribution(parm)(urng); } /** Returns the "mean" parameter of the distribution. */ RealType mean() const { return _mean; } /** Returns the smallest value that the distribution can produce. */ IntType min BOOST_PREVENT_MACRO_SUBSTITUTION() const { return 0; } /** Returns the largest value that the distribution can produce. */ IntType max BOOST_PREVENT_MACRO_SUBSTITUTION() const { return (std::numeric_limits<IntType>::max)(); } /** Returns the parameters of the distribution. */ param_type param() const { return param_type(_mean); } /** Sets parameters of the distribution. */ void param(const param_type& parm) { _mean = parm.mean(); init(); } /** * Effects: Subsequent uses of the distribution do not depend * on values produced by any engine prior to invoking reset. */ void reset() { } #ifndef BOOST_RANDOM_NO_STREAM_OPERATORS /** Writes the parameters of the distribution to a @c std::ostream. */ template<class CharT, class Traits> friend std::basic_ostream<CharT,Traits>& operator<<(std::basic_ostream<CharT,Traits>& os, const poisson_distribution& pd) { os << pd.param(); return os; } /** Reads the parameters of the distribution from a @c std::istream. */ template<class CharT, class Traits> friend std::basic_istream<CharT,Traits>& operator>>(std::basic_istream<CharT,Traits>& is, poisson_distribution& pd) { pd.read(is); return is; } #endif /** Returns true if the two distributions will produce the same sequence of values, given equal generators. */ friend bool operator==(const poisson_distribution& lhs, const poisson_distribution& rhs) { return lhs._mean == rhs._mean; } /** Returns true if the two distributions could produce different sequences of values, given equal generators. */ friend bool operator!=(const poisson_distribution& lhs, const poisson_distribution& rhs) { return !(lhs == rhs); } private: /// @cond show_private template<class CharT, class Traits> void read(std::basic_istream<CharT, Traits>& is) { param_type parm; if(is >> parm) { param(parm); } } bool use_inversion() const { return _mean < 10; } static RealType log_factorial(IntType k) { BOOST_ASSERT(k >= 0); BOOST_ASSERT(k < 10); return detail::poisson_table<RealType>::value[k]; } void init() { using std::sqrt; using std::exp; if(use_inversion()) { _exp_mean = exp(-_mean); } else { _ptrd.smu = sqrt(_mean); _ptrd.b = 0.931 + 2.53 * _ptrd.smu; _ptrd.a = -0.059 + 0.02483 * _ptrd.b; _ptrd.inv_alpha = 1.1239 + 1.1328 / (_ptrd.b - 3.4); _ptrd.v_r = 0.9277 - 3.6224 / (_ptrd.b - 2); } } template<class URNG> IntType generate(URNG& urng) const { using std::floor; using std::abs; using std::log; while(true) { RealType u; RealType v = uniform_01<RealType>()(urng); if(v <= 0.86 * _ptrd.v_r) { u = v / _ptrd.v_r - 0.43; return static_cast<IntType>(floor( (2*_ptrd.a/(0.5-abs(u)) + _ptrd.b)*u + _mean + 0.445)); } if(v >= _ptrd.v_r) { u = uniform_01<RealType>()(urng) - 0.5; } else { u = v/_ptrd.v_r - 0.93; u = ((u < 0)? -0.5 : 0.5) - u; v = uniform_01<RealType>()(urng) * _ptrd.v_r; } RealType us = 0.5 - abs(u); if(us < 0.013 && v > us) { continue; } RealType k = floor((2*_ptrd.a/us + _ptrd.b)*u+_mean+0.445); v = v*_ptrd.inv_alpha/(_ptrd.a/(us*us) + _ptrd.b); RealType log_sqrt_2pi = 0.91893853320467267; if(k >= 10) { if(log(v*_ptrd.smu) <= (k + 0.5)*log(_mean/k) - _mean - log_sqrt_2pi + k - (1/12. - (1/360. - 1/(1260.*k*k))/(k*k))/k) { return static_cast<IntType>(k); } } else if(k >= 0) { if(log(v) <= k*log(_mean) - _mean - log_factorial(static_cast<IntType>(k))) { return static_cast<IntType>(k); } } } } template<class URNG> IntType invert(URNG& urng) const { RealType p = _exp_mean; IntType x = 0; RealType u = uniform_01<RealType>()(urng); while(u > p) { u = u - p; ++x; p = _mean * p / x; } return x; } RealType _mean; union { // for ptrd struct { RealType v_r; RealType a; RealType b; RealType smu; RealType inv_alpha; } _ptrd; // for inversion RealType _exp_mean; }; /// @endcond }; } // namespace random using random::poisson_distribution; } // namespace boost #include <boost/random/detail/enable_warnings.hpp> #endif // BOOST_RANDOM_POISSON_DISTRIBUTION_HPP