1 | /* boost random/gamma_distribution.hpp header file |
---|
2 | * |
---|
3 | * Copyright Jens Maurer 2002 |
---|
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) |
---|
7 | * |
---|
8 | * See http://www.boost.org for most recent version including documentation. |
---|
9 | * |
---|
10 | * $Id: gamma_distribution.hpp,v 1.9 2004/07/27 03:43:32 dgregor Exp $ |
---|
11 | * |
---|
12 | */ |
---|
13 | |
---|
14 | #ifndef BOOST_RANDOM_GAMMA_DISTRIBUTION_HPP |
---|
15 | #define BOOST_RANDOM_GAMMA_DISTRIBUTION_HPP |
---|
16 | |
---|
17 | #include <cmath> |
---|
18 | #include <cassert> |
---|
19 | #include <boost/limits.hpp> |
---|
20 | #include <boost/static_assert.hpp> |
---|
21 | #include <boost/random/exponential_distribution.hpp> |
---|
22 | |
---|
23 | namespace boost { |
---|
24 | |
---|
25 | // Knuth |
---|
26 | template<class RealType = double> |
---|
27 | class gamma_distribution |
---|
28 | { |
---|
29 | public: |
---|
30 | typedef RealType input_type; |
---|
31 | typedef RealType result_type; |
---|
32 | |
---|
33 | #ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS |
---|
34 | BOOST_STATIC_ASSERT(!std::numeric_limits<RealType>::is_integer); |
---|
35 | #endif |
---|
36 | |
---|
37 | explicit gamma_distribution(const result_type& alpha = result_type(1)) |
---|
38 | : _exp(result_type(1)), _alpha(alpha) |
---|
39 | { |
---|
40 | assert(alpha > result_type(0)); |
---|
41 | init(); |
---|
42 | } |
---|
43 | |
---|
44 | // compiler-generated copy ctor and assignment operator are fine |
---|
45 | |
---|
46 | RealType alpha() const { return _alpha; } |
---|
47 | |
---|
48 | void reset() { _exp.reset(); } |
---|
49 | |
---|
50 | template<class Engine> |
---|
51 | result_type operator()(Engine& eng) |
---|
52 | { |
---|
53 | #ifndef BOOST_NO_STDC_NAMESPACE |
---|
54 | // allow for Koenig lookup |
---|
55 | using std::tan; using std::sqrt; using std::exp; using std::log; |
---|
56 | using std::pow; |
---|
57 | #endif |
---|
58 | if(_alpha == result_type(1)) { |
---|
59 | return _exp(eng); |
---|
60 | } else if(_alpha > result_type(1)) { |
---|
61 | // Can we have a boost::mathconst please? |
---|
62 | const result_type pi = result_type(3.14159265358979323846); |
---|
63 | for(;;) { |
---|
64 | result_type y = tan(pi * eng()); |
---|
65 | result_type x = sqrt(result_type(2)*_alpha-result_type(1))*y |
---|
66 | + _alpha-result_type(1); |
---|
67 | if(x <= result_type(0)) |
---|
68 | continue; |
---|
69 | if(eng() > |
---|
70 | (result_type(1)+y*y) * exp((_alpha-result_type(1)) |
---|
71 | *log(x/(_alpha-result_type(1))) |
---|
72 | - sqrt(result_type(2)*_alpha |
---|
73 | -result_type(1))*y)) |
---|
74 | continue; |
---|
75 | return x; |
---|
76 | } |
---|
77 | } else /* alpha < 1.0 */ { |
---|
78 | for(;;) { |
---|
79 | result_type u = eng(); |
---|
80 | result_type y = _exp(eng); |
---|
81 | result_type x, q; |
---|
82 | if(u < _p) { |
---|
83 | x = exp(-y/_alpha); |
---|
84 | q = _p*exp(-x); |
---|
85 | } else { |
---|
86 | x = result_type(1)+y; |
---|
87 | q = _p + (result_type(1)-_p) * pow(x, _alpha-result_type(1)); |
---|
88 | } |
---|
89 | if(u >= q) |
---|
90 | continue; |
---|
91 | return x; |
---|
92 | } |
---|
93 | } |
---|
94 | } |
---|
95 | |
---|
96 | #if !defined(BOOST_NO_OPERATORS_IN_NAMESPACE) && !defined(BOOST_NO_MEMBER_TEMPLATE_FRIENDS) |
---|
97 | template<class CharT, class Traits> |
---|
98 | friend std::basic_ostream<CharT,Traits>& |
---|
99 | operator<<(std::basic_ostream<CharT,Traits>& os, const gamma_distribution& gd) |
---|
100 | { |
---|
101 | os << gd._alpha; |
---|
102 | return os; |
---|
103 | } |
---|
104 | |
---|
105 | template<class CharT, class Traits> |
---|
106 | friend std::basic_istream<CharT,Traits>& |
---|
107 | operator>>(std::basic_istream<CharT,Traits>& is, gamma_distribution& gd) |
---|
108 | { |
---|
109 | is >> std::ws >> gd._alpha; |
---|
110 | gd.init(); |
---|
111 | return is; |
---|
112 | } |
---|
113 | #endif |
---|
114 | |
---|
115 | private: |
---|
116 | void init() |
---|
117 | { |
---|
118 | #ifndef BOOST_NO_STDC_NAMESPACE |
---|
119 | // allow for Koenig lookup |
---|
120 | using std::exp; |
---|
121 | #endif |
---|
122 | _p = exp(result_type(1)) / (_alpha + exp(result_type(1))); |
---|
123 | } |
---|
124 | |
---|
125 | exponential_distribution<RealType> _exp; |
---|
126 | result_type _alpha; |
---|
127 | // some data precomputed from the parameters |
---|
128 | result_type _p; |
---|
129 | }; |
---|
130 | |
---|
131 | } // namespace boost |
---|
132 | |
---|
133 | #endif // BOOST_RANDOM_GAMMA_DISTRIBUTION_HPP |
---|