1 /* boost random/binomial_distribution.hpp header file
3 * Copyright Steven Watanabe 2010
4 * Distributed under the Boost Software License, Version 1.0. (See
5 * accompanying file LICENSE_1_0.txt or copy at
6 * http://www.boost.org/LICENSE_1_0.txt)
8 * See http://www.boost.org for most recent version including documentation.
13 #ifndef BOOST_RANDOM_BINOMIAL_DISTRIBUTION_HPP_INCLUDED
14 #define BOOST_RANDOM_BINOMIAL_DISTRIBUTION_HPP_INCLUDED
16 #include <boost/config/no_tr1/cmath.hpp>
20 #include <boost/random/detail/config.hpp>
21 #include <boost/random/uniform_01.hpp>
23 #include <boost/random/detail/disable_warnings.hpp>
30 template<class RealType>
31 struct binomial_table {
32 static const RealType table[10];
35 template<class RealType>
36 const RealType binomial_table<RealType>::table[10] = {
52 * The binomial distribution is an integer valued distribution with
53 * two parameters, @c t and @c p. The values of the distribution
54 * are within the range [0,t].
56 * The distribution function is
57 * \f$\displaystyle P(k) = {t \choose k}p^k(1-p)^{t-k}\f$.
59 * The algorithm used is the BTRD algorithm described in
62 * "The generation of binomial random variates", Wolfgang Hormann,
63 * Journal of Statistical Computation and Simulation, Volume 46,
64 * Issue 1 & 2 April 1993 , pages 101 - 110
67 template<class IntType = int, class RealType = double>
68 class binomial_distribution {
70 typedef IntType result_type;
71 typedef RealType input_type;
75 typedef binomial_distribution distribution_type;
77 * Construct a param_type object. @c t and @c p
78 * are the parameters of the distribution.
80 * Requires: t >=0 && 0 <= p <= 1
82 explicit param_type(IntType t_arg = 1, RealType p_arg = RealType (0.5))
83 : _t(t_arg), _p(p_arg)
85 /** Returns the @c t parameter of the distribution. */
86 IntType t() const { return _t; }
87 /** Returns the @c p parameter of the distribution. */
88 RealType p() const { return _p; }
89 #ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
90 /** Writes the parameters of the distribution to a @c std::ostream. */
91 template<class CharT, class Traits>
92 friend std::basic_ostream<CharT,Traits>&
93 operator<<(std::basic_ostream<CharT,Traits>& os,
94 const param_type& parm)
96 os << parm._p << " " << parm._t;
100 /** Reads the parameters of the distribution from a @c std::istream. */
101 template<class CharT, class Traits>
102 friend std::basic_istream<CharT,Traits>&
103 operator>>(std::basic_istream<CharT,Traits>& is, param_type& parm)
105 is >> parm._p >> std::ws >> parm._t;
109 /** Returns true if the parameters have the same values. */
110 friend bool operator==(const param_type& lhs, const param_type& rhs)
112 return lhs._t == rhs._t && lhs._p == rhs._p;
114 /** Returns true if the parameters have different values. */
115 friend bool operator!=(const param_type& lhs, const param_type& rhs)
117 return !(lhs == rhs);
125 * Construct a @c binomial_distribution object. @c t and @c p
126 * are the parameters of the distribution.
128 * Requires: t >=0 && 0 <= p <= 1
130 explicit binomial_distribution(IntType t_arg = 1,
131 RealType p_arg = RealType(0.5))
132 : _t(t_arg), _p(p_arg)
138 * Construct an @c binomial_distribution object from the
141 explicit binomial_distribution(const param_type& parm)
142 : _t(parm.t()), _p(parm.p())
148 * Returns a random variate distributed according to the
149 * binomial distribution.
152 IntType operator()(URNG& urng) const
154 if(use_inversion()) {
156 return _t - invert(_t, 1-_p, urng);
158 return invert(_t, _p, urng);
160 } else if(0.5 < _p) {
161 return _t - generate(urng);
163 return generate(urng);
168 * Returns a random variate distributed according to the
169 * binomial distribution with parameters specified by @c param.
172 IntType operator()(URNG& urng, const param_type& parm) const
174 return binomial_distribution(parm)(urng);
177 /** Returns the @c t parameter of the distribution. */
178 IntType t() const { return _t; }
179 /** Returns the @c p parameter of the distribution. */
180 RealType p() const { return _p; }
182 /** Returns the smallest value that the distribution can produce. */
183 IntType min BOOST_PREVENT_MACRO_SUBSTITUTION() const { return 0; }
184 /** Returns the largest value that the distribution can produce. */
185 IntType max BOOST_PREVENT_MACRO_SUBSTITUTION() const { return _t; }
187 /** Returns the parameters of the distribution. */
188 param_type param() const { return param_type(_t, _p); }
189 /** Sets parameters of the distribution. */
190 void param(const param_type& parm)
198 * Effects: Subsequent uses of the distribution do not depend
199 * on values produced by any engine prior to invoking reset.
203 #ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
204 /** Writes the parameters of the distribution to a @c std::ostream. */
205 template<class CharT, class Traits>
206 friend std::basic_ostream<CharT,Traits>&
207 operator<<(std::basic_ostream<CharT,Traits>& os,
208 const binomial_distribution& bd)
214 /** Reads the parameters of the distribution from a @c std::istream. */
215 template<class CharT, class Traits>
216 friend std::basic_istream<CharT,Traits>&
217 operator>>(std::basic_istream<CharT,Traits>& is, binomial_distribution& bd)
224 /** Returns true if the two distributions will produce the same
225 sequence of values, given equal generators. */
226 friend bool operator==(const binomial_distribution& lhs,
227 const binomial_distribution& rhs)
229 return lhs._t == rhs._t && lhs._p == rhs._p;
231 /** Returns true if the two distributions could produce different
232 sequences of values, given equal generators. */
233 friend bool operator!=(const binomial_distribution& lhs,
234 const binomial_distribution& rhs)
236 return !(lhs == rhs);
241 /// @cond show_private
243 template<class CharT, class Traits>
244 void read(std::basic_istream<CharT, Traits>& is) {
251 bool use_inversion() const
253 // BTRD is safe when np >= 10
257 // computes the correction factor for the Stirling approximation
259 static RealType fc(IntType k)
261 if(k < 10) return detail::binomial_table<RealType>::table[k];
263 RealType ikp1 = RealType(1) / (k + 1);
264 return (RealType(1)/12
266 - (RealType(1)/1260)*(ikp1*ikp1))*(ikp1*ikp1))*ikp1;
275 RealType p = (0.5 < _p)? (1 - _p) : _p;
278 m = static_cast<IntType>((t+1)*p);
280 if(use_inversion()) {
281 _u.q_n = pow((1 - p), static_cast<RealType>(t));
284 _u.btrd.nr = (t+1)*_u.btrd.r;
285 _u.btrd.npq = t*p*(1-p);
286 RealType sqrt_npq = sqrt(_u.btrd.npq);
287 _u.btrd.b = 1.15 + 2.53 * sqrt_npq;
288 _u.btrd.a = -0.0873 + 0.0248*_u.btrd.b + 0.01*p;
289 _u.btrd.c = t*p + 0.5;
290 _u.btrd.alpha = (2.83 + 5.1/_u.btrd.b) * sqrt_npq;
291 _u.btrd.v_r = 0.92 - 4.2/_u.btrd.b;
292 _u.btrd.u_rv_r = 0.86*_u.btrd.v_r;
297 result_type generate(URNG& urng) const
305 RealType v = uniform_01<RealType>()(urng);
306 if(v <= _u.btrd.u_rv_r) {
307 u = v/_u.btrd.v_r - 0.43;
308 return static_cast<IntType>(floor(
309 (2*_u.btrd.a/(0.5 - abs(u)) + _u.btrd.b)*u + _u.btrd.c));
312 if(v >= _u.btrd.v_r) {
313 u = uniform_01<RealType>()(urng) - 0.5;
315 u = v/_u.btrd.v_r - 0.93;
316 u = ((u < 0)? -0.5 : 0.5) - u;
317 v = uniform_01<RealType>()(urng) * _u.btrd.v_r;
320 RealType us = 0.5 - abs(u);
321 IntType k = static_cast<IntType>(floor((2*_u.btrd.a/us + _u.btrd.b)*u + _u.btrd.c));
322 if(k < 0 || k > _t) continue;
323 v = v*_u.btrd.alpha/(_u.btrd.a/(us*us) + _u.btrd.b);
324 RealType km = abs(k - m);
331 f = f*(_u.btrd.nr/i - _u.btrd.r);
337 v = v*(_u.btrd.nr/i - _u.btrd.r);
343 // final acceptance/rejection
346 (km/_u.btrd.npq)*(((km/3. + 0.625)*km + 1./6)/_u.btrd.npq + 0.5);
347 RealType t = -km*km/(2*_u.btrd.npq);
348 if(v < t - rho) return k;
349 if(v > t + rho) continue;
351 IntType nm = _t - m + 1;
352 RealType h = (m + 0.5)*log((m + 1)/(_u.btrd.r*nm))
353 + fc(m) + fc(_t - m);
355 IntType nk = _t - k + 1;
356 if(v <= h + (_t+1)*log(static_cast<RealType>(nm)/nk)
357 + (k + 0.5)*log(nk*_u.btrd.r/(k+1))
370 IntType invert(IntType t, RealType p, URNG& urng) const
374 RealType a = (t + 1) * s;
376 RealType u = uniform_01<RealType>()(urng);
381 RealType r1 = ((a/x) - s) * r;
382 // If r gets too small then the round-off error
383 // becomes a problem. At this point, p(i) is
384 // decreasing exponentially, so if we just call
385 // it 0, it's close enough. Note that the
386 // minimum value of q_n is about 1e-7, so we
387 // may need to be a little careful to make sure that
388 // we don't terminate the first time through the loop
389 // for float. (Hence the test that r is decreasing)
390 if(r1 < std::numeric_limits<RealType>::epsilon() && r1 < r) {
427 // backwards compatibility
428 using random::binomial_distribution;
432 #include <boost/random/detail/enable_warnings.hpp>