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
28namespace boost {
29namespace random {
30
31namespace detail {
32
33template<class RealType>
34struct poisson_table {
35 static RealType value[10];
36};
37
38template<class RealType>
39RealType 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 */
67template<class IntType = int, class RealType = double>
68class poisson_distribution {
69public:
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
228private:
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
354using random::poisson_distribution;
355
356} // namespace boost
357
358#include <boost/random/detail/enable_warnings.hpp>
359
360#endif // BOOST_RANDOM_POISSON_DISTRIBUTION_HPP
361