1/* boost random/gamma_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$
12 *
13 */
14
15#ifndef BOOST_RANDOM_GAMMA_DISTRIBUTION_HPP
16#define BOOST_RANDOM_GAMMA_DISTRIBUTION_HPP
17
18#include <boost/config/no_tr1/cmath.hpp>
19#include <istream>
20#include <iosfwd>
21#include <boost/assert.hpp>
22#include <boost/limits.hpp>
23#include <boost/static_assert.hpp>
24#include <boost/random/detail/config.hpp>
25#include <boost/random/exponential_distribution.hpp>
26
27namespace boost {
28namespace random {
29
30// The algorithm is taken from Knuth
31
32/**
33 * The gamma distribution is a continuous distribution with two
34 * parameters alpha and beta. It produces values > 0.
35 *
36 * It has
37 * \f$\displaystyle p(x) = x^{\alpha-1}\frac{e^{-x/\beta}}{\beta^\alpha\Gamma(\alpha)}\f$.
38 */
39template<class RealType = double>
40class gamma_distribution
41{
42public:
43 typedef RealType input_type;
44 typedef RealType result_type;
45
46 class param_type
47 {
48 public:
49 typedef gamma_distribution distribution_type;
50
51 /**
52 * Constructs a @c param_type object from the "alpha" and "beta"
53 * parameters.
54 *
55 * Requires: alpha > 0 && beta > 0
56 */
57 param_type(const RealType& alpha_arg = RealType(1.0),
58 const RealType& beta_arg = RealType(1.0))
59 : _alpha(alpha_arg), _beta(beta_arg)
60 {
61 }
62
63 /** Returns the "alpha" parameter of the distribution. */
64 RealType alpha() const { return _alpha; }
65 /** Returns the "beta" parameter of the distribution. */
66 RealType beta() const { return _beta; }
67
68#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
69 /** Writes the parameters to a @c std::ostream. */
70 template<class CharT, class Traits>
71 friend std::basic_ostream<CharT, Traits>&
72 operator<<(std::basic_ostream<CharT, Traits>& os,
73 const param_type& parm)
74 {
75 os << parm._alpha << ' ' << parm._beta;
76 return os;
77 }
78
79 /** Reads the parameters from a @c std::istream. */
80 template<class CharT, class Traits>
81 friend std::basic_istream<CharT, Traits>&
82 operator>>(std::basic_istream<CharT, Traits>& is, param_type& parm)
83 {
84 is >> parm._alpha >> std::ws >> parm._beta;
85 return is;
86 }
87#endif
88
89 /** Returns true if the two sets of parameters are the same. */
90 friend bool operator==(const param_type& lhs, const param_type& rhs)
91 {
92 return lhs._alpha == rhs._alpha && lhs._beta == rhs._beta;
93 }
94 /** Returns true if the two sets fo parameters are different. */
95 friend bool operator!=(const param_type& lhs, const param_type& rhs)
96 {
97 return !(lhs == rhs);
98 }
99 private:
100 RealType _alpha;
101 RealType _beta;
102 };
103
104#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
105 BOOST_STATIC_ASSERT(!std::numeric_limits<RealType>::is_integer);
106#endif
107
108 /**
109 * Creates a new gamma_distribution with parameters "alpha" and "beta".
110 *
111 * Requires: alpha > 0 && beta > 0
112 */
113 explicit gamma_distribution(const result_type& alpha_arg = result_type(1.0),
114 const result_type& beta_arg = result_type(1.0))
115 : _exp(result_type(1)), _alpha(alpha_arg), _beta(beta_arg)
116 {
117 BOOST_ASSERT(_alpha > result_type(0));
118 BOOST_ASSERT(_beta > result_type(0));
119 init();
120 }
121
122 /** Constructs a @c gamma_distribution from its parameters. */
123 explicit gamma_distribution(const param_type& parm)
124 : _exp(result_type(1)), _alpha(parm.alpha()), _beta(parm.beta())
125 {
126 init();
127 }
128
129 // compiler-generated copy ctor and assignment operator are fine
130
131 /** Returns the "alpha" paramter of the distribution. */
132 RealType alpha() const { return _alpha; }
133 /** Returns the "beta" parameter of the distribution. */
134 RealType beta() const { return _beta; }
135 /** Returns the smallest value that the distribution can produce. */
136 RealType min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return 0; }
137 /* Returns the largest value that the distribution can produce. */
138 RealType max BOOST_PREVENT_MACRO_SUBSTITUTION () const
139 { return (std::numeric_limits<RealType>::infinity)(); }
140
141 /** Returns the parameters of the distribution. */
142 param_type param() const { return param_type(_alpha, _beta); }
143 /** Sets the parameters of the distribution. */
144 void param(const param_type& parm)
145 {
146 _alpha = parm.alpha();
147 _beta = parm.beta();
148 init();
149 }
150
151 /**
152 * Effects: Subsequent uses of the distribution do not depend
153 * on values produced by any engine prior to invoking reset.
154 */
155 void reset() { _exp.reset(); }
156
157 /**
158 * Returns a random variate distributed according to
159 * the gamma distribution.
160 */
161 template<class Engine>
162 result_type operator()(Engine& eng)
163 {
164#ifndef BOOST_NO_STDC_NAMESPACE
165 // allow for Koenig lookup
166 using std::tan; using std::sqrt; using std::exp; using std::log;
167 using std::pow;
168#endif
169 if(_alpha == result_type(1)) {
170 return _exp(eng) * _beta;
171 } else if(_alpha > result_type(1)) {
172 // Can we have a boost::mathconst please?
173 const result_type pi = result_type(3.14159265358979323846);
174 for(;;) {
175 result_type y = tan(pi * uniform_01<RealType>()(eng));
176 result_type x = sqrt(result_type(2)*_alpha-result_type(1))*y
177 + _alpha-result_type(1);
178 if(x <= result_type(0))
179 continue;
180 if(uniform_01<RealType>()(eng) >
181 (result_type(1)+y*y) * exp((_alpha-result_type(1))
182 *log(x/(_alpha-result_type(1)))
183 - sqrt(result_type(2)*_alpha
184 -result_type(1))*y))
185 continue;
186 return x * _beta;
187 }
188 } else /* alpha < 1.0 */ {
189 for(;;) {
190 result_type u = uniform_01<RealType>()(eng);
191 result_type y = _exp(eng);
192 result_type x, q;
193 if(u < _p) {
194 x = exp(-y/_alpha);
195 q = _p*exp(-x);
196 } else {
197 x = result_type(1)+y;
198 q = _p + (result_type(1)-_p) * pow(x,_alpha-result_type(1));
199 }
200 if(u >= q)
201 continue;
202 return x * _beta;
203 }
204 }
205 }
206
207 template<class URNG>
208 RealType operator()(URNG& urng, const param_type& parm) const
209 {
210 return gamma_distribution(parm)(urng);
211 }
212
213#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
214 /** Writes a @c gamma_distribution to a @c std::ostream. */
215 template<class CharT, class Traits>
216 friend std::basic_ostream<CharT,Traits>&
217 operator<<(std::basic_ostream<CharT,Traits>& os,
218 const gamma_distribution& gd)
219 {
220 os << gd.param();
221 return os;
222 }
223
224 /** Reads a @c gamma_distribution from a @c std::istream. */
225 template<class CharT, class Traits>
226 friend std::basic_istream<CharT,Traits>&
227 operator>>(std::basic_istream<CharT,Traits>& is, gamma_distribution& gd)
228 {
229 gd.read(is);
230 return is;
231 }
232#endif
233
234 /**
235 * Returns true if the two distributions will produce identical
236 * sequences of random variates given equal generators.
237 */
238 friend bool operator==(const gamma_distribution& lhs,
239 const gamma_distribution& rhs)
240 {
241 return lhs._alpha == rhs._alpha
242 && lhs._beta == rhs._beta
243 && lhs._exp == rhs._exp;
244 }
245
246 /**
247 * Returns true if the two distributions can produce different
248 * sequences of random variates, given equal generators.
249 */
250 friend bool operator!=(const gamma_distribution& lhs,
251 const gamma_distribution& rhs)
252 {
253 return !(lhs == rhs);
254 }
255
256private:
257 /// \cond hide_private_members
258
259 template<class CharT, class Traits>
260 void read(std::basic_istream<CharT, Traits>& is)
261 {
262 param_type parm;
263 if(is >> parm) {
264 param(parm);
265 }
266 }
267
268 void init()
269 {
270#ifndef BOOST_NO_STDC_NAMESPACE
271 // allow for Koenig lookup
272 using std::exp;
273#endif
274 _p = exp(result_type(1)) / (_alpha + exp(result_type(1)));
275 }
276 /// \endcond
277
278 exponential_distribution<RealType> _exp;
279 result_type _alpha;
280 result_type _beta;
281 // some data precomputed from the parameters
282 result_type _p;
283};
284
285
286} // namespace random
287
288using random::gamma_distribution;
289
290} // namespace boost
291
292#endif // BOOST_RANDOM_GAMMA_DISTRIBUTION_HPP
293