|
|
Boost-Commit : |
Subject: [Boost-commit] svn:boost r60130 - in sandbox/statistics/adaptive_rejection_sampling: boost/ars/test libs/ars/doc libs/ars/example
From: erwann.rogard_at_[hidden]
Date: 2010-03-03 14:03:05
Author: e_r
Date: 2010-03-03 14:03:04 EST (Wed, 03 Mar 2010)
New Revision: 60130
URL: http://svn.boost.org/trac/boost/changeset/60130
Log:
m
Added:
sandbox/statistics/adaptive_rejection_sampling/libs/ars/doc/convergence_test1.txt (contents, props changed)
Text files modified:
sandbox/statistics/adaptive_rejection_sampling/boost/ars/test/gamma_distribution.hpp | 94 +++++++++---------
sandbox/statistics/adaptive_rejection_sampling/boost/ars/test/normal_distribution.hpp | 94 +++++++++---------
sandbox/statistics/adaptive_rejection_sampling/boost/ars/test/standard_distribution.hpp | 199 ++++++++++++++++++++-------------------
sandbox/statistics/adaptive_rejection_sampling/libs/ars/doc/readme.txt | 63 ------------
sandbox/statistics/adaptive_rejection_sampling/libs/ars/example/standard_distribution.cpp | 22 ++-
5 files changed, 215 insertions(+), 257 deletions(-)
Modified: sandbox/statistics/adaptive_rejection_sampling/boost/ars/test/gamma_distribution.hpp
==============================================================================
--- sandbox/statistics/adaptive_rejection_sampling/boost/ars/test/gamma_distribution.hpp (original)
+++ sandbox/statistics/adaptive_rejection_sampling/boost/ars/test/gamma_distribution.hpp 2010-03-03 14:03:04 EST (Wed, 03 Mar 2010)
@@ -19,52 +19,54 @@
namespace ars{
namespace test{
-// Same as standard_distribution but for a specific distribution
-template<typename T>
-void gamma_distribution(
- T shape, // must be > 1
- T scale,
- T init_0, // must be >0 and < init_1
- T init_1, // must be > mode = (m-1) * theta
- unsigned n1, // 1e2
- unsigned n2, // 10
- unsigned n3, // 1
- unsigned n4, // 10
- unsigned n_max_reject,
- std::ostream& out
-)
-{
-
- using namespace boost;
- using namespace math;
- using namespace assign;
- typedef double val_;
- typedef ars::constant<val_> const_;
- typedef math::gamma_distribution<val_> mdist_t;
- typedef boost::mt19937 urng_;
-
- format f("Gamma(%1%,%2%) : "); f%shape%scale;
- out << f.str();
- mdist_t mdist(shape,scale);
- urng_ urng;
-
- const val_ inf_ = const_::inf_;
-
- ars::test::standard_distribution(
- mdist,
- static_cast<val_>(0), // x_min
- inf_, // x_max
- init_0,
- init_1,
- urng,
- n1,
- n2,
- n3,
- n4,
- n_max_reject,
- out
- );
-}
+struct gamma_distribution{
+
+ // Samples from gamma distribution using adaptive rejection sampling and
+ // outputs convergence statistics
+ template<typename T>
+ static void call(
+ T shape, // must be > 1
+ T scale,
+ T init_0, // must be >0 and < init_1
+ T init_1, // must be > mode = (m-1) * theta
+ unsigned n1, // 1e2
+ unsigned n2, // 10
+ unsigned n3, // 1
+ unsigned n4, // 10
+ unsigned n_max_reject,
+ std::ostream& out
+ )
+ {
+
+ using namespace boost;
+ using namespace math;
+ using namespace assign;
+ typedef double val_;
+ typedef ars::constant<val_> const_;
+ typedef math::gamma_distribution<val_> mdist_t;
+ typedef boost::mt19937 urng_;
+
+ mdist_t mdist(shape,scale);
+ urng_ urng;
+
+ const val_ inf_ = const_::inf_;
+
+ ars::test::standard_distribution::call(
+ mdist,
+ static_cast<val_>(0), // x_min
+ inf_, // x_max
+ init_0,
+ init_1,
+ urng,
+ n1,
+ n2,
+ n3,
+ n4,
+ n_max_reject,
+ out
+ );
+ }
+};
}//test
}// ars
Modified: sandbox/statistics/adaptive_rejection_sampling/boost/ars/test/normal_distribution.hpp
==============================================================================
--- sandbox/statistics/adaptive_rejection_sampling/boost/ars/test/normal_distribution.hpp (original)
+++ sandbox/statistics/adaptive_rejection_sampling/boost/ars/test/normal_distribution.hpp 2010-03-03 14:03:04 EST (Wed, 03 Mar 2010)
@@ -19,52 +19,54 @@
namespace ars{
namespace test{
-// Same as standard_distribution but for a specific distribution
-template<typename T>
-void normal_distribution(
- T mu,
- T sigma,
- T init_0, //must be < mu
- T init_1, //must be > mu
- unsigned n1, // 1e2
- unsigned n2, // 10
- unsigned n3, // 1
- unsigned n4, // 10
- unsigned n_max_reject,
- std::ostream& out
-)
-{
-
- using namespace boost;
- using namespace math;
- using namespace assign;
- typedef double value_t;
- typedef ars::constant<value_t> const_;
- typedef math::normal_distribution<value_t> mdist_t;
- typedef boost::mt19937 urng_t;
-
- const value_t inf_ = const_::inf_;
-
- format f("N(0,%1%) : "); f%sigma;
- out << f.str();
- mdist_t mdist(mu,sigma);
- urng_t urng;
-
- standard_distribution(
- mdist,
- inf_,
- inf_,
- init_0,
- init_1,
- urng,
- n1,
- n2,
- n3,
- n4,
- n_max_reject,
- out
- );
-}
+struct normal_distribution{
+
+ // Samples from normal distribution using adaptive rejection sampling and
+ // outputs convergence statistics
+ template<typename T>
+ static void call(
+ T mu,
+ T sigma,
+ T init_0, //must be < mu
+ T init_1, //must be > mu
+ unsigned n1, // 1e2
+ unsigned n2, // 10
+ unsigned n3, // 1
+ unsigned n4, // 10
+ unsigned n_max_reject,
+ std::ostream& out
+ )
+ {
+
+ using namespace boost;
+ using namespace math;
+ using namespace assign;
+ typedef double value_t;
+ typedef ars::constant<value_t> const_;
+ typedef math::normal_distribution<value_t> mdist_t;
+ typedef boost::mt19937 urng_t;
+
+ const value_t inf_ = const_::inf_;
+
+ mdist_t mdist(mu,sigma);
+ urng_t urng;
+
+ ars::test::standard_distribution::call(
+ mdist,
+ inf_,
+ inf_,
+ init_0,
+ init_1,
+ urng,
+ n1,
+ n2,
+ n3,
+ n4,
+ n_max_reject,
+ out
+ );
+ }
+};
}//test
}// ars
Modified: sandbox/statistics/adaptive_rejection_sampling/boost/ars/test/standard_distribution.hpp
==============================================================================
--- sandbox/statistics/adaptive_rejection_sampling/boost/ars/test/standard_distribution.hpp (original)
+++ sandbox/statistics/adaptive_rejection_sampling/boost/ars/test/standard_distribution.hpp 2010-03-03 14:03:04 EST (Wed, 03 Mar 2010)
@@ -51,103 +51,112 @@
namespace ars{
namespace test{
-// This function samples from distribution D and outputs convergence statistics
-template<typename D,typename U>
-void standard_distribution(
- const D& mdist, //e.g. D == math::normal_distribution<T>
- typename D::value_type x_min,
- typename D::value_type x_max,
- typename D::value_type init_0,
- typename D::value_type init_1,
- U& urng,
- unsigned n1, // # loops
- unsigned n2, // # subsamples per loop
- unsigned n3, // size of subsample
- unsigned n4, // At each loop, n2 *= n4
- unsigned n_max_reject,
- std::ostream& os
-){
-
- // The ars is re-initialized after each n3 sample.
- // n3 * n4 is the total size of the sample over which a KS is computed
-
- using namespace boost;
- using namespace math;
- using namespace assign;
- namespace dist = boost::statistics::detail::distribution;
- namespace ks = boost::statistics::detail::kolmogorov_smirnov;
- typedef std::string str_;
- typedef std::runtime_error err_;
- typedef typename dist::meta::value<D>::type val_;
- typedef std::vector<val_> vals_;
-
- typedef ars::function::adaptor<D> fun_t;
- typedef ars::proposal_sampler<val_,std::vector> ps_;
- typedef ars::sampler<ps_> ars_;
- typedef random::ref_distribution<ars_&> ref_ars_;
- typedef variate_generator<U&,ref_ars_> vg_ars_;
-
- typedef ks::tag::statistic<val_> tag_ks_;
- typedef boost::accumulators::stats<tag_ks_> acc_features_;
- typedef boost::accumulators::accumulator_set<val_,acc_features_> acc_;
-
- ars_ ars;
- ars.set_function(x_min, x_max, fun_t(mdist));
- {
- static const str_ str
- = "ars initialized every %1% with init_0 = %2% and init_1 = %3%";
- format f(str);
- f%n3%init_0%init_1;
- os << f.str() << std::endl;
- os << description(mdist) << std::endl;
+struct standard_distribution{
+
+ static void header(std::ostream& os){
+ os << "(a,b,c)" << std::endl;
+ os << "a : sample size" << std::endl
+ << "b : kolmogorov-smirnov statistic" << std::endl
+ << "c : number of rejections" << std::endl;
}
- long unsigned n_reject;
- os << "(sample size" << ','
- << "kolmogorov-smirnov statistic" << ','
- << "number of rejection steps)"
- << std::endl;
+
+ // Samples from distribution D using adaptive rejection sampling and outputs
+ // convergence statistics
+ template<typename D,typename U>
+ static void call(
+ const D& mdist, //e.g. D == math::normal_distribution<T>
+ typename D::value_type x_min,
+ typename D::value_type x_max,
+ typename D::value_type init_0,
+ typename D::value_type init_1,
+ U& urng,
+ unsigned n1, // # loops
+ unsigned n2, // # subsamples per loop
+ unsigned n3, // size of subsample
+ unsigned n4, // At each loop, n2 *= n4
+ unsigned n_max_reject,
+ std::ostream& os
+ ){
+
+ // The ars is re-initialized after each n3 sample.
+ // n3 * n4 is the total size of the sample over which a KS is computed
+
+ using namespace boost;
+ using namespace math;
+ using namespace assign;
+ namespace dist = boost::statistics::detail::distribution;
+ namespace ks = boost::statistics::detail::kolmogorov_smirnov;
+ typedef std::string str_;
+ typedef std::runtime_error err_;
+ typedef typename dist::meta::value<D>::type val_;
+ typedef std::vector<val_> vals_;
+
+ typedef ars::function::adaptor<D> fun_t;
+ typedef ars::proposal_sampler<val_,std::vector> ps_;
+ typedef ars::sampler<ps_> ars_;
+ typedef random::ref_distribution<ars_&> ref_ars_;
+ typedef variate_generator<U&,ref_ars_> vg_ars_;
+
+ typedef ks::tag::statistic<val_> tag_ks_;
+ typedef boost::accumulators::stats<tag_ks_> acc_features_;
+ typedef boost::accumulators::accumulator_set<val_,acc_features_> acc_;
+
+ ars_ ars;
+ ars.set_function(x_min, x_max, fun_t(mdist));
+ {
+ static const str_ str
+ = "Initialized every %1% draw(s) with x1 = %2% and x2 = %3%";
+ format f(str);
+ f%n3%init_0%init_1;
+ os << f.str() << std::endl;
+ os << description(mdist) << std::endl;
+ }
+ long unsigned n_reject;
- try{
- for(unsigned i1 = 0; i1<n1; i1++){
- n_reject = 0;
- acc_ acc;
- for(unsigned i2 = 0; i2<n2; i2++){
- try{
- ars.initialize(init_0,init_1);
- vg_ars_ vg_ars(urng,ref_ars_(ars));
- for(unsigned i = 0; i<n3; i++)
- {
- val_ x = vg_ars();
- acc(x);
- }
- // Without ref_ars_, n would be reset to 0
- n_reject += (
- (vg_ars.distribution()).distribution()
- ).n_reject();
- }catch(std::exception& e){
- format f("at i1 = %1%, i2 = %2% : %3%");
- f % i1 % i2 % e.what();
- throw std::runtime_error(f.str());
- }
- }
- val_ rate
- = static_cast<val_>(n_reject)/static_cast<val_>(n3*n2);
-
- os
- << '('
- << ks::extract::statistic<val_>(acc,mdist)
- << ','
- << rate
- << ')'
- << std::endl;
- n2 *= n4;
- }
- }catch(std::exception& e)
- {
- std::cerr << e.what() << std::endl;
- }
- os << std::endl;
-}
+ try{
+ for(unsigned i1 = 0; i1<n1; i1++){
+ n_reject = 0;
+ acc_ acc;
+ for(unsigned i2 = 0; i2<n2; i2++){
+ try{
+ ars.initialize(init_0,init_1);
+ vg_ars_ vg_ars(urng,ref_ars_(ars));
+ for(unsigned i = 0; i<n3; i++)
+ {
+ val_ x = vg_ars();
+ acc(x);
+ }
+ // Without ref_ars_, n would be reset to 0
+ n_reject += (
+ (vg_ars.distribution()).distribution()
+ ).n_reject();
+ }catch(std::exception& e){
+ format f("at i1 = %1%, i2 = %2% : %3%");
+ f % i1 % i2 % e.what();
+ throw std::runtime_error(f.str());
+ }
+ }
+ long int n = boost::accumulators::extract::count(acc);
+ val_ st = ks::extract::statistic<val_>(acc,mdist);
+ val_ rate = static_cast<val_>(n_reject)/static_cast<val_>(n3*n2);
+ os
+ << '('
+ << n
+ << ','
+ << st
+ << ','
+ << rate
+ << ')'
+ << std::endl;
+ n2 *= n4;
+ }
+ }catch(std::exception& e)
+ {
+ std::cerr << e.what() << std::endl;
+ }
+ os << std::endl;
+ }
+};
}//test
}// ars
Added: sandbox/statistics/adaptive_rejection_sampling/libs/ars/doc/convergence_test1.txt
==============================================================================
--- (empty file)
+++ sandbox/statistics/adaptive_rejection_sampling/libs/ars/doc/convergence_test1.txt 2010-03-03 14:03:04 EST (Wed, 03 Mar 2010)
@@ -0,0 +1,47 @@
+Macbook - 2.4 GHz Intel Core 2 Duo
+(a,b,c)
+a : sample size
+b : kolmogorov-smirnov statistic
+c : number of rejections
+
+Mac OS Leopard 10.6 - x86_64 - Release mode - gcc 4.2
+Initialized every 1 draw(s) with x1 = 102 and x2 = 102.01
+gamma(3,1)
+(10,0.16757,2.1)
+(100,0.0580157,2.46)
+(1000,0.0309302,2.318)
+(10000,0.00768814,2.3112)
+(100000,0.0021249,2.32135)
+
+Initialized every 1 draw(s) with x1 = 2.01 and x2 = 2.02
+gamma(3,1)
+(10,0.379978,2.3)
+(100,0.0560606,1.71)
+(1000,0.018001,1.683)
+(10000,0.00552136,1.6643)
+(100000,0.00235738,1.65932)
+
+Initialized every 1 draw(s) with x1 = -100 and x2 = 0.01
+normal(0,2)
+(10,0.153306,4.8)
+(100,0.0478689,4.29)
+(1000,0.0237086,4.679)
+(10000,0.00603333,4.6361)
+(100000,0.00306973,4.65347)
+
+Initialized every 1 draw(s) with x1 = -0.01 and x2 = 100
+normal(0,2)
+(10,0.275829,5.4)
+(100,0.0408603,4.57)
+(1000,0.0170958,4.683)
+(10000,0.00669968,4.6581)
+(100000,0.0023244,4.65068)
+
+Initialized every 1 draw(s) with x1 = -0.01 and x2 = 0.01
+normal(0,2)
+(10,0.173027,5.7)
+(100,0.104005,5.99)
+(1000,0.0356427,5.933)
+(10000,0.00754988,6.0014)
+(100000,0.00280912,5.98497)
+
Modified: sandbox/statistics/adaptive_rejection_sampling/libs/ars/doc/readme.txt
==============================================================================
--- sandbox/statistics/adaptive_rejection_sampling/libs/ars/doc/readme.txt (original)
+++ sandbox/statistics/adaptive_rejection_sampling/libs/ars/doc/readme.txt 2010-03-03 14:03:04 EST (Wed, 03 Mar 2010)
@@ -19,10 +19,6 @@
A tool for testing is provided for standard distributions, which makes use of
their known cdf to produce kolmogorov-smirnov statistics.
-[ Bugs ]
-
-Since Jan 20, convergence fails in the test.
-
[ Compiler ]
Tested Jan 20, 2010
@@ -41,10 +37,10 @@
/sandbox/statistics/random
[ History ]
-
+March 3rd, 2010 Fixed an apparent slow convergence in the output from
+ standard_distribution.cpp which was in fact due to a formatting defect.
Jan 20 2010 : small changes to test functions to adapt to the new
non_parametric::kolmogorov_smirnov functionality.
-
Jan 8 2010 : in proposal_sampler, found and fixed bug
//T t = max_element(
// begin(datas_),back_iter,
@@ -104,58 +100,3 @@
ars1d.h, started 22 April 2008 by tss
www.cs.cmu.edu/~tss/README_ars1d.txt
-[Output]
-
-Here's the output from main.cpp :
-
-[Session started at 2009-07-24 20:48:28 -0400.]
--> example_search_reflection
-N(0,2)x_0 = 100, x_1 = 100.01, , n = 14, p_0 : (-63.83,-509.284,15.9575), p_1 : (100.01,-1250.25,-25.0025)
-x_0 = -100.01, x_1 = -100, , n = 14, p_0 : (-100.01,-1250.25,25.0025), p_1 : (63.83,-509.284,-15.9575)
-x_0 = -0.02, x_1 = -0.01, , n = 2, p_0 : (-0.02,-5e-05,0.005), p_1 : (0.02,-5e-05,-0.005)
-<-
--> example_standard_distribution
-Gamma(3,1) : ars initialized every 1 with init_0 = 102 and init_1 = 102.01
-(Kolmogorov-statistic, # rejections per draw):
-((10,0.16757),2.1)
-((100,0.0580157),2.46)
-((1000,0.0309302),2.318)
-((10000,0.00768814),2.3112)
-((100000,0.0021249),2.32135)
-((1000000,0.000751711),2.32073)
-
-Gamma(3,1) : ars initialized every 1 with init_0 = 2.01 and init_1 = 2.02
-(Kolmogorov-statistic, # rejections per draw):
-((10,0.379978),2.3)
-((100,0.0560606),1.71)
-((1000,0.018001),1.683)
-((10000,0.00552136),1.6643)
-((100000,0.00235738),1.65932)
-((1000000,0.000615244),1.65827)
-
-N(0,2) : ars initialized every 1 with init_0 = -100 and init_1 = 0.01
-(Kolmogorov-statistic, # rejections per draw):
-((10,0.153306),4.8)
-((100,0.0478689),4.29)
-((1000,0.0237086),4.679)
-((10000,0.00603333),4.6361)
-((100000,0.00306973),4.65347)
-((1000000,0.000669519),4.65573)
-
-N(0,2) : ars initialized every 1 with init_0 = -0.01 and init_1 = 100
-(Kolmogorov-statistic, # rejections per draw):
-((10,0.275829),5.4)
-((100,0.0408603),4.57)
-((1000,0.0170958),4.683)
-((10000,0.00669968),4.6581)
-((100000,0.0023244),4.65068)
-((1000000,0.000614343),4.65194)
-
-N(0,2) : ars initialized every 1 with init_0 = -0.01 and init_1 = 0.01
-(Kolmogorov-statistic, # rejections per draw):
-((10,0.173027),5.7)
-((100,0.104005),5.99)
-((1000,0.0356427),5.933)
-((10000,0.00754988),6.0014)
-((100000,0.00280912),5.98497)
-((1000000,0.0013802),5.99028)
Modified: sandbox/statistics/adaptive_rejection_sampling/libs/ars/example/standard_distribution.cpp
==============================================================================
--- sandbox/statistics/adaptive_rejection_sampling/libs/ars/example/standard_distribution.cpp (original)
+++ sandbox/statistics/adaptive_rejection_sampling/libs/ars/example/standard_distribution.cpp 2010-03-03 14:03:04 EST (Wed, 03 Mar 2010)
@@ -26,19 +26,22 @@
typedef double value_;
- const unsigned n1 = 1e1; // # loops
- const unsigned n2 = 1e1; // # subsamples on first loop
- const unsigned n3 = 1e0; // size of subsample
- const unsigned n4 = 2e0; // At each loop, n2 *= n4
+ const unsigned n1 = 5e0; // # loops
+ const unsigned n2 = 1e1; // # subsamples on first loop
+ const unsigned n3 = 1e0; // size of subsample
+ const unsigned n4 = 1e1; // At each loop, n2 *= n4
const unsigned max_n_reject = 10;
+ ars::test::standard_distribution::header(out);
+
// The initial values are chosen to test the robustness of the
// implementation, within the range allowed by the algorithm.
{ // Domain = [0,inf)
value_ shape = 3.0;
value_ scale = 1.0;
value_ mode = (shape - 1.0) * scale; //2
- ars::test::gamma_distribution<value_>(
+ typedef ars::test::gamma_distribution test_;
+ test_::call(
shape,
scale,
mode + 100.0, //init0
@@ -50,7 +53,7 @@
max_n_reject,
out
);
- ars::test::gamma_distribution<value_>(
+ test_::call(
shape, //shape
scale, //scale
mode + 0.01,//init0
@@ -67,7 +70,8 @@
{ // Domain = (-inf,inf)
value_ mu = 0.0;
value_ sigma = 2.0;
- ars::test::normal_distribution<value_>(
+ typedef ars::test::normal_distribution test_;
+ test_::call(
mu,
sigma,
mu -100.0, //-100.0, //init0
@@ -80,7 +84,7 @@
out
);
- ars::test::normal_distribution<value_>(
+ test_::call(
mu,
sigma,
mu - 0.01, //init0
@@ -93,7 +97,7 @@
out
);
- ars::test::normal_distribution<value_>(
+ test_::call(
mu,
sigma,
-0.01, //init0
Boost-Commit list run by bdawes at acm.org, david.abrahams at rcn.com, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk