From 6add8ddfef4081e35ecfb7e67a457bfab97e1fe3 Mon Sep 17 00:00:00 2001 From: Howard Hinnant Date: Sat, 15 May 2010 21:36:23 +0000 Subject: [PATCH] Revisited [rand.dist.bern.bin] and [rand.dist.pois.poisson] with better algorithms git-svn-id: https://llvm.org/svn/llvm-project/libcxx/trunk@103886 91177308-0d34-0410-b5e6-96231b3b80d8 --- include/random | 716 ++++++++++-------- .../rand.dist.bern.bernoulli/eval.pass.cpp | 54 +- .../eval_param.pass.cpp | 56 +- .../rand.dist.bern.bin/eval.pass.cpp | 231 +++++- .../rand.dist.bern.bin/eval_param.pass.cpp | 88 ++- .../rand.dist.pois.poisson/eval.pass.cpp | 2 +- .../eval_param.pass.cpp | 2 +- 7 files changed, 766 insertions(+), 383 deletions(-) diff --git a/include/random b/include/random index 4045574b..63a44aae 100644 --- a/include/random +++ b/include/random @@ -3143,11 +3143,13 @@ public: { result_type __t_; double __p_; + double __pr_; + double __odds_ratio_; + result_type __r0_; public: typedef binomial_distribution distribution_type; - explicit param_type(result_type __t = 1, double __p = 0.5) - : __t_(__t), __p_(__p) {} + explicit param_type(result_type __t = 1, double __p = 0.5); result_type t() const {return __t_;} double p() const {return __p_;} @@ -3156,6 +3158,8 @@ public: {return __x.__t_ == __y.__t_ && __x.__p_ == __y.__p_;} friend bool operator!=(const param_type& __x, const param_type& __y) {return !(__x == __y);} + + friend class binomial_distribution; }; private: @@ -3191,17 +3195,56 @@ public: {return !(__x == __y);} }; +template +binomial_distribution<_IntType>::param_type::param_type(result_type __t, double __p) + : __t_(__t), __p_(__p) +{ + if (0 < __p_ && __p_ < 1) + { + __r0_ = static_cast((__t_ + 1) * __p_); + __pr_ = _STD::exp(_STD::lgamma(__t_ + 1.) - _STD::lgamma(__r0_ + 1.) - + _STD::lgamma(__t_ - __r0_ + 1.) + __r0_ * _STD::log(__p_) + + (__t_ - __r0_) * _STD::log(1 - __p_)); + __odds_ratio_ = __p_ / (1 - __p_); + } +} + template template _IntType -binomial_distribution<_IntType>::operator()(_URNG& __g, const param_type& __p) +binomial_distribution<_IntType>::operator()(_URNG& __g, const param_type& __pr) { - bernoulli_distribution __bd(__p.p()); - _IntType __r = 0; - _IntType __t = __p.t(); - for (_IntType __i = 0; __i < __t; ++__i) - __r += __bd(__g); - return __r; + if (__pr.__t_ == 0 || __pr.__p_ == 0) + return 0; + if (__pr.__p_ == 1) + return __pr.__t_; + uniform_real_distribution __gen; + double __u = __gen(__g) - __pr.__pr_; + if (__u < 0) + return __pr.__r0_; + double __pu = __pr.__pr_; + double __pd = __pu; + result_type __ru = __pr.__r0_; + result_type __rd = __ru; + while (true) + { + if (__rd >= 1) + { + __pd *= __rd / (__pr.__odds_ratio_ * (__pr.__t_ - __rd + 1)); + __u -= __pd; + if (__u < 0) + return __rd - 1; + } + --__rd; + ++__ru; + if (__ru <= __pr.__t_) + { + __pu *= (__pr.__t_ - __ru + 1) * __pr.__odds_ratio_ / __ru; + __u -= __pu; + if (__u < 0) + return __ru; + } + } } template @@ -3234,149 +3277,6 @@ operator>>(basic_istream<_CharT, _Traits>& __is, return __is; } -// poisson_distribution - -template -class poisson_distribution -{ -public: - // types - typedef _IntType result_type; - - class param_type - { - double __mean_; - double __sq_; - double __alxm_; - double __g_; - public: - typedef poisson_distribution distribution_type; - - explicit param_type(double __mean = 1.0); - - double mean() const {return __mean_;} - - friend bool operator==(const param_type& __x, const param_type& __y) - {return __x.__mean_ == __y.__mean_;} - friend bool operator!=(const param_type& __x, const param_type& __y) - {return !(__x == __y);} - - friend class poisson_distribution; - }; - -private: - param_type __p_; - -public: - // constructors and reset functions - explicit poisson_distribution(double __mean = 1.0) : __p_(__mean) {} - explicit poisson_distribution(const param_type& __p) : __p_(__p) {} - void reset() {} - - // generating functions - template result_type operator()(_URNG& __g) - {return (*this)(__g, __p_);} - template result_type operator()(_URNG& __g, const param_type& __p); - - // property functions - double mean() const {return __p_.mean();} - - param_type param() const {return __p_;} - void param(const param_type& __p) {__p_ = __p;} - - result_type min() const {return 0;} - result_type max() const {return numeric_limits::max();} - - friend bool operator==(const poisson_distribution& __x, - const poisson_distribution& __y) - {return __x.__p_ == __y.__p_;} - friend bool operator!=(const poisson_distribution& __x, - const poisson_distribution& __y) - {return !(__x == __y);} -}; - -template -poisson_distribution<_IntType>::param_type::param_type(double __mean) - : __mean_(__mean) -{ - if (__mean_ < 12.0) - { - __g_ = _STD::exp(-__mean_); - __alxm_ = 0; - __sq_ = 0; - } - else - { - __sq_ = _STD::sqrt(2.0 * __mean_); - __alxm_ = _STD::log(__mean_); - __g_ = __mean_ * __alxm_ - _STD::lgamma(__mean_ + 1); - } -} - -template -template -_IntType -poisson_distribution<_IntType>::operator()(_URNG& __g, const param_type& __p) -{ - result_type __x; - uniform_real_distribution __gen; - if (__p.__mean_ < 12.0) - { - __x = result_type(~0); - double __t = 1; - do - { - ++__x; - __t *= __gen(__g); - } while (__t > __p.__g_); - } - else - { - double __t; - const double __pi = 3.14159265358979323846264338328; - do - { - double _X; - double __y; - do - { - __y = _STD::tan(__pi * __gen(__g)); - _X = __p.__sq_ * __y + __p.__mean_; - } while (_X < 0); - __x = static_cast(_X); - __t = 0.9 * (1 + __y * __y) * _STD::exp(__x * __p.__alxm_ - - _STD::lgamma(__x + 1.0) - __p.__g_); - } while (__gen(__g) > __t); - } - return __x; -} - -template -basic_ostream<_CharT, _Traits>& -operator<<(basic_ostream<_CharT, _Traits>& __os, - const poisson_distribution<_IntType>& __x) -{ - __save_flags<_CharT, _Traits> _(__os); - __os.flags(ios_base::dec | ios_base::left); - return __os << __x.mean(); -} - -template -basic_istream<_CharT, _Traits>& -operator>>(basic_istream<_CharT, _Traits>& __is, - poisson_distribution<_IntType>& __x) -{ - typedef poisson_distribution<_IntType> _Eng; - typedef typename _Eng::param_type param_type; - __save_flags<_CharT, _Traits> _(__is); - __is.flags(ios_base::dec | ios_base::skipws); - double __mean; - __is >> __mean; - if (!__is.fail()) - __x.param(param_type(__mean)); - return __is; -} - // exponential_distribution template @@ -3477,6 +3377,370 @@ operator>>(basic_istream<_CharT, _Traits>& __is, return __is; } +// normal_distribution + +template +class normal_distribution +{ +public: + // types + typedef _RealType result_type; + + class param_type + { + result_type __mean_; + result_type __stddev_; + public: + typedef normal_distribution distribution_type; + + explicit param_type(result_type __mean = 0, result_type __stddev = 1) + : __mean_(__mean), __stddev_(__stddev) {} + + result_type mean() const {return __mean_;} + result_type stddev() const {return __stddev_;} + + friend bool operator==(const param_type& __x, const param_type& __y) + {return __x.__mean_ == __y.__mean_ && __x.__stddev_ == __y.__stddev_;} + friend bool operator!=(const param_type& __x, const param_type& __y) + {return !(__x == __y);} + }; + +private: + param_type __p_; + result_type _V_; + bool _V_hot_; + +public: + // constructors and reset functions + explicit normal_distribution(result_type __mean = 0, result_type __stddev = 1) + : __p_(param_type(__mean, __stddev)), _V_hot_(false) {} + explicit normal_distribution(const param_type& __p) + : __p_(__p), _V_hot_(false) {} + void reset() {_V_hot_ = false;} + + // generating functions + template result_type operator()(_URNG& __g) + {return (*this)(__g, __p_);} + template result_type operator()(_URNG& __g, const param_type& __p); + + // property functions + result_type mean() const {return __p_.mean();} + result_type stddev() const {return __p_.stddev();} + + param_type param() const {return __p_;} + void param(const param_type& __p) {__p_ = __p;} + + result_type min() const {return -numeric_limits::infinity();} + result_type max() const {return numeric_limits::infinity();} + + friend bool operator==(const normal_distribution& __x, + const normal_distribution& __y) + {return __x.__p_ == __y.__p_ && __x._V_hot_ == __y._V_hot_ && + (!__x._V_hot_ || __x._V_ == __y._V_);} + friend bool operator!=(const normal_distribution& __x, + const normal_distribution& __y) + {return !(__x == __y);} + + template + friend + basic_ostream<_CharT, _Traits>& + operator<<(basic_ostream<_CharT, _Traits>& __os, + const normal_distribution<_RT>& __x); + + template + friend + basic_istream<_CharT, _Traits>& + operator>>(basic_istream<_CharT, _Traits>& __is, + normal_distribution<_RT>& __x); +}; + +template +template +_RealType +normal_distribution<_RealType>::operator()(_URNG& __g, const param_type& __p) +{ + result_type _U; + if (_V_hot_) + { + _V_hot_ = false; + _U = _V_; + } + else + { + uniform_real_distribution _Uni(-1, 1); + result_type __u; + result_type __v; + result_type __s; + do + { + __u = _Uni(__g); + __v = _Uni(__g); + __s = __u * __u + __v * __v; + } while (__s > 1 || __s == 0); + result_type _F = _STD::sqrt(-2 * _STD::log(__s) / __s); + _V_ = __v * _F; + _V_hot_ = true; + _U = __u * _F; + } + return _U * __p.stddev() + __p.mean(); +} + +template +basic_ostream<_CharT, _Traits>& +operator<<(basic_ostream<_CharT, _Traits>& __os, + const normal_distribution<_RT>& __x) +{ + __save_flags<_CharT, _Traits> _(__os); + __os.flags(ios_base::dec | ios_base::left); + _CharT __sp = __os.widen(' '); + __os.fill(__sp); + __os << __x.mean() << __sp << __x.stddev() << __sp << __x._V_hot_; + if (__x._V_hot_) + __os << __sp << __x._V_; + return __os; +} + +template +basic_istream<_CharT, _Traits>& +operator>>(basic_istream<_CharT, _Traits>& __is, + normal_distribution<_RT>& __x) +{ + typedef normal_distribution<_RT> _Eng; + typedef typename _Eng::result_type result_type; + typedef typename _Eng::param_type param_type; + __save_flags<_CharT, _Traits> _(__is); + __is.flags(ios_base::dec | ios_base::skipws); + result_type __mean; + result_type __stddev; + result_type _V = 0; + bool _V_hot = false; + __is >> __mean >> __stddev >> _V_hot; + if (_V_hot) + __is >> _V; + if (!__is.fail()) + { + __x.param(param_type(__mean, __stddev)); + __x._V_hot_ = _V_hot; + __x._V_ = _V; + } + return __is; +} + +// poisson_distribution + +template +class poisson_distribution +{ +public: + // types + typedef _IntType result_type; + + class param_type + { + double __mean_; + double __s_; + double __d_; + double __l_; + double __omega_; + double __c0_; + double __c1_; + double __c2_; + double __c3_; + double __c_; + + public: + typedef poisson_distribution distribution_type; + + explicit param_type(double __mean = 1.0); + + double mean() const {return __mean_;} + + friend bool operator==(const param_type& __x, const param_type& __y) + {return __x.__mean_ == __y.__mean_;} + friend bool operator!=(const param_type& __x, const param_type& __y) + {return !(__x == __y);} + + friend class poisson_distribution; + }; + +private: + param_type __p_; + +public: + // constructors and reset functions + explicit poisson_distribution(double __mean = 1.0) : __p_(__mean) {} + explicit poisson_distribution(const param_type& __p) : __p_(__p) {} + void reset() {} + + // generating functions + template result_type operator()(_URNG& __g) + {return (*this)(__g, __p_);} + template result_type operator()(_URNG& __g, const param_type& __p); + + // property functions + double mean() const {return __p_.mean();} + + param_type param() const {return __p_;} + void param(const param_type& __p) {__p_ = __p;} + + result_type min() const {return 0;} + result_type max() const {return numeric_limits::max();} + + friend bool operator==(const poisson_distribution& __x, + const poisson_distribution& __y) + {return __x.__p_ == __y.__p_;} + friend bool operator!=(const poisson_distribution& __x, + const poisson_distribution& __y) + {return !(__x == __y);} +}; + +template +poisson_distribution<_IntType>::param_type::param_type(double __mean) + : __mean_(__mean) +{ + if (__mean_ < 10) + { + __s_ = 0; + __d_ = 0; + __l_ = _STD::exp(-__mean_); + __omega_ = 0; + __c3_ = 0; + __c2_ = 0; + __c1_ = 0; + __c0_ = 0; + __c_ = 0; + } + else + { + __s_ = _STD::sqrt(__mean_); + __d_ = 6 * __mean_ * __mean_; + __l_ = static_cast(__mean_ - 1.1484); + __omega_ = .3989423 / __s_; + double __b1_ = .4166667E-1 / __mean_; + double __b2_ = .3 * __b1_ * __b1_; + __c3_ = .1428571 * __b1_ * __b2_; + __c2_ = __b2_ - 15. * __c3_; + __c1_ = __b1_ - 6. * __b2_ + 45. * __c3_; + __c0_ = 1. - __b1_ + 3. * __b2_ - 15. * __c3_; + __c_ = .1069 / __mean_; + } +} + +template +template +_IntType +poisson_distribution<_IntType>::operator()(_URNG& __urng, const param_type& __pr) +{ + result_type __x; + uniform_real_distribution __urd; + if (__pr.__mean_ <= 10) + { + __x = 0; + for (double __p = __urd(__urng); __p > __pr.__l_; ++__x) + __p *= __urd(__urng); + } + else + { + double __difmuk; + double __g = __pr.__mean_ + __pr.__s_ * normal_distribution()(__urng); + double __u; + if (__g > 0) + { + __x = static_cast(__g); + if (__x >= __pr.__l_) + return __x; + __difmuk = __pr.__mean_ - __x; + __u = __urd(__urng); + if (__pr.__d_ * __u >= __difmuk * __difmuk * __difmuk) + return __x; + } + exponential_distribution __edist; + for (bool __using_exp_dist = false; true; __using_exp_dist = true) + { + double __e; + if (__using_exp_dist || __g < 0) + { + double __t; + do + { + __e = __edist(__urng); + __u = __urd(__urng); + __u += __u - 1; + __t = 1.8 + (__u < 0 ? -__e : __e); + } while (__t <= -.6744); + __x = __pr.__mean_ + __pr.__s_ * __t; + __difmuk = __pr.__mean_ - __x; + __using_exp_dist = true; + } + double __px; + double __py; + if (__x < 10) + { + const result_type __fac[] = {1, 1, 2, 6, 24, 120, 720, 5040, + 40320, 362880}; + __px = -__pr.__mean_; + __py = _STD::pow(__pr.__mean_, (double)__x) / __fac[__x]; + } + else + { + double __del = .8333333E-1 / __x; + __del -= 4.8 * __del * __del * __del; + double __v = __difmuk / __x; + if (_STD::abs(__v) > 0.25) + __px = __x * _STD::log(1 + __v) - __difmuk - __del; + else + __px = __x * __v * __v * (((((((.1250060 * __v + -.1384794) * + __v + .1421878) * __v + -.1661269) * __v + .2000118) * + __v + -.2500068) * __v + .3333333) * __v + -.5) - __del; + __py = .3989423 / _STD::sqrt(__x); + } + double __r = (0.5 - __difmuk) / __pr.__s_; + double __r2 = __r * __r; + double __fx = -0.5 * __r2; + double __fy = __pr.__omega_ * (((__pr.__c3_ * __r2 + __pr.__c2_) * + __r2 + __pr.__c1_) * __r2 + __pr.__c0_); + if (__using_exp_dist) + { + if (__pr.__c_ * _STD::abs(__u) <= __py * _STD::exp(__px + __e) - + __fy * _STD::exp(__fx + __e)) + break; + } + else + { + if (__fy - __u * __fy <= __py * _STD::exp(__px - __fx)) + break; + } + } + } + return __x; +} + +template +basic_ostream<_CharT, _Traits>& +operator<<(basic_ostream<_CharT, _Traits>& __os, + const poisson_distribution<_IntType>& __x) +{ + __save_flags<_CharT, _Traits> _(__os); + __os.flags(ios_base::dec | ios_base::left); + return __os << __x.mean(); +} + +template +basic_istream<_CharT, _Traits>& +operator>>(basic_istream<_CharT, _Traits>& __is, + poisson_distribution<_IntType>& __x) +{ + typedef poisson_distribution<_IntType> _Eng; + typedef typename _Eng::param_type param_type; + __save_flags<_CharT, _Traits> _(__is); + __is.flags(ios_base::dec | ios_base::skipws); + double __mean; + __is >> __mean; + if (!__is.fail()) + __x.param(param_type(__mean)); + return __is; +} + // gamma_distribution template @@ -3629,154 +3893,6 @@ operator>>(basic_istream<_CharT, _Traits>& __is, __x.param(param_type(__alpha, __beta)); return __is; } -// normal_distribution - -template -class normal_distribution -{ -public: - // types - typedef _RealType result_type; - - class param_type - { - result_type __mean_; - result_type __stddev_; - public: - typedef normal_distribution distribution_type; - - explicit param_type(result_type __mean = 0, result_type __stddev = 1) - : __mean_(__mean), __stddev_(__stddev) {} - - result_type mean() const {return __mean_;} - result_type stddev() const {return __stddev_;} - - friend bool operator==(const param_type& __x, const param_type& __y) - {return __x.__mean_ == __y.__mean_ && __x.__stddev_ == __y.__stddev_;} - friend bool operator!=(const param_type& __x, const param_type& __y) - {return !(__x == __y);} - }; - -private: - param_type __p_; - result_type _V_; - bool _V_hot_; - -public: - // constructors and reset functions - explicit normal_distribution(result_type __mean = 0, result_type __stddev = 1) - : __p_(param_type(__mean, __stddev)), _V_hot_(false) {} - explicit normal_distribution(const param_type& __p) - : __p_(__p), _V_hot_(false) {} - void reset() {_V_hot_ = false;} - - // generating functions - template result_type operator()(_URNG& __g) - {return (*this)(__g, __p_);} - template result_type operator()(_URNG& __g, const param_type& __p); - - // property functions - result_type mean() const {return __p_.mean();} - result_type stddev() const {return __p_.stddev();} - - param_type param() const {return __p_;} - void param(const param_type& __p) {__p_ = __p;} - - result_type min() const {return -numeric_limits::infinity();} - result_type max() const {return numeric_limits::infinity();} - - friend bool operator==(const normal_distribution& __x, - const normal_distribution& __y) - {return __x.__p_ == __y.__p_ && __x._V_hot_ == __y._V_hot_ && - (!__x._V_hot_ || __x._V_ == __y._V_);} - friend bool operator!=(const normal_distribution& __x, - const normal_distribution& __y) - {return !(__x == __y);} - - template - friend - basic_ostream<_CharT, _Traits>& - operator<<(basic_ostream<_CharT, _Traits>& __os, - const normal_distribution<_RT>& __x); - - template - friend - basic_istream<_CharT, _Traits>& - operator>>(basic_istream<_CharT, _Traits>& __is, - normal_distribution<_RT>& __x); -}; - -template -template -_RealType -normal_distribution<_RealType>::operator()(_URNG& __g, const param_type& __p) -{ - result_type _U; - if (_V_hot_) - { - _V_hot_ = false; - _U = _V_; - } - else - { - uniform_real_distribution _Uni(-1, 1); - result_type __u; - result_type __v; - result_type __s; - do - { - __u = _Uni(__g); - __v = _Uni(__g); - __s = __u * __u + __v * __v; - } while (__s > 1 || __s == 0); - result_type _F = _STD::sqrt(-2 * _STD::log(__s) / __s); - _V_ = __v * _F; - _V_hot_ = true; - _U = __u * _F; - } - return _U * __p.stddev() + __p.mean(); -} - -template -basic_ostream<_CharT, _Traits>& -operator<<(basic_ostream<_CharT, _Traits>& __os, - const normal_distribution<_RT>& __x) -{ - __save_flags<_CharT, _Traits> _(__os); - __os.flags(ios_base::dec | ios_base::left); - _CharT __sp = __os.widen(' '); - __os.fill(__sp); - __os << __x.mean() << __sp << __x.stddev() << __sp << __x._V_hot_; - if (__x._V_hot_) - __os << __sp << __x._V_; - return __os; -} - -template -basic_istream<_CharT, _Traits>& -operator>>(basic_istream<_CharT, _Traits>& __is, - normal_distribution<_RT>& __x) -{ - typedef normal_distribution<_RT> _Eng; - typedef typename _Eng::result_type result_type; - typedef typename _Eng::param_type param_type; - __save_flags<_CharT, _Traits> _(__is); - __is.flags(ios_base::dec | ios_base::skipws); - result_type __mean; - result_type __stddev; - result_type _V = 0; - bool _V_hot = false; - __is >> __mean >> __stddev >> _V_hot; - if (_V_hot) - __is >> _V; - if (!__is.fail()) - { - __x.param(param_type(__mean, __stddev)); - __x._V_hot_ = _V_hot; - __x._V_ = _V; - } - return __is; -} _LIBCPP_END_NAMESPACE_STD diff --git a/test/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.bernoulli/eval.pass.cpp b/test/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.bernoulli/eval.pass.cpp index 2c82a4fc..c8e9db17 100644 --- a/test/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.bernoulli/eval.pass.cpp +++ b/test/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.bernoulli/eval.pass.cpp @@ -14,22 +14,58 @@ // template result_type operator()(_URNG& g); #include +#include +#include #include +template +inline +T +sqr(T x) +{ + return x * x; +} + int main() { { typedef std::bernoulli_distribution D; - typedef std::minstd_rand0 G; + typedef std::minstd_rand G; G g; D d(.75); - int count = 0; - for (int i = 0; i < 10000; ++i) - { - bool u = d(g); - if (u) - ++count; - } - assert(count > 7400); + const int N = 100000; + std::vector u; + for (int i = 0; i < N; ++i) + u.push_back(d(g)); + double mean = std::accumulate(u.begin(), u.end(), + double(0)) / u.size(); + double var = 0; + for (int i = 0; i < u.size(); ++i) + var += sqr(u[i] - mean); + var /= u.size(); + double x_mean = d.p(); + double x_var = d.p()*(1-d.p()); + assert(std::abs(mean - x_mean) / x_mean < 0.01); + assert(std::abs(var - x_var) / x_var < 0.01); + } + { + typedef std::bernoulli_distribution D; + typedef std::minstd_rand G; + G g; + D d(.25); + const int N = 100000; + std::vector u; + for (int i = 0; i < N; ++i) + u.push_back(d(g)); + double mean = std::accumulate(u.begin(), u.end(), + double(0)) / u.size(); + double var = 0; + for (int i = 0; i < u.size(); ++i) + var += sqr(u[i] - mean); + var /= u.size(); + double x_mean = d.p(); + double x_var = d.p()*(1-d.p()); + assert(std::abs(mean - x_mean) / x_mean < 0.01); + assert(std::abs(var - x_var) / x_var < 0.01); } } diff --git a/test/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.bernoulli/eval_param.pass.cpp b/test/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.bernoulli/eval_param.pass.cpp index d952a0fd..03eacb61 100644 --- a/test/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.bernoulli/eval_param.pass.cpp +++ b/test/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.bernoulli/eval_param.pass.cpp @@ -14,24 +14,62 @@ // template result_type operator()(_URNG& g, const param_type& parm); #include +#include +#include #include +template +inline +T +sqr(T x) +{ + return x * x; +} + int main() { { typedef std::bernoulli_distribution D; typedef D::param_type P; - typedef std::minstd_rand0 G; + typedef std::minstd_rand G; G g; D d(.75); P p(.25); - int count = 0; - for (int i = 0; i < 10000; ++i) - { - bool u = d(g, p); - if (u) - ++count; - } - assert(count < 2600); + const int N = 100000; + std::vector u; + for (int i = 0; i < N; ++i) + u.push_back(d(g, p)); + double mean = std::accumulate(u.begin(), u.end(), + double(0)) / u.size(); + double var = 0; + for (int i = 0; i < u.size(); ++i) + var += sqr(u[i] - mean); + var /= u.size(); + double x_mean = p.p(); + double x_var = p.p()*(1-p.p()); + assert(std::abs(mean - x_mean) / x_mean < 0.01); + assert(std::abs(var - x_var) / x_var < 0.01); + } + { + typedef std::bernoulli_distribution D; + typedef D::param_type P; + typedef std::minstd_rand G; + G g; + D d(.25); + P p(.75); + const int N = 100000; + std::vector u; + for (int i = 0; i < N; ++i) + u.push_back(d(g, p)); + double mean = std::accumulate(u.begin(), u.end(), + double(0)) / u.size(); + double var = 0; + for (int i = 0; i < u.size(); ++i) + var += sqr(u[i] - mean); + var /= u.size(); + double x_mean = p.p(); + double x_var = p.p()*(1-p.p()); + assert(std::abs(mean - x_mean) / x_mean < 0.01); + assert(std::abs(var - x_var) / x_var < 0.01); } } diff --git a/test/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.bin/eval.pass.cpp b/test/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.bin/eval.pass.cpp index 707c0213..95e267c8 100644 --- a/test/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.bin/eval.pass.cpp +++ b/test/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.bin/eval.pass.cpp @@ -15,53 +15,218 @@ // template result_type operator()(_URNG& g); #include +#include +#include #include +template +inline +T +sqr(T x) +{ + return x * x; +} + int main() { { typedef std::binomial_distribution<> D; - typedef D::param_type P; - typedef std::minstd_rand0 G; + typedef std::minstd_rand G; G g; - D d(16, .25); - int count = 0; - int r = 0; - for (int i = 0; i < 100; ++i) - { - D::result_type u = d(g); - r += u; - } - assert(int(r/100. + .5) == 4); + D d(5, .75); + const int N = 100000; + std::vector u; + for (int i = 0; i < N; ++i) + u.push_back(d(g)); + double mean = std::accumulate(u.begin(), u.end(), + double(0)) / u.size(); + double var = 0; + for (int i = 0; i < u.size(); ++i) + var += sqr(u[i] - mean); + var /= u.size(); + double x_mean = d.t() * d.p(); + double x_var = x_mean*(1-d.p()); + assert(std::abs(mean - x_mean) / x_mean < 0.01); + assert(std::abs(var - x_var) / x_var < 0.01); } { typedef std::binomial_distribution<> D; - typedef D::param_type P; - typedef std::minstd_rand0 G; + typedef std::minstd_rand G; G g; - D d(16, .5); - int count = 0; - int r = 0; - for (int i = 0; i < 100; ++i) - { - D::result_type u = d(g); - r += u; - } - assert(int(r/100. + .5) == 8); + D d(30, .03125); + const int N = 100000; + std::vector u; + for (int i = 0; i < N; ++i) + u.push_back(d(g)); + double mean = std::accumulate(u.begin(), u.end(), + double(0)) / u.size(); + double var = 0; + for (int i = 0; i < u.size(); ++i) + var += sqr(u[i] - mean); + var /= u.size(); + double x_mean = d.t() * d.p(); + double x_var = x_mean*(1-d.p()); + assert(std::abs(mean - x_mean) / x_mean < 0.01); + assert(std::abs(var - x_var) / x_var < 0.01); } { typedef std::binomial_distribution<> D; - typedef D::param_type P; - typedef std::minstd_rand0 G; + typedef std::minstd_rand G; G g; - D d(16, .75); - int count = 0; - int r = 0; - for (int i = 0; i < 100; ++i) - { - D::result_type u = d(g); - r += u; - } - assert(int(r/100. + .5) == 12); + D d(40, .25); + const int N = 100000; + std::vector u; + for (int i = 0; i < N; ++i) + u.push_back(d(g)); + double mean = std::accumulate(u.begin(), u.end(), + double(0)) / u.size(); + double var = 0; + for (int i = 0; i < u.size(); ++i) + var += sqr(u[i] - mean); + var /= u.size(); + double x_mean = d.t() * d.p(); + double x_var = x_mean*(1-d.p()); + assert(std::abs(mean - x_mean) / x_mean < 0.01); + assert(std::abs(var - x_var) / x_var < 0.01); + } + { + typedef std::binomial_distribution<> D; + typedef std::minstd_rand G; + G g; + D d(40, 0); + const int N = 100000; + std::vector u; + for (int i = 0; i < N; ++i) + u.push_back(d(g)); + double mean = std::accumulate(u.begin(), u.end(), + double(0)) / u.size(); + double var = 0; + for (int i = 0; i < u.size(); ++i) + var += sqr(u[i] - mean); + var /= u.size(); + double x_mean = d.t() * d.p(); + double x_var = x_mean*(1-d.p()); + assert(mean == x_mean); + assert(var == x_var); + } + { + typedef std::binomial_distribution<> D; + typedef std::minstd_rand G; + G g; + D d(40, 1); + const int N = 100000; + std::vector u; + for (int i = 0; i < N; ++i) + u.push_back(d(g)); + double mean = std::accumulate(u.begin(), u.end(), + double(0)) / u.size(); + double var = 0; + for (int i = 0; i < u.size(); ++i) + var += sqr(u[i] - mean); + var /= u.size(); + double x_mean = d.t() * d.p(); + double x_var = x_mean*(1-d.p()); + assert(mean == x_mean); + assert(var == x_var); + } + { + typedef std::binomial_distribution<> D; + typedef std::minstd_rand G; + G g; + D d(400, 0.5); + const int N = 100000; + std::vector u; + for (int i = 0; i < N; ++i) + u.push_back(d(g)); + double mean = std::accumulate(u.begin(), u.end(), + double(0)) / u.size(); + double var = 0; + for (int i = 0; i < u.size(); ++i) + var += sqr(u[i] - mean); + var /= u.size(); + double x_mean = d.t() * d.p(); + double x_var = x_mean*(1-d.p()); + assert(std::abs(mean - x_mean) / x_mean < 0.01); + assert(std::abs(var - x_var) / x_var < 0.01); + } + { + typedef std::binomial_distribution<> D; + typedef std::minstd_rand G; + G g; + D d(1, 0.5); + const int N = 100000; + std::vector u; + for (int i = 0; i < N; ++i) + u.push_back(d(g)); + double mean = std::accumulate(u.begin(), u.end(), + double(0)) / u.size(); + double var = 0; + for (int i = 0; i < u.size(); ++i) + var += sqr(u[i] - mean); + var /= u.size(); + double x_mean = d.t() * d.p(); + double x_var = x_mean*(1-d.p()); + assert(std::abs(mean - x_mean) / x_mean < 0.01); + assert(std::abs(var - x_var) / x_var < 0.01); + } + { + typedef std::binomial_distribution<> D; + typedef std::minstd_rand G; + G g; + D d(0, 0.005); + const int N = 100000; + std::vector u; + for (int i = 0; i < N; ++i) + u.push_back(d(g)); + double mean = std::accumulate(u.begin(), u.end(), + double(0)) / u.size(); + double var = 0; + for (int i = 0; i < u.size(); ++i) + var += sqr(u[i] - mean); + var /= u.size(); + double x_mean = d.t() * d.p(); + double x_var = x_mean*(1-d.p()); + assert(mean == x_mean); + assert(var == x_var); + } + { + typedef std::binomial_distribution<> D; + typedef std::minstd_rand G; + G g; + D d(0, 0); + const int N = 100000; + std::vector u; + for (int i = 0; i < N; ++i) + u.push_back(d(g)); + double mean = std::accumulate(u.begin(), u.end(), + double(0)) / u.size(); + double var = 0; + for (int i = 0; i < u.size(); ++i) + var += sqr(u[i] - mean); + var /= u.size(); + double x_mean = d.t() * d.p(); + double x_var = x_mean*(1-d.p()); + assert(mean == x_mean); + assert(var == x_var); + } + { + typedef std::binomial_distribution<> D; + typedef std::minstd_rand G; + G g; + D d(0, 1); + const int N = 100000; + std::vector u; + for (int i = 0; i < N; ++i) + u.push_back(d(g)); + double mean = std::accumulate(u.begin(), u.end(), + double(0)) / u.size(); + double var = 0; + for (int i = 0; i < u.size(); ++i) + var += sqr(u[i] - mean); + var /= u.size(); + double x_mean = d.t() * d.p(); + double x_var = x_mean*(1-d.p()); + assert(mean == x_mean); + assert(var == x_var); } } diff --git a/test/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.bin/eval_param.pass.cpp b/test/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.bin/eval_param.pass.cpp index 74139d4d..4bb9ea10 100644 --- a/test/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.bin/eval_param.pass.cpp +++ b/test/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.bin/eval_param.pass.cpp @@ -15,56 +15,84 @@ // template result_type operator()(_URNG& g, const param_type& parm); #include +#include +#include #include +template +inline +T +sqr(T x) +{ + return x * x; +} + int main() { { typedef std::binomial_distribution<> D; typedef D::param_type P; - typedef std::minstd_rand0 G; + typedef std::minstd_rand G; G g; D d(16, .75); - P p(16, .25); - int count = 0; - int r = 0; - for (int i = 0; i < 100; ++i) - { - D::result_type u = d(g, p); - r += u; - } - assert(int(r/100. + .5) == 4); + P p(5, .75); + const int N = 100000; + std::vector u; + for (int i = 0; i < N; ++i) + u.push_back(d(g, p)); + double mean = std::accumulate(u.begin(), u.end(), + double(0)) / u.size(); + double var = 0; + for (int i = 0; i < u.size(); ++i) + var += sqr(u[i] - mean); + var /= u.size(); + double x_mean = p.t() * p.p(); + double x_var = x_mean*(1-p.p()); + assert(std::abs(mean - x_mean) / x_mean < 0.01); + assert(std::abs(var - x_var) / x_var < 0.01); } { typedef std::binomial_distribution<> D; typedef D::param_type P; - typedef std::minstd_rand0 G; + typedef std::minstd_rand G; G g; D d(16, .75); - P p(16, .5); - int count = 0; - int r = 0; - for (int i = 0; i < 100; ++i) - { - D::result_type u = d(g, p); - r += u; - } - assert(int(r/100. + .5) == 8); + P p(30, .03125); + const int N = 100000; + std::vector u; + for (int i = 0; i < N; ++i) + u.push_back(d(g, p)); + double mean = std::accumulate(u.begin(), u.end(), + double(0)) / u.size(); + double var = 0; + for (int i = 0; i < u.size(); ++i) + var += sqr(u[i] - mean); + var /= u.size(); + double x_mean = p.t() * p.p(); + double x_var = x_mean*(1-p.p()); + assert(std::abs(mean - x_mean) / x_mean < 0.01); + assert(std::abs(var - x_var) / x_var < 0.01); } { typedef std::binomial_distribution<> D; typedef D::param_type P; - typedef std::minstd_rand0 G; + typedef std::minstd_rand G; G g; D d(16, .75); - P p(16, .75); - int count = 0; - int r = 0; - for (int i = 0; i < 100; ++i) - { - D::result_type u = d(g, p); - r += u; - } - assert(int(r/100. + .5) == 12); + P p(40, .25); + const int N = 100000; + std::vector u; + for (int i = 0; i < N; ++i) + u.push_back(d(g, p)); + double mean = std::accumulate(u.begin(), u.end(), + double(0)) / u.size(); + double var = 0; + for (int i = 0; i < u.size(); ++i) + var += sqr(u[i] - mean); + var /= u.size(); + double x_mean = p.t() * p.p(); + double x_var = x_mean*(1-p.p()); + assert(std::abs(mean - x_mean) / x_mean < 0.01); + assert(std::abs(var - x_var) / x_var < 0.01); } } diff --git a/test/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.poisson/eval.pass.cpp b/test/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.poisson/eval.pass.cpp index f03f8871..4c8a1fdd 100644 --- a/test/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.poisson/eval.pass.cpp +++ b/test/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.poisson/eval.pass.cpp @@ -72,7 +72,7 @@ int main() typedef std::minstd_rand G; G g; D d(20); - const int N = 10000; + const int N = 100000; std::vector u; for (int i = 0; i < N; ++i) u.push_back(d(g)); diff --git a/test/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.poisson/eval_param.pass.cpp b/test/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.poisson/eval_param.pass.cpp index b7226a3f..0c094bda 100644 --- a/test/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.poisson/eval_param.pass.cpp +++ b/test/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.poisson/eval_param.pass.cpp @@ -78,7 +78,7 @@ int main() G g; D d(2); P p(20); - const int N = 10000; + const int N = 100000; std::vector u; for (int i = 0; i < N; ++i) u.push_back(d(g, p));