PageRenderTime 51ms CodeModel.GetById 25ms app.highlight 22ms RepoModel.GetById 2ms app.codeStats 0ms

/Src/Dependencies/Boost/boost/random/poisson_distribution.hpp

http://hadesmem.googlecode.com/
C++ Header | 360 lines | 243 code | 45 blank | 72 comment | 19 complexity | 5ab7e1b21f5414cf3ba6691500274b70 MD5 | raw file
  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: poisson_distribution.hpp 71018 2011-04-05 21:27:52Z steven_watanabe $
 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