boost/random/binomial_distribution.hpp

0001 /* boost random/binomial_distribution.hpp header file
0002  *
0003  * Copyright Steven Watanabe 2010
0004  * Distributed under the Boost Software License, Version 1.0. (See
0005  * accompanying file LICENSE_1_0.txt or copy at
0006  * http://www.boost.org/LICENSE_1_0.txt)
0007  *
0008  * See http://www.boost.org for most recent version including documentation.
0009  *
0010  * $Id$
0011  */
0012
0013 #ifndef BOOST_RANDOM_BINOMIAL_DISTRIBUTION_HPP_INCLUDED
0014 #define BOOST_RANDOM_BINOMIAL_DISTRIBUTION_HPP_INCLUDED
0015
0016 #include <boost/config/no_tr1/cmath.hpp>
0017 #include <cstdlib>
0018 #include <iosfwd>
0019
0020 #include <boost/random/detail/config.hpp>
0021 #include <boost/random/uniform_01.hpp>
0022
0023 #include <boost/random/detail/disable_warnings.hpp>
0024
0025 namespace boost {
0026 namespace random {
0027
0028 namespace detail {
0029
0030 template<class RealType>
0031 struct binomial_table {
0032     static const RealType table[10];
0033 };
0034
0035 template<class RealType>
0036 const RealType binomial_table<RealType>::table[10] = {
0037     0.08106146679532726,
0038     0.04134069595540929,
0039     0.02767792568499834,
0040     0.02079067210376509,
0041     0.01664469118982119,
0042     0.01387612882307075,
0043     0.01189670994589177,
0044     0.01041126526197209,
0045     0.009255462182712733,
0046     0.008330563433362871
0047 };
0048
0049 }
0050
0051 /**
0052  * The binomial distribution is an integer valued distribution with
0053  * two parameters, @c t and @c p.  The values of the distribution
0054  * are within the range [0,t].
0055  *
0056  * The distribution function is
0057  * \f$\displaystyle P(k) = {t \choose k}p^k(1-p)^{t-k}\f$.
0058  *
0059  * The algorithm used is the BTRD algorithm described in
0060  *
0061  *  @blockquote
0062  *  "The generation of binomial random variates", Wolfgang Hormann,
0063  *  Journal of Statistical Computation and Simulation, Volume 46,
0064  *  Issue 1 & 2 April 1993 , pages 101 - 110
0065  *  @endblockquote
0066  */
0067 template<class IntType = int, class RealType = double>
0068 class binomial_distribution {
0069 public:
0070     typedef IntType result_type;
0071     typedef RealType input_type;
0072
0073     class param_type {
0074     public:
0075         typedef binomial_distribution distribution_type;
0076         /**
0077          * Construct a param_type object.  @c t and @c p
0078          * are the parameters of the distribution.
0079          *
0080          * Requires: t >=0 && 0 <= p <= 1
0081          */
0082         explicit param_type(IntType t_arg = 1, RealType p_arg = RealType (0.5))
0083           : _t(t_arg), _p(p_arg)
0084         {}
0085         /** Returns the @c t parameter of the distribution. */
0086         IntType t() const { return _t; }
0087         /** Returns the @c p parameter of the distribution. */
0088         RealType p() const { return _p; }
0089 #ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
0090         /** Writes the parameters of the distribution to a @c std::ostream. */
0091         template<class CharT, class Traits>
0092         friend std::basic_ostream<CharT,Traits>&
0093         operator<<(std::basic_ostream<CharT,Traits>& os,
0094                    const param_type& parm)
0095         {
0096             os << parm._p << " " << parm._t;
0097             return os;
0098         }
0099
0100         /** Reads the parameters of the distribution from a @c std::istream. */
0101         template<class CharT, class Traits>
0102         friend std::basic_istream<CharT,Traits>&
0103         operator>>(std::basic_istream<CharT,Traits>& is, param_type& parm)
0104         {
0105             is >> parm._p >> std::ws >> parm._t;
0106             return is;
0107         }
0108 #endif
0109         /** Returns true if the parameters have the same values. */
0110         friend bool operator==(const param_type& lhs, const param_type& rhs)
0111         {
0112             return lhs._t == rhs._t && lhs._p == rhs._p;
0113         }
0114         /** Returns true if the parameters have different values. */
0115         friend bool operator!=(const param_type& lhs, const param_type& rhs)
0116         {
0117             return !(lhs == rhs);
0118         }
0119     private:
0120         IntType _t;
0121         RealType _p;
0122     };
0123
0124     /**
0125      * Construct a @c binomial_distribution object. @c t and @c p
0126      * are the parameters of the distribution.
0127      *
0128      * Requires: t >=0 && 0 <= p <= 1
0129      */
0130     explicit binomial_distribution(IntType t_arg = 1,
0131                                    RealType p_arg = RealType(0.5))
0132       : _t(t_arg), _p(p_arg)
0133     {
0134         init();
0135     }
0136
0137     /**
0138      * Construct an @c binomial_distribution object from the
0139      * parameters.
0140      */
0141     explicit binomial_distribution(const param_type& parm)
0142       : _t(parm.t()), _p(parm.p())
0143     {
0144         init();
0145     }
0146
0147     /**
0148      * Returns a random variate distributed according to the
0149      * binomial distribution.
0150      */
0151     template<class URNG>
0152     IntType operator()(URNG& urng) const
0153     {
0154         if(use_inversion()) {
0155             if(0.5 < _p) {
0156                 return _t - invert(_t, 1-_p, urng);
0157             } else {
0158                 return invert(_t, _p, urng);
0159             }
0160         } else if(0.5 < _p) {
0161             return _t - generate(urng);
0162         } else {
0163             return generate(urng);
0164         }
0165     }
0166
0167     /**
0168      * Returns a random variate distributed according to the
0169      * binomial distribution with parameters specified by @c param.
0170      */
0171     template<class URNG>
0172     IntType operator()(URNG& urng, const param_type& parm) const
0173     {
0174         return binomial_distribution(parm)(urng);
0175     }
0176
0177     /** Returns the @c t parameter of the distribution. */
0178     IntType t() const { return _t; }
0179     /** Returns the @c p parameter of the distribution. */
0180     RealType p() const { return _p; }
0181
0182     /** Returns the smallest value that the distribution can produce. */
0183     IntType min BOOST_PREVENT_MACRO_SUBSTITUTION() const { return 0; }
0184     /** Returns the largest value that the distribution can produce. */
0185     IntType max BOOST_PREVENT_MACRO_SUBSTITUTION() const { return _t; }
0186
0187     /** Returns the parameters of the distribution. */
0188     param_type param() const { return param_type(_t, _p); }
0189     /** Sets parameters of the distribution. */
0190     void param(const param_type& parm)
0191     {
0192         _t = parm.t();
0193         _p = parm.p();
0194         init();
0195     }
0196
0197     /**
0198      * Effects: Subsequent uses of the distribution do not depend
0199      * on values produced by any engine prior to invoking reset.
0200      */
0201     void reset() { }
0202
0203 #ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
0204     /** Writes the parameters of the distribution to a @c std::ostream. */
0205     template<class CharT, class Traits>
0206     friend std::basic_ostream<CharT,Traits>&
0207     operator<<(std::basic_ostream<CharT,Traits>& os,
0208                const binomial_distribution& bd)
0209     {
0210         os << bd.param();
0211         return os;
0212     }
0213
0214     /** Reads the parameters of the distribution from a @c std::istream. */
0215     template<class CharT, class Traits>
0216     friend std::basic_istream<CharT,Traits>&
0217     operator>>(std::basic_istream<CharT,Traits>& is, binomial_distribution& bd)
0218     {
0219         bd.read(is);
0220         return is;
0221     }
0222 #endif
0223
0224     /** Returns true if the two distributions will produce the same
0225         sequence of values, given equal generators. */
0226     friend bool operator==(const binomial_distribution& lhs,
0227                            const binomial_distribution& rhs)
0228     {
0229         return lhs._t == rhs._t && lhs._p == rhs._p;
0230     }
0231     /** Returns true if the two distributions could produce different
0232         sequences of values, given equal generators. */
0233     friend bool operator!=(const binomial_distribution& lhs,
0234                            const binomial_distribution& rhs)
0235     {
0236         return !(lhs == rhs);
0237     }
0238
0239 private:
0240
0241     /// @cond show_private
0242
0243     template<class CharT, class Traits>
0244     void read(std::basic_istream<CharT, Traits>& is) {
0245         param_type parm;
0246         if(is >> parm) {
0247             param(parm);
0248         }
0249     }
0250
0251     bool use_inversion() const
0252     {
0253         // BTRD is safe when np >= 10
0254         return m < 11;
0255     }
0256
0257     // computes the correction factor for the Stirling approximation
0258     // for log(k!)
0259     static RealType fc(IntType k)
0260     {
0261         if(k < 10) return detail::binomial_table<RealType>::table[k];
0262         else {
0263             RealType ikp1 = RealType(1) / (k + 1);
0264             return (RealType(1)/12
0265                  - (RealType(1)/360
0266                  - (RealType(1)/1260)*(ikp1*ikp1))*(ikp1*ikp1))*ikp1;
0267         }
0268     }
0269
0270     void init()
0271     {
0272         using std::sqrt;
0273         using std::pow;
0274
0275         RealType p = (0.5 < _p)? (1 - _p) : _p;
0276         IntType t = _t;
0277
0278         m = static_cast<IntType>((t+1)*p);
0279
0280         if(use_inversion()) {
0281             _u.q_n = pow((1 - p), static_cast<RealType>(t));
0282         } else {
0283             _u.btrd.r = p/(1-p);
0284             _u.btrd.nr = (t+1)*_u.btrd.r;
0285             _u.btrd.npq = t*p*(1-p);
0286             RealType sqrt_npq = sqrt(_u.btrd.npq);
0287             _u.btrd.b = 1.15 + 2.53 * sqrt_npq;
0288             _u.btrd.a = -0.0873 + 0.0248*_u.btrd.b + 0.01*p;
0289             _u.btrd.c = t*p + 0.5;
0290             _u.btrd.alpha = (2.83 + 5.1/_u.btrd.b) * sqrt_npq;
0291             _u.btrd.v_r = 0.92 - 4.2/_u.btrd.b;
0292             _u.btrd.u_rv_r = 0.86*_u.btrd.v_r;
0293         }
0294     }
0295
0296     template<class URNG>
0297     result_type generate(URNG& urng) const
0298     {
0299         using std::floor;
0300         using std::abs;
0301         using std::log;
0302
0303         while(true) {
0304             RealType u;
0305             RealType v = uniform_01<RealType>()(urng);
0306             if(v <= _u.btrd.u_rv_r) {
0307                 u = v/_u.btrd.v_r - 0.43;
0308                 return static_cast<IntType>(floor(
0309                     (2*_u.btrd.a/(0.5 - abs(u)) + _u.btrd.b)*u + _u.btrd.c));
0310             }
0311
0312             if(v >= _u.btrd.v_r) {
0313                 u = uniform_01<RealType>()(urng) - 0.5;
0314             } else {
0315                 u = v/_u.btrd.v_r - 0.93;
0316                 u = ((u < 0)? -0.5 : 0.5) - u;
0317                 v = uniform_01<RealType>()(urng) * _u.btrd.v_r;
0318             }
0319
0320             RealType us = 0.5 - abs(u);
0321             IntType k = static_cast<IntType>(floor((2*_u.btrd.a/us + _u.btrd.b)*u + _u.btrd.c));
0322             if(k < 0 || k > _t) continue;
0323             v = v*_u.btrd.alpha/(_u.btrd.a/(us*us) + _u.btrd.b);
0324             RealType km = abs(k - m);
0325             if(km <= 15) {
0326                 RealType f = 1;
0327                 if(m < k) {
0328                     IntType i = m;
0329                     do {
0330                         ++i;
0331                         f = f*(_u.btrd.nr/i - _u.btrd.r);
0332                     } while(i != k);
0333                 } else if(m > k) {
0334                     IntType i = k;
0335                     do {
0336                         ++i;
0337                         v = v*(_u.btrd.nr/i - _u.btrd.r);
0338                     } while(i != m);
0339                 }
0340                 if(v <= f) return k;
0341                 else continue;
0342             } else {
0343                 // final acceptance/rejection
0344                 v = log(v);
0345                 RealType rho =
0346                     (km/_u.btrd.npq)*(((km/3. + 0.625)*km + 1./6)/_u.btrd.npq + 0.5);
0347                 RealType t = -km*km/(2*_u.btrd.npq);
0348                 if(v < t - rho) return k;
0349                 if(v > t + rho) continue;
0350
0351                 IntType nm = _t - m + 1;
0352                 RealType h = (m + 0.5)*log((m + 1)/(_u.btrd.r*nm))
0353                            + fc(m) + fc(_t - m);
0354
0355                 IntType nk = _t - k + 1;
0356                 if(v <= h + (_t+1)*log(static_cast<RealType>(nm)/nk)
0357                           + (k + 0.5)*log(nk*_u.btrd.r/(k+1))
0358                           - fc(k)
0359                           - fc(_t - k))
0360                 {
0361                     return k;
0362                 } else {
0363                     continue;
0364                 }
0365             }
0366         }
0367     }
0368
0369     template<class URNG>
0370     IntType invert(IntType t, RealType p, URNG& urng) const
0371     {
0372         RealType q = 1 - p;
0373         RealType s = p / q;
0374         RealType a = (t + 1) * s;
0375         RealType r = _u.q_n;
0376         RealType u = uniform_01<RealType>()(urng);
0377         IntType x = 0;
0378         while(u > r) {
0379             u = u - r;
0380             ++x;
0381             RealType r1 = ((a/x) - s) * r;
0382             // If r gets too small then the round-off error
0383             // becomes a problem.  At this point, p(i) is
0384             // decreasing exponentially, so if we just call
0385             // it 0, it's close enough.  Note that the
0386             // minimum value of q_n is about 1e-7, so we
0387             // may need to be a little careful to make sure that
0388             // we don't terminate the first time through the loop
0389             // for float.  (Hence the test that r is decreasing)
0390             if(r1 < std::numeric_limits<RealType>::epsilon() && r1 < r) {
0391                 break;
0392             }
0393             r = r1;
0394         }
0395         return x;
0396     }
0397
0398     // parameters
0399     IntType _t;
0400     RealType _p;
0401
0402     // common data
0403     IntType m;
0404
0405     union {
0406         // for btrd
0407         struct {
0408             RealType r;
0409             RealType nr;
0410             RealType npq;
0411             RealType b;
0412             RealType a;
0413             RealType c;
0414             RealType alpha;
0415             RealType v_r;
0416             RealType u_rv_r;
0417         } btrd;
0418         // for inversion
0419         RealType q_n;
0420     } _u;
0421
0422     /// @endcond
0423 };
0424
0425 }
0426
0427 // backwards compatibility
0428 using random::binomial_distribution;
0429
0430 }
0431
0432 #include <boost/random/detail/enable_warnings.hpp>
0433
0434 #endif