]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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 | ||
27 | namespace boost { | |
28 | namespace 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 | */ | |
39 | template<class RealType = double> | |
40 | class gamma_distribution | |
41 | { | |
42 | public: | |
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 | ||
256 | private: | |
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 | ||
288 | using random::gamma_distribution; | |
289 | ||
290 | } // namespace boost | |
291 | ||
292 | #endif // BOOST_RANDOM_GAMMA_DISTRIBUTION_HPP |