boost/random/poisson_distribution.hpp

   1 /* boost random/poisson_distribution.hpp header file
   2  *
   3  * Copyright Jens Maurer 2002
   4  * Copyright Steven Watanabe 2010
   5  * Distributed under the Boost Software License, Version 1.0. (See
   6  * accompanying file LICENSE_1_0.txt or copy at
   7  * http://www.boost.org/LICENSE_1_0.txt)
   8  *
   9  * See http://www.boost.org for most recent version including documentation.
  10  *
  11  * $Id: poisson_distribution.hpp 71018 2011-04-05 21:27:52Z steven_watanabe $
  12  *
  13  */
  14
  15 #ifndef BOOST_RANDOM_POISSON_DISTRIBUTION_HPP
  16 #define BOOST_RANDOM_POISSON_DISTRIBUTION_HPP
  17
  18 #include <boost/config/no_tr1/cmath.hpp>
  19 #include <cstdlib>
  20 #include <iosfwd>
  21 #include <boost/assert.hpp>
  22 #include <boost/limits.hpp>
  23 #include <boost/random/uniform_01.hpp>
  24 #include <boost/random/detail/config.hpp>
  25
  26 #include <boost/random/detail/disable_warnings.hpp>
  27
  28 namespace boost {
  29 namespace random {
  30
  31 namespace detail {
  32
  33 template<class RealType>
  34 struct poisson_table {
  35     static RealType value[10];
  36 };
  37
  38 template<class RealType>
  39 RealType poisson_table<RealType>::value[10] = {
  40     0.0,
  41     0.0,
  42     0.69314718055994529,
  43     1.7917594692280550,
  44     3.1780538303479458,
  45     4.7874917427820458,
  46     6.5792512120101012,
  47     8.5251613610654147,
  48     10.604602902745251,
  49     12.801827480081469
  50 };
  51
  52 }
  53
  54 /**
  55  * An instantiation of the class template @c poisson_distribution is a
  56  * model of \random_distribution.  The poisson distribution has
  57  * \f$p(i) = \frac{e^{-\lambda}\lambda^i}{i!}\f$
  58  *
  59  * This implementation is based on the PTRD algorithm described
  60  *
  61  *  @blockquote
  62  *  "The transformed rejection method for generating Poisson random variables",
  63  *  Wolfgang Hormann, Insurance: Mathematics and Economics
  64  *  Volume 12, Issue 1, February 1993, Pages 39-45
  65  *  @endblockquote
  66  */
  67 template<class IntType = int, class RealType = double>
  68 class poisson_distribution {
  69 public:
  70     typedef IntType result_type;
  71     typedef RealType input_type;
  72
  73     class param_type {
  74     public:
  75         typedef poisson_distribution distribution_type;
  76         /**
  77          * Construct a param_type object with the parameter "mean"
  78          *
  79          * Requires: mean > 0
  80          */
  81         explicit param_type(RealType mean_arg = RealType(1))
  82           : _mean(mean_arg)
  83         {
  84             BOOST_ASSERT(_mean > 0);
  85         }
  86         /* Returns the "mean" parameter of the distribution. */
  87         RealType mean() const { return _mean; }
  88 #ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
  89         /** Writes the parameters of the distribution to a @c std::ostream. */
  90         template<class CharT, class Traits>
  91         friend std::basic_ostream<CharT, Traits>&
  92         operator<<(std::basic_ostream<CharT, Traits>& os,
  93                    const param_type& parm)
  94         {
  95             os << parm._mean;
  96             return os;
  97         }
  98
  99         /** Reads the parameters of the distribution from a @c std::istream. */
 100         template<class CharT, class Traits>
 101         friend std::basic_istream<CharT, Traits>&
 102         operator>>(std::basic_istream<CharT, Traits>& is, param_type& parm)
 103         {
 104             is >> parm._mean;
 105             return is;
 106         }
 107 #endif
 108         /** Returns true if the parameters have the same values. */
 109         friend bool operator==(const param_type& lhs, const param_type& rhs)
 110         {
 111             return lhs._mean == rhs._mean;
 112         }
 113         /** Returns true if the parameters have different values. */
 114         friend bool operator!=(const param_type& lhs, const param_type& rhs)
 115         {
 116             return !(lhs == rhs);
 117         }
 118     private:
 119         RealType _mean;
 120     };
 121
 122     /**
 123      * Constructs a @c poisson_distribution with the parameter @c mean.
 124      *
 125      * Requires: mean > 0
 126      */
 127     explicit poisson_distribution(RealType mean_arg = RealType(1))
 128       : _mean(mean_arg)
 129     {
 130         BOOST_ASSERT(_mean > 0);
 131         init();
 132     }
 133
 134     /**
 135      * Construct an @c poisson_distribution object from the
 136      * parameters.
 137      */
 138     explicit poisson_distribution(const param_type& parm)
 139       : _mean(parm.mean())
 140     {
 141         init();
 142     }
 143
 144     /**
 145      * Returns a random variate distributed according to the
 146      * poisson distribution.
 147      */
 148     template<class URNG>
 149     IntType operator()(URNG& urng) const
 150     {
 151         if(use_inversion()) {
 152             return invert(urng);
 153         } else {
 154             return generate(urng);
 155         }
 156     }
 157
 158     /**
 159      * Returns a random variate distributed according to the
 160      * poisson distribution with parameters specified by param.
 161      */
 162     template<class URNG>
 163     IntType operator()(URNG& urng, const param_type& parm) const
 164     {
 165         return poisson_distribution(parm)(urng);
 166     }
 167
 168     /** Returns the "mean" parameter of the distribution. */
 169     RealType mean() const { return _mean; }
 170
 171     /** Returns the smallest value that the distribution can produce. */
 172     IntType min BOOST_PREVENT_MACRO_SUBSTITUTION() const { return 0; }
 173     /** Returns the largest value that the distribution can produce. */
 174     IntType max BOOST_PREVENT_MACRO_SUBSTITUTION() const
 175     { return (std::numeric_limits<IntType>::max)(); }
 176
 177     /** Returns the parameters of the distribution. */
 178     param_type param() const { return param_type(_mean); }
 179     /** Sets parameters of the distribution. */
 180     void param(const param_type& parm)
 181     {
 182         _mean = parm.mean();
 183         init();
 184     }
 185
 186     /**
 187      * Effects: Subsequent uses of the distribution do not depend
 188      * on values produced by any engine prior to invoking reset.
 189      */
 190     void reset() { }
 191
 192 #ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
 193     /** Writes the parameters of the distribution to a @c std::ostream. */
 194     template<class CharT, class Traits>
 195     friend std::basic_ostream<CharT,Traits>&
 196     operator<<(std::basic_ostream<CharT,Traits>& os,
 197                const poisson_distribution& pd)
 198     {
 199         os << pd.param();
 200         return os;
 201     }
 202
 203     /** Reads the parameters of the distribution from a @c std::istream. */
 204     template<class CharT, class Traits>
 205     friend std::basic_istream<CharT,Traits>&
 206     operator>>(std::basic_istream<CharT,Traits>& is, poisson_distribution& pd)
 207     {
 208         pd.read(is);
 209         return is;
 210     }
 211 #endif
 212
 213     /** Returns true if the two distributions will produce the same
 214         sequence of values, given equal generators. */
 215     friend bool operator==(const poisson_distribution& lhs,
 216                            const poisson_distribution& rhs)
 217     {
 218         return lhs._mean == rhs._mean;
 219     }
 220     /** Returns true if the two distributions could produce different
 221         sequences of values, given equal generators. */
 222     friend bool operator!=(const poisson_distribution& lhs,
 223                            const poisson_distribution& rhs)
 224     {
 225         return !(lhs == rhs);
 226     }
 227
 228 private:
 229
 230     /// @cond show_private
 231
 232     template<class CharT, class Traits>
 233     void read(std::basic_istream<CharT, Traits>& is) {
 234         param_type parm;
 235         if(is >> parm) {
 236             param(parm);
 237         }
 238     }
 239
 240     bool use_inversion() const
 241     {
 242         return _mean < 10;
 243     }
 244
 245     static RealType log_factorial(IntType k)
 246     {
 247         BOOST_ASSERT(k >= 0);
 248         BOOST_ASSERT(k < 10);
 249         return detail::poisson_table<RealType>::value[k];
 250     }
 251
 252     void init()
 253     {
 254         using std::sqrt;
 255         using std::exp;
 256
 257         if(use_inversion()) {
 258             _exp_mean = exp(-_mean);
 259         } else {
 260             _ptrd.smu = sqrt(_mean);
 261             _ptrd.b = 0.931 + 2.53 * _ptrd.smu;
 262             _ptrd.a = -0.059 + 0.02483 * _ptrd.b;
 263             _ptrd.inv_alpha = 1.1239 + 1.1328 / (_ptrd.b - 3.4);
 264             _ptrd.v_r = 0.9277 - 3.6224 / (_ptrd.b - 2);
 265         }
 266     }
 267
 268     template<class URNG>
 269     IntType generate(URNG& urng) const
 270     {
 271         using std::floor;
 272         using std::abs;
 273         using std::log;
 274
 275         while(true) {
 276             RealType u;
 277             RealType v = uniform_01<RealType>()(urng);
 278             if(v <= 0.86 * _ptrd.v_r) {
 279                 u = v / _ptrd.v_r - 0.43;
 280                 return static_cast<IntType>(floor(
 281                     (2*_ptrd.a/(0.5-abs(u)) + _ptrd.b)*u + _mean + 0.445));
 282             }
 283
 284             if(v >= _ptrd.v_r) {
 285                 u = uniform_01<RealType>()(urng) - 0.5;
 286             } else {
 287                 u = v/_ptrd.v_r - 0.93;
 288                 u = ((u < 0)? -0.5 : 0.5) - u;
 289                 v = uniform_01<RealType>()(urng) * _ptrd.v_r;
 290             }
 291
 292             RealType us = 0.5 - abs(u);
 293             if(us < 0.013 && v > us) {
 294                 continue;
 295             }
 296
 297             RealType k = floor((2*_ptrd.a/us + _ptrd.b)*u+_mean+0.445);
 298             v = v*_ptrd.inv_alpha/(_ptrd.a/(us*us) + _ptrd.b);
 299
 300             RealType log_sqrt_2pi = 0.91893853320467267;
 301
 302             if(k >= 10) {
 303                 if(log(v*_ptrd.smu) <= (k + 0.5)*log(_mean/k)
 304                                - _mean
 305                                - log_sqrt_2pi
 306                                + k
 307                                - (1/12. - (1/360. - 1/(1260.*k*k))/(k*k))/k) {
 308                     return static_cast<IntType>(k);
 309                 }
 310             } else if(k >= 0) {
 311                 if(log(v) <= k*log(_mean)
 312                            - _mean
 313                            - log_factorial(static_cast<IntType>(k))) {
 314                     return static_cast<IntType>(k);
 315                 }
 316             }
 317         }
 318     }
 319
 320     template<class URNG>
 321     IntType invert(URNG& urng) const
 322     {
 323         RealType p = _exp_mean;
 324         IntType x = 0;
 325         RealType u = uniform_01<RealType>()(urng);
 326         while(u > p) {
 327             u = u - p;
 328             ++x;
 329             p = _mean * p / x;
 330         }
 331         return x;
 332     }
 333
 334     RealType _mean;
 335
 336     union {
 337         // for ptrd
 338         struct {
 339             RealType v_r;
 340             RealType a;
 341             RealType b;
 342             RealType smu;
 343             RealType inv_alpha;
 344         } _ptrd;
 345         // for inversion
 346         RealType _exp_mean;
 347     };
 348
 349     /// @endcond
 350 };
 351
 352 } // namespace random
 353
 354 using random::poisson_distribution;
 355
 356 } // namespace boost
 357
 358 #include <boost/random/detail/enable_warnings.hpp>
 359
 360 #endif // BOOST_RANDOM_POISSON_DISTRIBUTION_HPP