boost/random/poisson_distribution.hpp

0001 /* boost random/poisson_distribution.hpp header file
0002  *
0003  * Copyright Jens Maurer 2002
0004  * Copyright Steven Watanabe 2010
0005  * Distributed under the Boost Software License, Version 1.0. (See
0006  * accompanying file LICENSE_1_0.txt or copy at
0007  * http://www.boost.org/LICENSE_1_0.txt)
0008  *
0009  * See http://www.boost.org for most recent version including documentation.
0010  *
0011  * $Id$
0012  *
0013  */
0014
0015 #ifndef BOOST_RANDOM_POISSON_DISTRIBUTION_HPP
0016 #define BOOST_RANDOM_POISSON_DISTRIBUTION_HPP
0017
0018 #include <boost/config/no_tr1/cmath.hpp>
0019 #include <cstdlib>
0020 #include <iosfwd>
0021 #include <boost/assert.hpp>
0022 #include <boost/limits.hpp>
0023 #include <boost/random/uniform_01.hpp>
0024 #include <boost/random/detail/config.hpp>
0025
0026 #include <boost/random/detail/disable_warnings.hpp>
0027
0028 namespace boost {
0029 namespace random {
0030
0031 namespace detail {
0032
0033 template<class RealType>
0034 struct poisson_table {
0035     static RealType value[10];
0036 };
0037
0038 template<class RealType>
0039 RealType poisson_table<RealType>::value[10] = {
0040     0.0,
0041     0.0,
0042     0.69314718055994529,
0043     1.7917594692280550,
0044     3.1780538303479458,
0045     4.7874917427820458,
0046     6.5792512120101012,
0047     8.5251613610654147,
0048     10.604602902745251,
0049     12.801827480081469
0050 };
0051
0052 }
0053
0054 /**
0055  * An instantiation of the class template @c poisson_distribution is a
0056  * model of \random_distribution.  The poisson distribution has
0057  * \f$p(i) = \frac{e^{-\lambda}\lambda^i}{i!}\f$
0058  *
0059  * This implementation is based on the PTRD algorithm described
0060  *
0061  *  @blockquote
0062  *  "The transformed rejection method for generating Poisson random variables",
0063  *  Wolfgang Hormann, Insurance: Mathematics and Economics
0064  *  Volume 12, Issue 1, February 1993, Pages 39-45
0065  *  @endblockquote
0066  */
0067 template<class IntType = int, class RealType = double>
0068 class poisson_distribution {
0069 public:
0070     typedef IntType result_type;
0071     typedef RealType input_type;
0072
0073     class param_type {
0074     public:
0075         typedef poisson_distribution distribution_type;
0076         /**
0077          * Construct a param_type object with the parameter "mean"
0078          *
0079          * Requires: mean > 0
0080          */
0081         explicit param_type(RealType mean_arg = RealType(1))
0082           : _mean(mean_arg)
0083         {
0084             BOOST_ASSERT(_mean > 0);
0085         }
0086         /* Returns the "mean" parameter of the distribution. */
0087         RealType mean() const { return _mean; }
0088 #ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
0089         /** Writes the parameters of the distribution to a @c std::ostream. */
0090         template<class CharT, class Traits>
0091         friend std::basic_ostream<CharT, Traits>&
0092         operator<<(std::basic_ostream<CharT, Traits>& os,
0093                    const param_type& parm)
0094         {
0095             os << parm._mean;
0096             return os;
0097         }
0098
0099         /** Reads the parameters of the distribution from a @c std::istream. */
0100         template<class CharT, class Traits>
0101         friend std::basic_istream<CharT, Traits>&
0102         operator>>(std::basic_istream<CharT, Traits>& is, param_type& parm)
0103         {
0104             is >> parm._mean;
0105             return is;
0106         }
0107 #endif
0108         /** Returns true if the parameters have the same values. */
0109         friend bool operator==(const param_type& lhs, const param_type& rhs)
0110         {
0111             return lhs._mean == rhs._mean;
0112         }
0113         /** Returns true if the parameters have different values. */
0114         friend bool operator!=(const param_type& lhs, const param_type& rhs)
0115         {
0116             return !(lhs == rhs);
0117         }
0118     private:
0119         RealType _mean;
0120     };
0121
0122     /**
0123      * Constructs a @c poisson_distribution with the parameter @c mean.
0124      *
0125      * Requires: mean > 0
0126      */
0127     explicit poisson_distribution(RealType mean_arg = RealType(1))
0128       : _mean(mean_arg)
0129     {
0130         BOOST_ASSERT(_mean > 0);
0131         init();
0132     }
0133
0134     /**
0135      * Construct an @c poisson_distribution object from the
0136      * parameters.
0137      */
0138     explicit poisson_distribution(const param_type& parm)
0139       : _mean(parm.mean())
0140     {
0141         init();
0142     }
0143
0144     /**
0145      * Returns a random variate distributed according to the
0146      * poisson distribution.
0147      */
0148     template<class URNG>
0149     IntType operator()(URNG& urng) const
0150     {
0151         if(use_inversion()) {
0152             return invert(urng);
0153         } else {
0154             return generate(urng);
0155         }
0156     }
0157
0158     /**
0159      * Returns a random variate distributed according to the
0160      * poisson distribution with parameters specified by param.
0161      */
0162     template<class URNG>
0163     IntType operator()(URNG& urng, const param_type& parm) const
0164     {
0165         return poisson_distribution(parm)(urng);
0166     }
0167
0168     /** Returns the "mean" parameter of the distribution. */
0169     RealType mean() const { return _mean; }
0170
0171     /** Returns the smallest value that the distribution can produce. */
0172     IntType min BOOST_PREVENT_MACRO_SUBSTITUTION() const { return 0; }
0173     /** Returns the largest value that the distribution can produce. */
0174     IntType max BOOST_PREVENT_MACRO_SUBSTITUTION() const
0175     { return (std::numeric_limits<IntType>::max)(); }
0176
0177     /** Returns the parameters of the distribution. */
0178     param_type param() const { return param_type(_mean); }
0179     /** Sets parameters of the distribution. */
0180     void param(const param_type& parm)
0181     {
0182         _mean = parm.mean();
0183         init();
0184     }
0185
0186     /**
0187      * Effects: Subsequent uses of the distribution do not depend
0188      * on values produced by any engine prior to invoking reset.
0189      */
0190     void reset() { }
0191
0192 #ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
0193     /** Writes the parameters of the distribution to a @c std::ostream. */
0194     template<class CharT, class Traits>
0195     friend std::basic_ostream<CharT,Traits>&
0196     operator<<(std::basic_ostream<CharT,Traits>& os,
0197                const poisson_distribution& pd)
0198     {
0199         os << pd.param();
0200         return os;
0201     }
0202
0203     /** Reads the parameters of the distribution from a @c std::istream. */
0204     template<class CharT, class Traits>
0205     friend std::basic_istream<CharT,Traits>&
0206     operator>>(std::basic_istream<CharT,Traits>& is, poisson_distribution& pd)
0207     {
0208         pd.read(is);
0209         return is;
0210     }
0211 #endif
0212
0213     /** Returns true if the two distributions will produce the same
0214         sequence of values, given equal generators. */
0215     friend bool operator==(const poisson_distribution& lhs,
0216                            const poisson_distribution& rhs)
0217     {
0218         return lhs._mean == rhs._mean;
0219     }
0220     /** Returns true if the two distributions could produce different
0221         sequences of values, given equal generators. */
0222     friend bool operator!=(const poisson_distribution& lhs,
0223                            const poisson_distribution& rhs)
0224     {
0225         return !(lhs == rhs);
0226     }
0227
0228 private:
0229
0230     /// @cond show_private
0231
0232     template<class CharT, class Traits>
0233     void read(std::basic_istream<CharT, Traits>& is) {
0234         param_type parm;
0235         if(is >> parm) {
0236             param(parm);
0237         }
0238     }
0239
0240     bool use_inversion() const
0241     {
0242         return _mean < 10;
0243     }
0244
0245     static RealType log_factorial(IntType k)
0246     {
0247         BOOST_ASSERT(k >= 0);
0248         BOOST_ASSERT(k < 10);
0249         return detail::poisson_table<RealType>::value[k];
0250     }
0251
0252     void init()
0253     {
0254         using std::sqrt;
0255         using std::exp;
0256
0257         if(use_inversion()) {
0258             _u._exp_mean = exp(-_mean);
0259         } else {
0260             _u._ptrd.smu = sqrt(_mean);
0261             _u._ptrd.b = 0.931 + 2.53 * _u._ptrd.smu;
0262             _u._ptrd.a = -0.059 + 0.02483 * _u._ptrd.b;
0263             _u._ptrd.inv_alpha = 1.1239 + 1.1328 / (_u._ptrd.b - 3.4);
0264             _u._ptrd.v_r = 0.9277 - 3.6224 / (_u._ptrd.b - 2);
0265         }
0266     }
0267
0268     template<class URNG>
0269     IntType generate(URNG& urng) const
0270     {
0271         using std::floor;
0272         using std::abs;
0273         using std::log;
0274
0275         while(true) {
0276             RealType u;
0277             RealType v = uniform_01<RealType>()(urng);
0278             if(v <= 0.86 * _u._ptrd.v_r) {
0279                 u = v / _u._ptrd.v_r - 0.43;
0280                 return static_cast<IntType>(floor(
0281                     (2*_u._ptrd.a/(0.5-abs(u)) + _u._ptrd.b)*u + _mean + 0.445));
0282             }
0283
0284             if(v >= _u._ptrd.v_r) {
0285                 u = uniform_01<RealType>()(urng) - 0.5;
0286             } else {
0287                 u = v/_u._ptrd.v_r - 0.93;
0288                 u = ((u < 0)? -0.5 : 0.5) - u;
0289                 v = uniform_01<RealType>()(urng) * _u._ptrd.v_r;
0290             }
0291
0292             RealType us = 0.5 - abs(u);
0293             if(us < 0.013 && v > us) {
0294                 continue;
0295             }
0296
0297             RealType k = floor((2*_u._ptrd.a/us + _u._ptrd.b)*u+_mean+0.445);
0298             v = v*_u._ptrd.inv_alpha/(_u._ptrd.a/(us*us) + _u._ptrd.b);
0299
0300             RealType log_sqrt_2pi = 0.91893853320467267;
0301
0302             if(k >= 10) {
0303                 if(log(v*_u._ptrd.smu) <= (k + 0.5)*log(_mean/k)
0304                                - _mean
0305                                - log_sqrt_2pi
0306                                + k
0307                                - (1/12. - (1/360. - 1/(1260.*k*k))/(k*k))/k) {
0308                     return static_cast<IntType>(k);
0309                 }
0310             } else if(k >= 0) {
0311                 if(log(v) <= k*log(_mean)
0312                            - _mean
0313                            - log_factorial(static_cast<IntType>(k))) {
0314                     return static_cast<IntType>(k);
0315                 }
0316             }
0317         }
0318     }
0319
0320     template<class URNG>
0321     IntType invert(URNG& urng) const
0322     {
0323         RealType p = _u._exp_mean;
0324         IntType x = 0;
0325         RealType u = uniform_01<RealType>()(urng);
0326         while(u > p) {
0327             u = u - p;
0328             ++x;
0329             p = _mean * p / x;
0330         }
0331         return x;
0332     }
0333
0334     RealType _mean;
0335
0336     union {
0337         // for ptrd
0338         struct {
0339             RealType v_r;
0340             RealType a;
0341             RealType b;
0342             RealType smu;
0343             RealType inv_alpha;
0344         } _ptrd;
0345         // for inversion
0346         RealType _exp_mean;
0347     } _u;
0348
0349     /// @endcond
0350 };
0351
0352 } // namespace random
0353
0354 using random::poisson_distribution;
0355
0356 } // namespace boost
0357
0358 #include <boost/random/detail/enable_warnings.hpp>
0359
0360 #endif // BOOST_RANDOM_POISSON_DISTRIBUTION_HPP