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