Subject: [Boost-commit] svn:boost r63126 - in trunk: boost/random libs/random/doc libs/random/test
From: steven_at_[hidden]
Date: 2010-06-19 23:24:30
Author: steven_watanabe
Date: 2010-06-19 23:24:28 EDT (Sat, 19 Jun 2010)
New Revision: 63126
URL: http://svn.boost.org/trac/boost/changeset/63126
Log:
Better implementation of the poisson distribution. Fixes #1540.
Added:
trunk/libs/random/test/test_poisson.cpp (contents, props changed)
trunk/libs/random/test/test_poisson_distribution.cpp (contents, props changed)
Text files modified:
trunk/boost/random/poisson_distribution.hpp | 377 ++++++++++++++++++++++++++++++++-------
trunk/libs/random/doc/random.qbk | 2
trunk/libs/random/test/Jamfile.v2 | 2
3 files changed, 309 insertions(+), 72 deletions(-)
Modified: trunk/boost/random/poisson_distribution.hpp
==============================================================================
--- trunk/boost/random/poisson_distribution.hpp (original)
+++ trunk/boost/random/poisson_distribution.hpp 2010-06-19 23:24:28 EDT (Sat, 19 Jun 2010)
@@ -1,6 +1,7 @@
/* 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)
@@ -15,102 +16,336 @@
#define BOOST_RANDOM_POISSON_DISTRIBUTION_HPP
#include <boost/config/no_tr1/cmath.hpp>
+#include <cstdlib>
#include <cassert>
-#include <iostream>
+#include <iosfwd>
#include <boost/limits.hpp>
-#include <boost/static_assert.hpp>
+#include <boost/random/uniform_01.hpp>
#include <boost/random/detail/config.hpp>
namespace boost {
+namespace random {
-// Knuth
+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
-{
+class poisson_distribution {
public:
- typedef RealType input_type;
- typedef IntType result_type;
+ typedef IntType result_type;
+ typedef RealType input_type;
- /**
- * Constructs a @c poisson_distribution with the parameter @c mean.
- *
- * Requires: mean > 0
- */
- explicit poisson_distribution(const RealType& mean_arg = RealType(1))
- : _mean(mean_arg)
- {
-#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
- // MSVC fails BOOST_STATIC_ASSERT with std::numeric_limits at class scope
- BOOST_STATIC_ASSERT(std::numeric_limits<IntType>::is_integer);
- BOOST_STATIC_ASSERT(!std::numeric_limits<RealType>::is_integer);
+ 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)
+ {
+ 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)
+ {
+ 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 parm.
+ */
+ template<class URNG>
+ IntType operator()(URNG& urng, const param_type& parm) const
+ {
+ return poisson_distribution(parm)(urng);
+ }
- assert(_mean > RealType(0));
- init();
- }
-
- // compiler-generated copy ctor and assignment operator are fine
-
- /**
- * Returns: the "mean" parameter of the distribution.
- */
- RealType mean() const { return _mean; }
- void reset() { }
-
- template<class Engine>
- result_type operator()(Engine& eng)
- {
- // TODO: This is O(_mean), but it should be O(log(_mean)) for large _mean
- RealType product = RealType(1);
- for(result_type m = 0; ; ++m) {
- product *= eng();
- if(product <= _exp_mean)
- return m;
+ /** 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
- template<class CharT, class Traits>
- friend std::basic_ostream<CharT,Traits>&
- operator<<(std::basic_ostream<CharT,Traits>& os, const poisson_distribution& pd)
- {
- os << pd._mean;
- return os;
- }
-
- template<class CharT, class Traits>
- friend std::basic_istream<CharT,Traits>&
- operator>>(std::basic_istream<CharT,Traits>& is, poisson_distribution& pd)
- {
- is >> std::ws >> pd._mean;
- pd.init();
- return is;
- }
+ /** 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)
+ {
+ param_type parm;
+ if(is >> parm) {
+ pd.param(parm);
+ }
+ 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 hide_private_members
- void init()
- {
-#ifndef BOOST_NO_STDC_NAMESPACE
- // allow for Koenig lookup
- using std::exp;
-#endif
- _exp_mean = exp(-_mean);
- }
- /// \endcond
-
- RealType _mean;
- // some precomputed data from the parameters
- RealType _exp_mean;
+
+ /// @cond
+
+ bool use_inversion() const
+ {
+ return _mean < 10;
+ }
+
+ static RealType log_factorial(IntType k)
+ {
+ assert(k >= 0);
+ 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
#endif // BOOST_RANDOM_POISSON_DISTRIBUTION_HPP
Modified: trunk/libs/random/doc/random.qbk
==============================================================================
--- trunk/libs/random/doc/random.qbk (original)
+++ trunk/libs/random/doc/random.qbk 2010-06-19 23:24:28 EDT (Sat, 19 Jun 2010)
@@ -65,7 +65,7 @@
[def __binomial_distribution [classref boost::random::binomial_distribution binomial_distribution]]
[def __cauchy_distribution [classref boost::cauchy_distribution cauchy_distribution]]
[def __gamma_distribution [classref boost::gamma_distribution gamma_distribution]]
-[def __poisson_distribution [classref boost::poisson_distribution poisson_distribution]]
+[def __poisson_distribution [classref boost::random::poisson_distribution poisson_distribution]]
[def __geometric_distribution [classref boost::geometric_distribution geometric_distribution]]
[def __triangle_distribution [classref boost::triangle_distribution triangle_distribution]]
[def __exponential_distribution [classref boost::exponential_distribution exponential_distribution]]
Modified: trunk/libs/random/test/Jamfile.v2
==============================================================================
--- trunk/libs/random/test/Jamfile.v2 (original)
+++ trunk/libs/random/test/Jamfile.v2 2010-06-19 23:24:28 EDT (Sat, 19 Jun 2010)
@@ -43,6 +43,8 @@
run test_binomial.cpp ;
run test_binomial_distribution.cpp /boost//unit_test_framework ;
+run test_poisson.cpp ;
+run test_poisson_distribution.cpp /boost//unit_test_framework ;
# run nondet_random_speed.cpp ;
# run random_device.cpp ;
Added: trunk/libs/random/test/test_poisson.cpp
==============================================================================
--- (empty file)
+++ trunk/libs/random/test/test_poisson.cpp 2010-06-19 23:24:28 EDT (Sat, 19 Jun 2010)
@@ -0,0 +1,122 @@
+/* test_poisson.cpp
+ *
+ * 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)
+ *
+ * $Id$
+ *
+ */
+
+#include <boost/random/poisson_distribution.hpp>
+#include <boost/random/uniform_real.hpp>
+#include <boost/random/mersenne_twister.hpp>
+#include <boost/math/distributions/poisson.hpp>
+#include <boost/lexical_cast.hpp>
+#include <boost/exception/diagnostic_information.hpp>
+#include <vector>
+#include <iostream>
+#include <numeric>
+
+#include "chi_squared_test.hpp"
+
+bool do_test(double mean, long long max) {
+ std::cout << "running poisson(" << mean << ")" << " " << max << " times: " << std::flush;
+
+ int max_value = static_cast<int>(mean * 4);
+ std::vector<double> expected(max_value+1);
+ {
+ boost::math::poisson dist(mean);
+ for(int i = 0; i <= max_value; ++i) {
+ expected[i] = pdf(dist, i);
+ }
+ expected.back() += 1 - cdf(dist, max_value);
+ }
+
+ boost::random::poisson_distribution<int, double> dist(mean);
+ boost::mt19937 gen;
+ std::vector<long long> results(max_value + 1);
+ for(long long i = 0; i < max; ++i) {
+ ++results[std::min(dist(gen), max_value)];
+ }
+
+ long long sum = std::accumulate(results.begin(), results.end(), 0ll);
+ if(sum != max) {
+ std::cout << "*** Failed: incorrect total: " << sum << " ***" << std::endl;
+ return false;
+ }
+ double chsqr = chi_squared_test(results, expected, max);
+
+ bool result = chsqr < 0.99;
+ const char* err = result? "" : "*";
+ std::cout << std::setprecision(17) << chsqr << err << std::endl;
+
+ std::cout << std::setprecision(6);
+
+ return result;
+}
+
+bool do_tests(int repeat, double max_mean, long long trials) {
+ boost::mt19937 gen;
+ boost::uniform_real<> rdist(1e-15, max_mean);
+ int errors = 0;
+ for(int i = 0; i < repeat; ++i) {
+ if(!do_test(rdist(gen), trials)) {
+ ++errors;
+ }
+ }
+ if(errors != 0) {
+ std::cout << "*** " << errors << " errors detected ***" << std::endl;
+ }
+ return errors == 0;
+}
+
+int usage() {
+ std::cerr << "Usage: test_poisson -r <repeat> -m <max mean> -t <trials>" << std::endl;
+ return 2;
+}
+
+template<class T>
+bool handle_option(int& argc, char**& argv, char opt, T& value) {
+ if(argv[0][1] == opt && argc > 1) {
+ --argc;
+ ++argv;
+ value = boost::lexical_cast<T>(argv[0]);
+ return true;
+ } else {
+ return false;
+ }
+}
+
+int main(int argc, char** argv) {
+ int repeat = 10;
+ double max_mean = 100000;
+ long long trials = 1000000ll;
+
+ if(argc > 0) {
+ --argc;
+ ++argv;
+ }
+ while(argc > 0) {
+ if(argv[0][0] != '-') return usage();
+ else if(!handle_option(argc, argv, 'r', repeat)
+ && !handle_option(argc, argv, 'm', max_mean)
+ && !handle_option(argc, argv, 't', trials)) {
+ return usage();
+ }
+ --argc;
+ ++argv;
+ }
+
+ try {
+ if(do_tests(repeat, max_mean, trials)) {
+ return 0;
+ } else {
+ return EXIT_FAILURE;
+ }
+ } catch(...) {
+ std::cerr << boost::current_exception_diagnostic_information() << std::endl;
+ return EXIT_FAILURE;
+ }
+}
Added: trunk/libs/random/test/test_poisson_distribution.cpp
==============================================================================
--- (empty file)
+++ trunk/libs/random/test/test_poisson_distribution.cpp 2010-06-19 23:24:28 EDT (Sat, 19 Jun 2010)
@@ -0,0 +1,104 @@
+/* test_poisson_distribution.cpp
+ *
+ * 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)
+ *
+ * $Id$
+ *
+ */
+
+#include <boost/random/poisson_distribution.hpp>
+#include <boost/random/linear_congruential.hpp>
+#include <sstream>
+
+#define BOOST_TEST_MAIN
+#include <boost/test/unit_test.hpp>
+
+BOOST_AUTO_TEST_CASE(test_constructors) {
+ boost::random::poisson_distribution<> dist;
+ BOOST_CHECK_EQUAL(dist.mean(), 1.0);
+ boost::random::poisson_distribution<> dist_one(7.5);
+ BOOST_CHECK_EQUAL(dist_one.mean(), 7.5);
+ boost::random::poisson_distribution<> copy(dist);
+ BOOST_CHECK_EQUAL(dist, copy);
+ boost::random::poisson_distribution<> copy_one(dist_one);
+ BOOST_CHECK_EQUAL(dist_one, copy_one);
+}
+
+BOOST_AUTO_TEST_CASE(test_param) {
+ boost::random::poisson_distribution<> dist(7.5);
+ boost::random::poisson_distribution<>::param_type param = dist.param();
+ BOOST_CHECK_EQUAL(param.mean(), 7.5);
+ boost::random::poisson_distribution<> copy1(param);
+ BOOST_CHECK_EQUAL(dist, copy1);
+ boost::random::poisson_distribution<> copy2;
+ copy2.param(param);
+ BOOST_CHECK_EQUAL(dist, copy2);
+
+ boost::random::poisson_distribution<>::param_type param_copy = param;
+ BOOST_CHECK_EQUAL(param, param_copy);
+ BOOST_CHECK(param == param_copy);
+ BOOST_CHECK(!(param != param_copy));
+ boost::random::poisson_distribution<>::param_type param_default;
+ BOOST_CHECK_EQUAL(param_default.mean(), 1.0);
+ BOOST_CHECK(param != param_default);
+ BOOST_CHECK(!(param == param_default));
+ boost::random::poisson_distribution<>::param_type param_one(7.5);
+ BOOST_CHECK_EQUAL(param_one.mean(), 7.5);
+ BOOST_CHECK(param_default != param_one);
+ BOOST_CHECK(!(param_default == param_one));
+}
+
+BOOST_AUTO_TEST_CASE(test_min_max) {
+ boost::random::poisson_distribution<> dist;
+ BOOST_CHECK_EQUAL((dist.min)(), 0);
+ BOOST_CHECK_EQUAL((dist.max)(), (std::numeric_limits<int>::max)());
+ boost::random::poisson_distribution<> dist_one(7.5);
+ BOOST_CHECK_EQUAL((dist_one.min)(), 0);
+ BOOST_CHECK_EQUAL((dist_one.max)(), (std::numeric_limits<int>::max)());
+}
+
+BOOST_AUTO_TEST_CASE(test_comparison) {
+ boost::random::poisson_distribution<> dist;
+ boost::random::poisson_distribution<> dist_copy(dist);
+ boost::random::poisson_distribution<> dist_one(7.5);
+ boost::random::poisson_distribution<> dist_one_copy(dist_one);
+ BOOST_CHECK(dist == dist_copy);
+ BOOST_CHECK(!(dist != dist_copy));
+ BOOST_CHECK(dist_one == dist_one_copy);
+ BOOST_CHECK(!(dist_one != dist_one_copy));
+ BOOST_CHECK(dist != dist_one);
+ BOOST_CHECK(!(dist == dist_one));
+}
+
+BOOST_AUTO_TEST_CASE(test_streaming) {
+ boost::random::poisson_distribution<> dist(7.5);
+ std::stringstream stream;
+ stream << dist;
+ boost::random::poisson_distribution<> restored_dist;
+ stream >> restored_dist;
+ BOOST_CHECK_EQUAL(dist, restored_dist);
+}
+
+BOOST_AUTO_TEST_CASE(test_generation) {
+ boost::minstd_rand0 gen;
+ boost::random::poisson_distribution<> dist;
+ boost::random::poisson_distribution<> dist_1000(1000);
+ for(int i = 0; i < 10; ++i) {
+ // Basic sanity checks. Technically these tests are incorrect,
+ // since the range of a poisson_distribution is [0, inf), but
+ // the probability that there's an error is very small.
+ int value = dist(gen);
+ BOOST_CHECK_GE(value, 0);
+ BOOST_CHECK_LE(value, 10);
+ int value_two = dist_1000(gen);
+ BOOST_CHECK_GE(value_two, 10);
+ int value_param = dist_1000(gen, dist.param());
+ BOOST_CHECK_GE(value_param, 0);
+ BOOST_CHECK_LE(value_param, 10);
+ int value_two_param = dist(gen, dist_1000.param());
+ BOOST_CHECK_GE(value_two_param, 10);
+ }
+}