boost/boost/random/poisson_distribution.hpp - nest-cam/4320010/boost - Git at Google

 /* boost random/poisson_distribution.hpp header file
  *
  * Copyright Jens Maurer 2002
  * Copyright Steven Watanabe 2010
  * Distributed under the Boost Software License, Version 1.0. (See
  * accompanying file LICENSE_1_0.txt or copy at
  * http://www.boost.org/LICENSE_1_0.txt)
  *
  * See http://www.boost.org for most recent version including documentation.
  *
  * $Id$
  *
  */

 #ifndef BOOST_RANDOM_POISSON_DISTRIBUTION_HPP
 #define BOOST_RANDOM_POISSON_DISTRIBUTION_HPP

 #include <boost/config/no_tr1/cmath.hpp>
 #include <cstdlib>
 #include <iosfwd>
 #include <boost/assert.hpp>
 #include <boost/limits.hpp>
 #include <boost/random/uniform_01.hpp>
 #include <boost/random/detail/config.hpp>

 #include <boost/random/detail/disable_warnings.hpp>

 namespace boost {
 namespace random {

 namespace detail {

 template<class RealType>
 struct poisson_table {
     static RealType value[10];
 };

 template<class RealType>
 RealType poisson_table<RealType>::value[10] = {
     0.0,
     0.0,
     0.69314718055994529,
     1.7917594692280550,
     3.1780538303479458,
     4.7874917427820458,
     6.5792512120101012,
     8.5251613610654147,
     10.604602902745251,
     12.801827480081469
 };

 }

 /**
  * An instantiation of the class template @c poisson_distribution is a
  * model of \random_distribution.  The poisson distribution has
  * \f$p(i) = \frac{e^{-\lambda}\lambda^i}{i!}\f$
  *
  * This implementation is based on the PTRD algorithm described
  *
  *  @blockquote
  *  "The transformed rejection method for generating Poisson random variables",
  *  Wolfgang Hormann, Insurance: Mathematics and Economics
  *  Volume 12, Issue 1, February 1993, Pages 39-45
  *  @endblockquote
  */
 template<class IntType = int, class RealType = double>
 class poisson_distribution {
 public:
     typedef IntType result_type;
     typedef RealType input_type;

     class param_type {
     public:
         typedef poisson_distribution distribution_type;
         /**
          * Construct a param_type object with the parameter "mean"
          *
          * Requires: mean > 0
          */
         explicit param_type(RealType mean_arg = RealType(1))
           : _mean(mean_arg)
         {
             BOOST_ASSERT(_mean > 0);
         }
         /* Returns the "mean" parameter of the distribution. */
         RealType mean() const { return _mean; }
 #ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
         /** Writes the parameters of the distribution to a @c std::ostream. */
         template<class CharT, class Traits>
         friend std::basic_ostream<CharT, Traits>&
         operator<<(std::basic_ostream<CharT, Traits>& os,
                    const param_type& parm)
         {
             os << parm._mean;
             return os;
         }

         /** Reads the parameters of the distribution from a @c std::istream. */
         template<class CharT, class Traits>
         friend std::basic_istream<CharT, Traits>&
         operator>>(std::basic_istream<CharT, Traits>& is, param_type& parm)
         {
             is >> parm._mean;
             return is;
         }
 #endif
         /** Returns true if the parameters have the same values. */
         friend bool operator==(const param_type& lhs, const param_type& rhs)
         {
             return lhs._mean == rhs._mean;
         }
         /** Returns true if the parameters have different values. */
         friend bool operator!=(const param_type& lhs, const param_type& rhs)
         {
             return !(lhs == rhs);
         }
     private:
         RealType _mean;
     };

     /**
      * Constructs a @c poisson_distribution with the parameter @c mean.
      *
      * Requires: mean > 0
      */
     explicit poisson_distribution(RealType mean_arg = RealType(1))
       : _mean(mean_arg)
     {
         BOOST_ASSERT(_mean > 0);
         init();
     }

     /**
      * Construct an @c poisson_distribution object from the
      * parameters.
      */
     explicit poisson_distribution(const param_type& parm)
       : _mean(parm.mean())
     {
         init();
     }

     /**
      * Returns a random variate distributed according to the
      * poisson distribution.
      */
     template<class URNG>
     IntType operator()(URNG& urng) const
     {
         if(use_inversion()) {
             return invert(urng);
         } else {
             return generate(urng);
         }
     }

     /**
      * Returns a random variate distributed according to the
      * poisson distribution with parameters specified by param.
      */
     template<class URNG>
     IntType operator()(URNG& urng, const param_type& parm) const
     {
         return poisson_distribution(parm)(urng);
     }

     /** Returns the "mean" parameter of the distribution. */
     RealType mean() const { return _mean; }

     /** Returns the smallest value that the distribution can produce. */
     IntType min BOOST_PREVENT_MACRO_SUBSTITUTION() const { return 0; }
     /** Returns the largest value that the distribution can produce. */
     IntType max BOOST_PREVENT_MACRO_SUBSTITUTION() const
     { return (std::numeric_limits<IntType>::max)(); }

     /** Returns the parameters of the distribution. */
     param_type param() const { return param_type(_mean); }
     /** Sets parameters of the distribution. */
     void param(const param_type& parm)
     {
         _mean = parm.mean();
         init();
     }

     /**
      * Effects: Subsequent uses of the distribution do not depend
      * on values produced by any engine prior to invoking reset.
      */
     void reset() { }

 #ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
     /** Writes the parameters of the distribution to a @c std::ostream. */
     template<class CharT, class Traits>
     friend std::basic_ostream<CharT,Traits>&
     operator<<(std::basic_ostream<CharT,Traits>& os,
                const poisson_distribution& pd)
     {
         os << pd.param();
         return os;
     }

     /** Reads the parameters of the distribution from a @c std::istream. */
     template<class CharT, class Traits>
     friend std::basic_istream<CharT,Traits>&
     operator>>(std::basic_istream<CharT,Traits>& is, poisson_distribution& pd)
     {
         pd.read(is);
         return is;
     }
 #endif

     /** Returns true if the two distributions will produce the same
         sequence of values, given equal generators. */
     friend bool operator==(const poisson_distribution& lhs,
                            const poisson_distribution& rhs)
     {
         return lhs._mean == rhs._mean;
     }
     /** Returns true if the two distributions could produce different
         sequences of values, given equal generators. */
     friend bool operator!=(const poisson_distribution& lhs,
                            const poisson_distribution& rhs)
     {
         return !(lhs == rhs);
     }

 private:

     /// @cond show_private

     template<class CharT, class Traits>
     void read(std::basic_istream<CharT, Traits>& is) {
         param_type parm;
         if(is >> parm) {
             param(parm);
         }
     }

     bool use_inversion() const
     {
         return _mean < 10;
     }

     static RealType log_factorial(IntType k)
     {
         BOOST_ASSERT(k >= 0);
         BOOST_ASSERT(k < 10);
         return detail::poisson_table<RealType>::value[k];
     }

     void init()
     {
         using std::sqrt;
         using std::exp;

         if(use_inversion()) {
             _exp_mean = exp(-_mean);
         } else {
             _ptrd.smu = sqrt(_mean);
             _ptrd.b = 0.931 + 2.53 * _ptrd.smu;
             _ptrd.a = -0.059 + 0.02483 * _ptrd.b;
             _ptrd.inv_alpha = 1.1239 + 1.1328 / (_ptrd.b - 3.4);
             _ptrd.v_r = 0.9277 - 3.6224 / (_ptrd.b - 2);
         }
     }

     template<class URNG>
     IntType generate(URNG& urng) const
     {
         using std::floor;
         using std::abs;
         using std::log;

         while(true) {
             RealType u;
             RealType v = uniform_01<RealType>()(urng);
             if(v <= 0.86 * _ptrd.v_r) {
                 u = v / _ptrd.v_r - 0.43;
                 return static_cast<IntType>(floor(
                     (2*_ptrd.a/(0.5-abs(u)) + _ptrd.b)*u + _mean + 0.445));
             }

             if(v >= _ptrd.v_r) {
                 u = uniform_01<RealType>()(urng) - 0.5;
             } else {
                 u = v/_ptrd.v_r - 0.93;
                 u = ((u < 0)? -0.5 : 0.5) - u;
                 v = uniform_01<RealType>()(urng) * _ptrd.v_r;
             }

             RealType us = 0.5 - abs(u);
             if(us < 0.013 && v > us) {
                 continue;
             }

             RealType k = floor((2*_ptrd.a/us + _ptrd.b)*u+_mean+0.445);
             v = v*_ptrd.inv_alpha/(_ptrd.a/(us*us) + _ptrd.b);

             RealType log_sqrt_2pi = 0.91893853320467267;

             if(k >= 10) {
                 if(log(v*_ptrd.smu) <= (k + 0.5)*log(_mean/k)
                                - _mean
                                - log_sqrt_2pi
                                + k
                                - (1/12. - (1/360. - 1/(1260.*k*k))/(k*k))/k) {
                     return static_cast<IntType>(k);
                 }
             } else if(k >= 0) {
                 if(log(v) <= k*log(_mean)
                            - _mean
                            - log_factorial(static_cast<IntType>(k))) {
                     return static_cast<IntType>(k);
                 }
             }
         }
     }

     template<class URNG>
     IntType invert(URNG& urng) const
     {
         RealType p = _exp_mean;
         IntType x = 0;
         RealType u = uniform_01<RealType>()(urng);
         while(u > p) {
             u = u - p;
             ++x;
             p = _mean * p / x;
         }
         return x;
     }

     RealType _mean;

     union {
         // for ptrd
         struct {
             RealType v_r;
             RealType a;
             RealType b;
             RealType smu;
             RealType inv_alpha;
         } _ptrd;
         // for inversion
         RealType _exp_mean;
     };

     /// @endcond
 };

 } // namespace random

 using random::poisson_distribution;

 } // namespace boost

 #include <boost/random/detail/enable_warnings.hpp>

 #endif // BOOST_RANDOM_POISSON_DISTRIBUTION_HPP
	/* boost random/poisson_distribution.hpp header file
	*
	* Copyright Jens Maurer 2002
	* Copyright Steven Watanabe 2010
	* Distributed under the Boost Software License, Version 1.0. (See
	* accompanying file LICENSE_1_0.txt or copy at
	* http://www.boost.org/LICENSE_1_0.txt)
	*
	* See http://www.boost.org for most recent version including documentation.
	*
	* $Id$
	*
	*/

	#ifndef BOOST_RANDOM_POISSON_DISTRIBUTION_HPP
	#define BOOST_RANDOM_POISSON_DISTRIBUTION_HPP

	#include <boost/config/no_tr1/cmath.hpp>
	#include <cstdlib>
	#include <iosfwd>
	#include <boost/assert.hpp>
	#include <boost/limits.hpp>
	#include <boost/random/uniform_01.hpp>
	#include <boost/random/detail/config.hpp>

	#include <boost/random/detail/disable_warnings.hpp>

	namespace boost {
	namespace random {

	namespace detail {

	template<class RealType>
	struct poisson_table {
	static RealType value[10];
	};

	template<class RealType>
	RealType poisson_table<RealType>::value[10] = {
	0.0,
	0.0,
	0.69314718055994529,
	1.7917594692280550,
	3.1780538303479458,
	4.7874917427820458,
	6.5792512120101012,
	8.5251613610654147,
	10.604602902745251,
	12.801827480081469
	};

	}

	/**
	* An instantiation of the class template @c poisson_distribution is a
	* model of \random_distribution. The poisson distribution has
	* \f$p(i) = \frac{e^{-\lambda}\lambda^i}{i!}\f$
	*
	* This implementation is based on the PTRD algorithm described
	*
	* @blockquote
	* "The transformed rejection method for generating Poisson random variables",
	* Wolfgang Hormann, Insurance: Mathematics and Economics
	* Volume 12, Issue 1, February 1993, Pages 39-45
	* @endblockquote
	*/
	template<class IntType = int, class RealType = double>
	class poisson_distribution {
	public:
	typedef IntType result_type;
	typedef RealType input_type;

	class param_type {
	public:
	typedef poisson_distribution distribution_type;
	/**
	* Construct a param_type object with the parameter "mean"
	*
	* Requires: mean > 0
	*/
	explicit param_type(RealType mean_arg = RealType(1))
	: _mean(mean_arg)
	{
	BOOST_ASSERT(_mean > 0);
	}
	/* Returns the "mean" parameter of the distribution. */
	RealType mean() const { return _mean; }
	#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
	/** Writes the parameters of the distribution to a @c std::ostream. */
	template<class CharT, class Traits>
	friend std::basic_ostream<CharT, Traits>&
	operator<<(std::basic_ostream<CharT, Traits>& os,
	const param_type& parm)
	{
	os << parm._mean;
	return os;
	}

	/** Reads the parameters of the distribution from a @c std::istream. */
	template<class CharT, class Traits>
	friend std::basic_istream<CharT, Traits>&
	operator>>(std::basic_istream<CharT, Traits>& is, param_type& parm)
	{
	is >> parm._mean;
	return is;
	}
	#endif
	/** Returns true if the parameters have the same values. */
	friend bool operator==(const param_type& lhs, const param_type& rhs)
	{
	return lhs._mean == rhs._mean;
	}
	/** Returns true if the parameters have different values. */
	friend bool operator!=(const param_type& lhs, const param_type& rhs)
	{
	return !(lhs == rhs);
	}
	private:
	RealType _mean;
	};

	/**
	* Constructs a @c poisson_distribution with the parameter @c mean.
	*
	* Requires: mean > 0
	*/
	explicit poisson_distribution(RealType mean_arg = RealType(1))
	: _mean(mean_arg)
	{
	BOOST_ASSERT(_mean > 0);
	init();
	}

	/**
	* Construct an @c poisson_distribution object from the
	* parameters.
	*/
	explicit poisson_distribution(const param_type& parm)
	: _mean(parm.mean())
	{
	init();
	}

	/**
	* Returns a random variate distributed according to the
	* poisson distribution.
	*/
	template<class URNG>
	IntType operator()(URNG& urng) const
	{
	if(use_inversion()) {
	return invert(urng);
	} else {
	return generate(urng);
	}
	}

	/**
	* Returns a random variate distributed according to the
	* poisson distribution with parameters specified by param.
	*/
	template<class URNG>
	IntType operator()(URNG& urng, const param_type& parm) const
	{
	return poisson_distribution(parm)(urng);
	}

	/** Returns the "mean" parameter of the distribution. */
	RealType mean() const { return _mean; }

	/** Returns the smallest value that the distribution can produce. */
	IntType min BOOST_PREVENT_MACRO_SUBSTITUTION() const { return 0; }
	/** Returns the largest value that the distribution can produce. */
	IntType max BOOST_PREVENT_MACRO_SUBSTITUTION() const
	{ return (std::numeric_limits<IntType>::max)(); }

	/** Returns the parameters of the distribution. */
	param_type param() const { return param_type(_mean); }
	/** Sets parameters of the distribution. */
	void param(const param_type& parm)
	{
	_mean = parm.mean();
	init();
	}

	/**
	* Effects: Subsequent uses of the distribution do not depend
	* on values produced by any engine prior to invoking reset.
	*/
	void reset() { }

	#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
	/** Writes the parameters of the distribution to a @c std::ostream. */
	template<class CharT, class Traits>
	friend std::basic_ostream<CharT,Traits>&
	operator<<(std::basic_ostream<CharT,Traits>& os,
	const poisson_distribution& pd)
	{
	os << pd.param();
	return os;
	}

	/** Reads the parameters of the distribution from a @c std::istream. */
	template<class CharT, class Traits>
	friend std::basic_istream<CharT,Traits>&
	operator>>(std::basic_istream<CharT,Traits>& is, poisson_distribution& pd)
	{
	pd.read(is);
	return is;
	}
	#endif

	/** Returns true if the two distributions will produce the same
	sequence of values, given equal generators. */
	friend bool operator==(const poisson_distribution& lhs,
	const poisson_distribution& rhs)
	{
	return lhs._mean == rhs._mean;
	}
	/** Returns true if the two distributions could produce different
	sequences of values, given equal generators. */
	friend bool operator!=(const poisson_distribution& lhs,
	const poisson_distribution& rhs)
	{
	return !(lhs == rhs);
	}

	private:

	/// @cond show_private

	template<class CharT, class Traits>
	void read(std::basic_istream<CharT, Traits>& is) {
	param_type parm;
	if(is >> parm) {
	param(parm);
	}
	}

	bool use_inversion() const
	{
	return _mean < 10;
	}

	static RealType log_factorial(IntType k)
	{
	BOOST_ASSERT(k >= 0);
	BOOST_ASSERT(k < 10);
	return detail::poisson_table<RealType>::value[k];
	}

	void init()
	{
	using std::sqrt;
	using std::exp;

	if(use_inversion()) {
	_exp_mean = exp(-_mean);
	} else {
	_ptrd.smu = sqrt(_mean);
	_ptrd.b = 0.931 + 2.53 * _ptrd.smu;
	_ptrd.a = -0.059 + 0.02483 * _ptrd.b;
	_ptrd.inv_alpha = 1.1239 + 1.1328 / (_ptrd.b - 3.4);
	_ptrd.v_r = 0.9277 - 3.6224 / (_ptrd.b - 2);
	}
	}

	template<class URNG>
	IntType generate(URNG& urng) const
	{
	using std::floor;
	using std::abs;
	using std::log;

	while(true) {
	RealType u;
	RealType v = uniform_01<RealType>()(urng);
	if(v <= 0.86 * _ptrd.v_r) {
	u = v / _ptrd.v_r - 0.43;
	return static_cast<IntType>(floor(
	(2_ptrd.a/(0.5-abs(u)) + _ptrd.b)u + _mean + 0.445));
	}

	if(v >= _ptrd.v_r) {
	u = uniform_01<RealType>()(urng) - 0.5;
	} else {
	u = v/_ptrd.v_r - 0.93;
	u = ((u < 0)? -0.5 : 0.5) - u;
	v = uniform_01<RealType>()(urng) * _ptrd.v_r;
	}

	RealType us = 0.5 - abs(u);
	if(us < 0.013 && v > us) {
	continue;
	}

	RealType k = floor((2_ptrd.a/us + _ptrd.b)u+_mean+0.445);
	v = v_ptrd.inv_alpha/(_ptrd.a/(usus) + _ptrd.b);

	RealType log_sqrt_2pi = 0.91893853320467267;

	if(k >= 10) {
	if(log(v_ptrd.smu) <= (k + 0.5)log(_mean/k)
	- _mean
	- log_sqrt_2pi
	+ k
	- (1/12. - (1/360. - 1/(1260.kk))/(k*k))/k) {
	return static_cast<IntType>(k);
	}
	} else if(k >= 0) {
	if(log(v) <= k*log(_mean)
	- _mean
	- log_factorial(static_cast<IntType>(k))) {
	return static_cast<IntType>(k);
	}
	}
	}
	}

	template<class URNG>
	IntType invert(URNG& urng) const
	{
	RealType p = _exp_mean;
	IntType x = 0;
	RealType u = uniform_01<RealType>()(urng);
	while(u > p) {
	u = u - p;
	++x;
	p = _mean * p / x;
	}
	return x;
	}

	RealType _mean;

	union {
	// for ptrd
	struct {
	RealType v_r;
	RealType a;
	RealType b;
	RealType smu;
	RealType inv_alpha;
	} _ptrd;
	// for inversion
	RealType _exp_mean;
	};

	/// @endcond
	};

	} // namespace random

	using random::poisson_distribution;

	} // namespace boost

	#include <boost/random/detail/enable_warnings.hpp>

	#endif // BOOST_RANDOM_POISSON_DISTRIBUTION_HPP