]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | /* boost random/poisson_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_POISSON_DISTRIBUTION_HPP | |
16 | #define BOOST_RANDOM_POISSON_DISTRIBUTION_HPP | |
17 | ||
18 | #include <boost/config/no_tr1/cmath.hpp> | |
19 | #include <cstdlib> | |
20 | #include <iosfwd> | |
21 | #include <boost/assert.hpp> | |
22 | #include <boost/limits.hpp> | |
23 | #include <boost/random/uniform_01.hpp> | |
24 | #include <boost/random/detail/config.hpp> | |
25 | ||
26 | #include <boost/random/detail/disable_warnings.hpp> | |
27 | ||
28 | namespace boost { | |
29 | namespace random { | |
30 | ||
31 | namespace detail { | |
32 | ||
33 | template<class RealType> | |
34 | struct poisson_table { | |
35 | static RealType value[10]; | |
36 | }; | |
37 | ||
38 | template<class RealType> | |
39 | RealType poisson_table<RealType>::value[10] = { | |
40 | 0.0, | |
41 | 0.0, | |
42 | 0.69314718055994529, | |
43 | 1.7917594692280550, | |
44 | 3.1780538303479458, | |
45 | 4.7874917427820458, | |
46 | 6.5792512120101012, | |
47 | 8.5251613610654147, | |
48 | 10.604602902745251, | |
49 | 12.801827480081469 | |
50 | }; | |
51 | ||
52 | } | |
53 | ||
54 | /** | |
55 | * An instantiation of the class template @c poisson_distribution is a | |
56 | * model of \random_distribution. The poisson distribution has | |
57 | * \f$p(i) = \frac{e^{-\lambda}\lambda^i}{i!}\f$ | |
58 | * | |
59 | * This implementation is based on the PTRD algorithm described | |
60 | * | |
61 | * @blockquote | |
62 | * "The transformed rejection method for generating Poisson random variables", | |
63 | * Wolfgang Hormann, Insurance: Mathematics and Economics | |
64 | * Volume 12, Issue 1, February 1993, Pages 39-45 | |
65 | * @endblockquote | |
66 | */ | |
67 | template<class IntType = int, class RealType = double> | |
68 | class poisson_distribution { | |
69 | public: | |
70 | typedef IntType result_type; | |
71 | typedef RealType input_type; | |
72 | ||
73 | class param_type { | |
74 | public: | |
75 | typedef poisson_distribution distribution_type; | |
76 | /** | |
77 | * Construct a param_type object with the parameter "mean" | |
78 | * | |
79 | * Requires: mean > 0 | |
80 | */ | |
81 | explicit param_type(RealType mean_arg = RealType(1)) | |
82 | : _mean(mean_arg) | |
83 | { | |
84 | BOOST_ASSERT(_mean > 0); | |
85 | } | |
86 | /* Returns the "mean" parameter of the distribution. */ | |
87 | RealType mean() const { return _mean; } | |
88 | #ifndef BOOST_RANDOM_NO_STREAM_OPERATORS | |
89 | /** Writes the parameters of the distribution to a @c std::ostream. */ | |
90 | template<class CharT, class Traits> | |
91 | friend std::basic_ostream<CharT, Traits>& | |
92 | operator<<(std::basic_ostream<CharT, Traits>& os, | |
93 | const param_type& parm) | |
94 | { | |
95 | os << parm._mean; | |
96 | return os; | |
97 | } | |
98 | ||
99 | /** Reads the parameters of the distribution from a @c std::istream. */ | |
100 | template<class CharT, class Traits> | |
101 | friend std::basic_istream<CharT, Traits>& | |
102 | operator>>(std::basic_istream<CharT, Traits>& is, param_type& parm) | |
103 | { | |
104 | is >> parm._mean; | |
105 | return is; | |
106 | } | |
107 | #endif | |
108 | /** Returns true if the parameters have the same values. */ | |
109 | friend bool operator==(const param_type& lhs, const param_type& rhs) | |
110 | { | |
111 | return lhs._mean == rhs._mean; | |
112 | } | |
113 | /** Returns true if the parameters have different values. */ | |
114 | friend bool operator!=(const param_type& lhs, const param_type& rhs) | |
115 | { | |
116 | return !(lhs == rhs); | |
117 | } | |
118 | private: | |
119 | RealType _mean; | |
120 | }; | |
121 | ||
122 | /** | |
123 | * Constructs a @c poisson_distribution with the parameter @c mean. | |
124 | * | |
125 | * Requires: mean > 0 | |
126 | */ | |
127 | explicit poisson_distribution(RealType mean_arg = RealType(1)) | |
128 | : _mean(mean_arg) | |
129 | { | |
130 | BOOST_ASSERT(_mean > 0); | |
131 | init(); | |
132 | } | |
133 | ||
134 | /** | |
135 | * Construct an @c poisson_distribution object from the | |
136 | * parameters. | |
137 | */ | |
138 | explicit poisson_distribution(const param_type& parm) | |
139 | : _mean(parm.mean()) | |
140 | { | |
141 | init(); | |
142 | } | |
143 | ||
144 | /** | |
145 | * Returns a random variate distributed according to the | |
146 | * poisson distribution. | |
147 | */ | |
148 | template<class URNG> | |
149 | IntType operator()(URNG& urng) const | |
150 | { | |
151 | if(use_inversion()) { | |
152 | return invert(urng); | |
153 | } else { | |
154 | return generate(urng); | |
155 | } | |
156 | } | |
157 | ||
158 | /** | |
159 | * Returns a random variate distributed according to the | |
160 | * poisson distribution with parameters specified by param. | |
161 | */ | |
162 | template<class URNG> | |
163 | IntType operator()(URNG& urng, const param_type& parm) const | |
164 | { | |
165 | return poisson_distribution(parm)(urng); | |
166 | } | |
167 | ||
168 | /** Returns the "mean" parameter of the distribution. */ | |
169 | RealType mean() const { return _mean; } | |
170 | ||
171 | /** Returns the smallest value that the distribution can produce. */ | |
172 | IntType min BOOST_PREVENT_MACRO_SUBSTITUTION() const { return 0; } | |
173 | /** Returns the largest value that the distribution can produce. */ | |
174 | IntType max BOOST_PREVENT_MACRO_SUBSTITUTION() const | |
175 | { return (std::numeric_limits<IntType>::max)(); } | |
176 | ||
177 | /** Returns the parameters of the distribution. */ | |
178 | param_type param() const { return param_type(_mean); } | |
179 | /** Sets parameters of the distribution. */ | |
180 | void param(const param_type& parm) | |
181 | { | |
182 | _mean = parm.mean(); | |
183 | init(); | |
184 | } | |
185 | ||
186 | /** | |
187 | * Effects: Subsequent uses of the distribution do not depend | |
188 | * on values produced by any engine prior to invoking reset. | |
189 | */ | |
190 | void reset() { } | |
191 | ||
192 | #ifndef BOOST_RANDOM_NO_STREAM_OPERATORS | |
193 | /** Writes the parameters of the distribution to a @c std::ostream. */ | |
194 | template<class CharT, class Traits> | |
195 | friend std::basic_ostream<CharT,Traits>& | |
196 | operator<<(std::basic_ostream<CharT,Traits>& os, | |
197 | const poisson_distribution& pd) | |
198 | { | |
199 | os << pd.param(); | |
200 | return os; | |
201 | } | |
202 | ||
203 | /** Reads the parameters of the distribution from a @c std::istream. */ | |
204 | template<class CharT, class Traits> | |
205 | friend std::basic_istream<CharT,Traits>& | |
206 | operator>>(std::basic_istream<CharT,Traits>& is, poisson_distribution& pd) | |
207 | { | |
208 | pd.read(is); | |
209 | return is; | |
210 | } | |
211 | #endif | |
212 | ||
213 | /** Returns true if the two distributions will produce the same | |
214 | sequence of values, given equal generators. */ | |
215 | friend bool operator==(const poisson_distribution& lhs, | |
216 | const poisson_distribution& rhs) | |
217 | { | |
218 | return lhs._mean == rhs._mean; | |
219 | } | |
220 | /** Returns true if the two distributions could produce different | |
221 | sequences of values, given equal generators. */ | |
222 | friend bool operator!=(const poisson_distribution& lhs, | |
223 | const poisson_distribution& rhs) | |
224 | { | |
225 | return !(lhs == rhs); | |
226 | } | |
227 | ||
228 | private: | |
229 | ||
230 | /// @cond show_private | |
231 | ||
232 | template<class CharT, class Traits> | |
233 | void read(std::basic_istream<CharT, Traits>& is) { | |
234 | param_type parm; | |
235 | if(is >> parm) { | |
236 | param(parm); | |
237 | } | |
238 | } | |
239 | ||
240 | bool use_inversion() const | |
241 | { | |
242 | return _mean < 10; | |
243 | } | |
244 | ||
245 | static RealType log_factorial(IntType k) | |
246 | { | |
247 | BOOST_ASSERT(k >= 0); | |
248 | BOOST_ASSERT(k < 10); | |
249 | return detail::poisson_table<RealType>::value[k]; | |
250 | } | |
251 | ||
252 | void init() | |
253 | { | |
254 | using std::sqrt; | |
255 | using std::exp; | |
256 | ||
257 | if(use_inversion()) { | |
258 | _exp_mean = exp(-_mean); | |
259 | } else { | |
260 | _ptrd.smu = sqrt(_mean); | |
261 | _ptrd.b = 0.931 + 2.53 * _ptrd.smu; | |
262 | _ptrd.a = -0.059 + 0.02483 * _ptrd.b; | |
263 | _ptrd.inv_alpha = 1.1239 + 1.1328 / (_ptrd.b - 3.4); | |
264 | _ptrd.v_r = 0.9277 - 3.6224 / (_ptrd.b - 2); | |
265 | } | |
266 | } | |
267 | ||
268 | template<class URNG> | |
269 | IntType generate(URNG& urng) const | |
270 | { | |
271 | using std::floor; | |
272 | using std::abs; | |
273 | using std::log; | |
274 | ||
275 | while(true) { | |
276 | RealType u; | |
277 | RealType v = uniform_01<RealType>()(urng); | |
278 | if(v <= 0.86 * _ptrd.v_r) { | |
279 | u = v / _ptrd.v_r - 0.43; | |
280 | return static_cast<IntType>(floor( | |
281 | (2*_ptrd.a/(0.5-abs(u)) + _ptrd.b)*u + _mean + 0.445)); | |
282 | } | |
283 | ||
284 | if(v >= _ptrd.v_r) { | |
285 | u = uniform_01<RealType>()(urng) - 0.5; | |
286 | } else { | |
287 | u = v/_ptrd.v_r - 0.93; | |
288 | u = ((u < 0)? -0.5 : 0.5) - u; | |
289 | v = uniform_01<RealType>()(urng) * _ptrd.v_r; | |
290 | } | |
291 | ||
292 | RealType us = 0.5 - abs(u); | |
293 | if(us < 0.013 && v > us) { | |
294 | continue; | |
295 | } | |
296 | ||
297 | RealType k = floor((2*_ptrd.a/us + _ptrd.b)*u+_mean+0.445); | |
298 | v = v*_ptrd.inv_alpha/(_ptrd.a/(us*us) + _ptrd.b); | |
299 | ||
300 | RealType log_sqrt_2pi = 0.91893853320467267; | |
301 | ||
302 | if(k >= 10) { | |
303 | if(log(v*_ptrd.smu) <= (k + 0.5)*log(_mean/k) | |
304 | - _mean | |
305 | - log_sqrt_2pi | |
306 | + k | |
307 | - (1/12. - (1/360. - 1/(1260.*k*k))/(k*k))/k) { | |
308 | return static_cast<IntType>(k); | |
309 | } | |
310 | } else if(k >= 0) { | |
311 | if(log(v) <= k*log(_mean) | |
312 | - _mean | |
313 | - log_factorial(static_cast<IntType>(k))) { | |
314 | return static_cast<IntType>(k); | |
315 | } | |
316 | } | |
317 | } | |
318 | } | |
319 | ||
320 | template<class URNG> | |
321 | IntType invert(URNG& urng) const | |
322 | { | |
323 | RealType p = _exp_mean; | |
324 | IntType x = 0; | |
325 | RealType u = uniform_01<RealType>()(urng); | |
326 | while(u > p) { | |
327 | u = u - p; | |
328 | ++x; | |
329 | p = _mean * p / x; | |
330 | } | |
331 | return x; | |
332 | } | |
333 | ||
334 | RealType _mean; | |
335 | ||
336 | union { | |
337 | // for ptrd | |
338 | struct { | |
339 | RealType v_r; | |
340 | RealType a; | |
341 | RealType b; | |
342 | RealType smu; | |
343 | RealType inv_alpha; | |
344 | } _ptrd; | |
345 | // for inversion | |
346 | RealType _exp_mean; | |
347 | }; | |
348 | ||
349 | /// @endcond | |
350 | }; | |
351 | ||
352 | } // namespace random | |
353 | ||
354 | using random::poisson_distribution; | |
355 | ||
356 | } // namespace boost | |
357 | ||
358 | #include <boost/random/detail/enable_warnings.hpp> | |
359 | ||
360 | #endif // BOOST_RANDOM_POISSON_DISTRIBUTION_HPP |