Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
Boost165GammaDistribution.hpp
Go to the documentation of this file.
1
17/* boost random/gamma_distribution.hpp header file
18 *
19 * Copyright Jens Maurer 2002
20 * Copyright Steven Watanabe 2010
21 * Distributed under the Boost Software License, Version 1.0. (See
22 * accompanying file LICENSE_1_0.txt or copy at
23 * http://www.boost.org/LICENSE_1_0.txt)
24 *
25 * See http://www.boost.org for most recent version including documentation.
26 *
27 * $Id$
28 *
29 */
30
31#ifndef BOOST_165_RANDOM_GAMMA_DISTRIBUTION_HPP
32#define BOOST_165_RANDOM_GAMMA_DISTRIBUTION_HPP
33
34#include <boost/assert.hpp>
35#include <boost/config/no_tr1/cmath.hpp>
36#include <boost/limits.hpp>
37#include <boost/random/detail/config.hpp>
38//#include <boost/random/exponential_distribution.hpp> // swapped for chaste copy of exponential v1.65
39#include <boost/static_assert.hpp>
40#include <iosfwd>
41#include <istream>
43
44namespace boost
45{
46namespace random
47{
48
49 // The algorithm is taken from Knuth
50
58 template <class RealType = double>
60 {
61 public:
62 typedef RealType input_type;
63 typedef RealType result_type;
64
66 {
67 public:
69
76 param_type(const RealType& alpha_arg = RealType(1.0),
77 const RealType& beta_arg = RealType(1.0))
78 : _alpha(alpha_arg), _beta(beta_arg)
79 {
80 }
81
83 RealType alpha() const { return _alpha; }
85 RealType beta() const { return _beta; }
86
87#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
89 template <class CharT, class Traits>
90 friend std::basic_ostream<CharT, Traits>&
91 operator<<(std::basic_ostream<CharT, Traits>& os,
92 const param_type& parm)
93 {
94 os << parm._alpha << ' ' << parm._beta;
95 return os;
96 }
97
99 template <class CharT, class Traits>
100 friend std::basic_istream<CharT, Traits>&
101 operator>>(std::basic_istream<CharT, Traits>& is, param_type& parm)
102 {
103 is >> parm._alpha >> std::ws >> parm._beta;
104 return is;
105 }
106#endif
107
109 friend bool operator==(const param_type& lhs, const param_type& rhs)
110 {
111 return lhs._alpha == rhs._alpha && lhs._beta == rhs._beta;
112 }
114 friend bool operator!=(const param_type& lhs, const param_type& rhs)
115 {
116 return !(lhs == rhs);
117 }
118
119 private:
120 RealType _alpha;
121 RealType _beta;
122 };
123
124#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
125 BOOST_STATIC_ASSERT(!std::numeric_limits<RealType>::is_integer);
126#endif
127
133 explicit gamma_distribution_v165(const result_type& alpha_arg = result_type(1.0),
134 const result_type& beta_arg = result_type(1.0))
135 : _exp(result_type(1)), _alpha(alpha_arg), _beta(beta_arg)
136 {
137 BOOST_ASSERT(_alpha > result_type(0));
138 BOOST_ASSERT(_beta > result_type(0));
139 init();
140 }
141
144 : _exp(result_type(1)), _alpha(parm.alpha()), _beta(parm.beta())
145 {
146 init();
147 }
148
149 // compiler-generated copy ctor and assignment operator are fine
150
152 RealType alpha() const { return _alpha; }
154 RealType beta() const { return _beta; }
156 RealType min BOOST_PREVENT_MACRO_SUBSTITUTION() const { return 0; }
157 /* Returns the largest value that the distribution can produce. */
158 RealType max BOOST_PREVENT_MACRO_SUBSTITUTION() const
159 {
160 return (std::numeric_limits<RealType>::infinity)();
161 }
162
164 param_type param() const { return param_type(_alpha, _beta); }
166 void param(const param_type& parm)
167 {
168 _alpha = parm.alpha();
169 _beta = parm.beta();
170 init();
171 }
172
177 void reset() { _exp.reset(); }
178
183 template <class Engine>
184 result_type operator()(Engine& eng)
185 {
186#ifndef BOOST_NO_STDC_NAMESPACE
187 // allow for Koenig lookup
188 using std::tan;
189 using std::sqrt;
190 using std::exp;
191 using std::log;
192 using std::pow;
193#endif
194 if (_alpha == result_type(1))
195 {
196 return _exp(eng) * _beta;
197 }
198 else if (_alpha > result_type(1))
199 {
200 // Can we have a boost::mathconst please?
201 const result_type pi = result_type(3.14159265358979323846);
202 for (;;)
203 {
204 result_type y = tan(pi * uniform_01<RealType>()(eng));
205 result_type x = sqrt(result_type(2) * _alpha - result_type(1)) * y
206 + _alpha - result_type(1);
207 if (x <= result_type(0))
208 continue;
209 if (uniform_01<RealType>()(eng) > (result_type(1) + y * y) * exp((_alpha - result_type(1))
210 * log(x / (_alpha - result_type(1)))
211 - sqrt(result_type(2) * _alpha
212 - result_type(1))
213 * y))
214 continue;
215 return x * _beta;
216 }
217 }
218 else /* alpha < 1.0 */
219 {
220 for (;;)
221 {
222 result_type u = uniform_01<RealType>()(eng);
223 result_type y = _exp(eng);
224 result_type x, q;
225 if (u < _p)
226 {
227 x = exp(-y / _alpha);
228 q = _p * exp(-x);
229 }
230 else
231 {
232 x = result_type(1) + y;
233 q = _p + (result_type(1) - _p) * pow(x, _alpha - result_type(1));
234 }
235 if (u >= q)
236 continue;
237 return x * _beta;
238 }
239 }
240 }
241
242 template <class URNG>
243 RealType operator()(URNG& urng, const param_type& parm) const
244 {
245 return gamma_distribution_v165(parm)(urng);
246 }
247
248#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
250 template <class CharT, class Traits>
251 friend std::basic_ostream<CharT, Traits>&
252 operator<<(std::basic_ostream<CharT, Traits>& os,
253 const gamma_distribution_v165& gd)
254 {
255 os << gd.param();
256 return os;
257 }
258
260 template <class CharT, class Traits>
261 friend std::basic_istream<CharT, Traits>&
262 operator>>(std::basic_istream<CharT, Traits>& is, gamma_distribution_v165& gd)
263 {
264 gd.read(is);
265 return is;
266 }
267#endif
268
273 friend bool operator==(const gamma_distribution_v165& lhs,
274 const gamma_distribution_v165& rhs)
275 {
276 return lhs._alpha == rhs._alpha
277 && lhs._beta == rhs._beta
278 && lhs._exp == rhs._exp;
279 }
280
285 friend bool operator!=(const gamma_distribution_v165& lhs,
286 const gamma_distribution_v165& rhs)
287 {
288 return !(lhs == rhs);
289 }
290
291 private:
293
294 template <class CharT, class Traits>
295 void read(std::basic_istream<CharT, Traits>& is)
296 {
297 param_type parm;
298 if (is >> parm)
299 {
300 param(parm);
301 }
302 }
303
304 void init()
305 {
306#ifndef BOOST_NO_STDC_NAMESPACE
307 // allow for Koenig lookup
308 using std::exp;
309#endif
310 _p = exp(result_type(1)) / (_alpha + exp(result_type(1)));
311 }
313
314 exponential_distribution_v165<RealType> _exp;
315 result_type _alpha;
316 result_type _beta;
317 // some data precomputed from the parameters
318 result_type _p;
319 };
320
321} // namespace random
322
323} // namespace boost
324
325#endif // BOOST_165_RANDOM_GAMMA_DISTRIBUTION_HPP
friend std::basic_ostream< CharT, Traits > & operator<<(std::basic_ostream< CharT, Traits > &os, const param_type &parm)
friend bool operator!=(const param_type &lhs, const param_type &rhs)
friend std::basic_istream< CharT, Traits > & operator>>(std::basic_istream< CharT, Traits > &is, param_type &parm)
friend bool operator==(const param_type &lhs, const param_type &rhs)
param_type(const RealType &alpha_arg=RealType(1.0), const RealType &beta_arg=RealType(1.0))
friend bool operator==(const gamma_distribution_v165 &lhs, const gamma_distribution_v165 &rhs)
friend std::basic_ostream< CharT, Traits > & operator<<(std::basic_ostream< CharT, Traits > &os, const gamma_distribution_v165 &gd)
friend bool operator!=(const gamma_distribution_v165 &lhs, const gamma_distribution_v165 &rhs)
gamma_distribution_v165(const result_type &alpha_arg=result_type(1.0), const result_type &beta_arg=result_type(1.0))
friend std::basic_istream< CharT, Traits > & operator>>(std::basic_istream< CharT, Traits > &is, gamma_distribution_v165 &gd)