Chris@16: /* boost random/binomial_distribution.hpp header file
Chris@16:  *
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: #ifndef BOOST_RANDOM_BINOMIAL_DISTRIBUTION_HPP_INCLUDED
Chris@16: #define BOOST_RANDOM_BINOMIAL_DISTRIBUTION_HPP_INCLUDED
Chris@16: 
Chris@16: #include <boost/config/no_tr1/cmath.hpp>
Chris@16: #include <cstdlib>
Chris@16: #include <iosfwd>
Chris@16: 
Chris@16: #include <boost/random/detail/config.hpp>
Chris@16: #include <boost/random/uniform_01.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 binomial_table {
Chris@16:     static const RealType table[10];
Chris@16: };
Chris@16: 
Chris@16: template<class RealType>
Chris@16: const RealType binomial_table<RealType>::table[10] = {
Chris@16:     0.08106146679532726,
Chris@16:     0.04134069595540929,
Chris@16:     0.02767792568499834,
Chris@16:     0.02079067210376509,
Chris@16:     0.01664469118982119,
Chris@16:     0.01387612882307075,
Chris@16:     0.01189670994589177,
Chris@16:     0.01041126526197209,
Chris@16:     0.009255462182712733,
Chris@16:     0.008330563433362871
Chris@16: };
Chris@16: 
Chris@16: }
Chris@16: 
Chris@16: /**
Chris@16:  * The binomial distribution is an integer valued distribution with
Chris@16:  * two parameters, @c t and @c p.  The values of the distribution
Chris@16:  * are within the range [0,t].
Chris@16:  *
Chris@16:  * The distribution function is
Chris@16:  * \f$\displaystyle P(k) = {t \choose k}p^k(1-p)^{t-k}\f$.
Chris@16:  *
Chris@16:  * The algorithm used is the BTRD algorithm described in
Chris@16:  *
Chris@16:  *  @blockquote
Chris@16:  *  "The generation of binomial random variates", Wolfgang Hormann,
Chris@16:  *  Journal of Statistical Computation and Simulation, Volume 46,
Chris@16:  *  Issue 1 & 2 April 1993 , pages 101 - 110
Chris@16:  *  @endblockquote
Chris@16:  */
Chris@16: template<class IntType = int, class RealType = double>
Chris@16: class binomial_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 binomial_distribution distribution_type;
Chris@16:         /**
Chris@16:          * Construct a param_type object.  @c t and @c p
Chris@16:          * are the parameters of the distribution.
Chris@16:          *
Chris@16:          * Requires: t >=0 && 0 <= p <= 1
Chris@16:          */
Chris@16:         explicit param_type(IntType t_arg = 1, RealType p_arg = RealType (0.5))
Chris@16:           : _t(t_arg), _p(p_arg)
Chris@16:         {}
Chris@16:         /** Returns the @c t parameter of the distribution. */
Chris@16:         IntType t() const { return _t; }
Chris@16:         /** Returns the @c p parameter of the distribution. */
Chris@16:         RealType p() const { return _p; }
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._p << " " << parm._t;
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._p >> std::ws >> parm._t;
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._t == rhs._t && lhs._p == rhs._p;
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:         IntType _t;
Chris@16:         RealType _p;
Chris@16:     };
Chris@16:     
Chris@16:     /**
Chris@16:      * Construct a @c binomial_distribution object. @c t and @c p
Chris@16:      * are the parameters of the distribution.
Chris@16:      *
Chris@16:      * Requires: t >=0 && 0 <= p <= 1
Chris@16:      */
Chris@16:     explicit binomial_distribution(IntType t_arg = 1,
Chris@16:                                    RealType p_arg = RealType(0.5))
Chris@16:       : _t(t_arg), _p(p_arg)
Chris@16:     {
Chris@16:         init();
Chris@16:     }
Chris@16:     
Chris@16:     /**
Chris@16:      * Construct an @c binomial_distribution object from the
Chris@16:      * parameters.
Chris@16:      */
Chris@16:     explicit binomial_distribution(const param_type& parm)
Chris@16:       : _t(parm.t()), _p(parm.p())
Chris@16:     {
Chris@16:         init();
Chris@16:     }
Chris@16:     
Chris@16:     /**
Chris@16:      * Returns a random variate distributed according to the
Chris@16:      * binomial distribution.
Chris@16:      */
Chris@16:     template<class URNG>
Chris@16:     IntType operator()(URNG& urng) const
Chris@16:     {
Chris@16:         if(use_inversion()) {
Chris@16:             if(0.5 < _p) {
Chris@16:                 return _t - invert(_t, 1-_p, urng);
Chris@16:             } else {
Chris@16:                 return invert(_t, _p, urng);
Chris@16:             }
Chris@16:         } else if(0.5 < _p) {
Chris@16:             return _t - generate(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:      * binomial distribution with parameters specified by @c param.
Chris@16:      */
Chris@16:     template<class URNG>
Chris@16:     IntType operator()(URNG& urng, const param_type& parm) const
Chris@16:     {
Chris@16:         return binomial_distribution(parm)(urng);
Chris@16:     }
Chris@16: 
Chris@16:     /** Returns the @c t parameter of the distribution. */
Chris@16:     IntType t() const { return _t; }
Chris@16:     /** Returns the @c p parameter of the distribution. */
Chris@16:     RealType p() const { return _p; }
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 { return _t; }
Chris@16: 
Chris@16:     /** Returns the parameters of the distribution. */
Chris@16:     param_type param() const { return param_type(_t, _p); }
Chris@16:     /** Sets parameters of the distribution. */
Chris@16:     void param(const param_type& parm)
Chris@16:     {
Chris@16:         _t = parm.t();
Chris@16:         _p = parm.p();
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 binomial_distribution& bd)
Chris@16:     {
Chris@16:         os << bd.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, binomial_distribution& bd)
Chris@16:     {
Chris@16:         bd.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 binomial_distribution& lhs,
Chris@16:                            const binomial_distribution& rhs)
Chris@16:     {
Chris@16:         return lhs._t == rhs._t && lhs._p == rhs._p;
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 binomial_distribution& lhs,
Chris@16:                            const binomial_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:         // BTRD is safe when np >= 10
Chris@16:         return m < 11;
Chris@16:     }
Chris@16: 
Chris@16:     // computes the correction factor for the Stirling approximation
Chris@16:     // for log(k!)
Chris@16:     static RealType fc(IntType k)
Chris@16:     {
Chris@16:         if(k < 10) return detail::binomial_table<RealType>::table[k];
Chris@16:         else {
Chris@16:             RealType ikp1 = RealType(1) / (k + 1);
Chris@16:             return (RealType(1)/12
Chris@16:                  - (RealType(1)/360
Chris@16:                  - (RealType(1)/1260)*(ikp1*ikp1))*(ikp1*ikp1))*ikp1;
Chris@16:         }
Chris@16:     }
Chris@16: 
Chris@16:     void init()
Chris@16:     {
Chris@16:         using std::sqrt;
Chris@16:         using std::pow;
Chris@16: 
Chris@16:         RealType p = (0.5 < _p)? (1 - _p) : _p;
Chris@16:         IntType t = _t;
Chris@16:         
Chris@16:         m = static_cast<IntType>((t+1)*p);
Chris@16: 
Chris@16:         if(use_inversion()) {
Chris@16:             q_n = pow((1 - p), static_cast<RealType>(t));
Chris@16:         } else {
Chris@16:             btrd.r = p/(1-p);
Chris@16:             btrd.nr = (t+1)*btrd.r;
Chris@16:             btrd.npq = t*p*(1-p);
Chris@16:             RealType sqrt_npq = sqrt(btrd.npq);
Chris@16:             btrd.b = 1.15 + 2.53 * sqrt_npq;
Chris@16:             btrd.a = -0.0873 + 0.0248*btrd.b + 0.01*p;
Chris@16:             btrd.c = t*p + 0.5;
Chris@16:             btrd.alpha = (2.83 + 5.1/btrd.b) * sqrt_npq;
Chris@16:             btrd.v_r = 0.92 - 4.2/btrd.b;
Chris@16:             btrd.u_rv_r = 0.86*btrd.v_r;
Chris@16:         }
Chris@16:     }
Chris@16: 
Chris@16:     template<class URNG>
Chris@16:     result_type 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 <= btrd.u_rv_r) {
Chris@16:                 RealType u = v/btrd.v_r - 0.43;
Chris@16:                 return static_cast<IntType>(floor(
Chris@16:                     (2*btrd.a/(0.5 - abs(u)) + btrd.b)*u + btrd.c));
Chris@16:             }
Chris@16: 
Chris@16:             if(v >= btrd.v_r) {
Chris@16:                 u = uniform_01<RealType>()(urng) - 0.5;
Chris@16:             } else {
Chris@16:                 u = v/btrd.v_r - 0.93;
Chris@16:                 u = ((u < 0)? -0.5 : 0.5) - u;
Chris@16:                 v = uniform_01<RealType>()(urng) * btrd.v_r;
Chris@16:             }
Chris@16: 
Chris@16:             RealType us = 0.5 - abs(u);
Chris@16:             IntType k = static_cast<IntType>(floor((2*btrd.a/us + btrd.b)*u + btrd.c));
Chris@16:             if(k < 0 || k > _t) continue;
Chris@16:             v = v*btrd.alpha/(btrd.a/(us*us) + btrd.b);
Chris@16:             RealType km = abs(k - m);
Chris@16:             if(km <= 15) {
Chris@16:                 RealType f = 1;
Chris@16:                 if(m < k) {
Chris@16:                     IntType i = m;
Chris@16:                     do {
Chris@16:                         ++i;
Chris@16:                         f = f*(btrd.nr/i - btrd.r);
Chris@16:                     } while(i != k);
Chris@16:                 } else if(m > k) {
Chris@16:                     IntType i = k;
Chris@16:                     do {
Chris@16:                         ++i;
Chris@16:                         v = v*(btrd.nr/i - btrd.r);
Chris@16:                     } while(i != m);
Chris@16:                 }
Chris@16:                 if(v <= f) return k;
Chris@16:                 else continue;
Chris@16:             } else {
Chris@16:                 // final acceptance/rejection
Chris@16:                 v = log(v);
Chris@16:                 RealType rho =
Chris@16:                     (km/btrd.npq)*(((km/3. + 0.625)*km + 1./6)/btrd.npq + 0.5);
Chris@16:                 RealType t = -km*km/(2*btrd.npq);
Chris@16:                 if(v < t - rho) return k;
Chris@16:                 if(v > t + rho) continue;
Chris@16: 
Chris@16:                 IntType nm = _t - m + 1;
Chris@16:                 RealType h = (m + 0.5)*log((m + 1)/(btrd.r*nm))
Chris@16:                            + fc(m) + fc(_t - m);
Chris@16: 
Chris@16:                 IntType nk = _t - k + 1;
Chris@16:                 if(v <= h + (_t+1)*log(static_cast<RealType>(nm)/nk)
Chris@16:                           + (k + 0.5)*log(nk*btrd.r/(k+1))
Chris@16:                           - fc(k)
Chris@16:                           - fc(_t - k))
Chris@16:                 {
Chris@16:                     return k;
Chris@16:                 } else {
Chris@16:                     continue;
Chris@16:                 }
Chris@16:             }
Chris@16:         }
Chris@16:     }
Chris@16: 
Chris@16:     template<class URNG>
Chris@16:     IntType invert(IntType t, RealType p, URNG& urng) const
Chris@16:     {
Chris@16:         RealType q = 1 - p;
Chris@16:         RealType s = p / q;
Chris@16:         RealType a = (t + 1) * s;
Chris@16:         RealType r = q_n;
Chris@16:         RealType u = uniform_01<RealType>()(urng);
Chris@16:         IntType x = 0;
Chris@16:         while(u > r) {
Chris@16:             u = u - r;
Chris@16:             ++x;
Chris@101:             RealType r1 = ((a/x) - s) * r;
Chris@101:             // If r gets too small then the round-off error
Chris@101:             // becomes a problem.  At this point, p(i) is
Chris@101:             // decreasing exponentially, so if we just call
Chris@101:             // it 0, it's close enough.  Note that the
Chris@101:             // minimum value of q_n is about 1e-7, so we
Chris@101:             // may need to be a little careful to make sure that
Chris@101:             // we don't terminate the first time through the loop
Chris@101:             // for float.  (Hence the test that r is decreasing)
Chris@101:             if(r1 < std::numeric_limits<RealType>::epsilon() && r1 < r) {
Chris@101:                 break;
Chris@101:             }
Chris@101:             r = r1;
Chris@16:         }
Chris@16:         return x;
Chris@16:     }
Chris@16: 
Chris@16:     // parameters
Chris@16:     IntType _t;
Chris@16:     RealType _p;
Chris@16: 
Chris@16:     // common data
Chris@16:     IntType m;
Chris@16: 
Chris@16:     union {
Chris@16:         // for btrd
Chris@16:         struct {
Chris@16:             RealType r;
Chris@16:             RealType nr;
Chris@16:             RealType npq;
Chris@16:             RealType b;
Chris@16:             RealType a;
Chris@16:             RealType c;
Chris@16:             RealType alpha;
Chris@16:             RealType v_r;
Chris@16:             RealType u_rv_r;
Chris@16:         } btrd;
Chris@16:         // for inversion
Chris@16:         RealType q_n;
Chris@16:     };
Chris@16: 
Chris@16:     /// @endcond
Chris@16: };
Chris@16: 
Chris@16: }
Chris@16: 
Chris@16: // backwards compatibility
Chris@16: using random::binomial_distribution;
Chris@16: 
Chris@16: }
Chris@16: 
Chris@16: #include <boost/random/detail/enable_warnings.hpp>
Chris@16: 
Chris@16: #endif