|
|
Boost-Commit : |
Subject: [Boost-commit] svn:boost r52001 - in sandbox/conditionally_specified_distribution/adaptive_rejection_sampling: boost/adaptive_rejection_sampling boost/adaptive_rejection_sampling/example libs/adaptive_rejection_sampling/src libs/adaptive_rejection_sampling/src/doc libs/adaptive_rejection_sampling/src/example
From: erwann.rogard_at_[hidden]
Date: 2009-03-26 18:18:08
Author: e_r
Date: 2009-03-26 18:18:07 EDT (Thu, 26 Mar 2009)
New Revision: 52001
URL: http://svn.boost.org/trac/boost/changeset/52001
Log:
Fixed bugs and more execption handling
Text files modified:
sandbox/conditionally_specified_distribution/adaptive_rejection_sampling/boost/adaptive_rejection_sampling/approximation.hpp | 10
sandbox/conditionally_specified_distribution/adaptive_rejection_sampling/boost/adaptive_rejection_sampling/approximation_impl.hpp | 904 +++++++++++++++++++++++----------------
sandbox/conditionally_specified_distribution/adaptive_rejection_sampling/boost/adaptive_rejection_sampling/example/standard_gaussian.hpp | 7
sandbox/conditionally_specified_distribution/adaptive_rejection_sampling/boost/adaptive_rejection_sampling/point.hpp | 74 +-
sandbox/conditionally_specified_distribution/adaptive_rejection_sampling/boost/adaptive_rejection_sampling/sampler.hpp | 13
sandbox/conditionally_specified_distribution/adaptive_rejection_sampling/libs/adaptive_rejection_sampling/src/doc/readme.txt | 6
sandbox/conditionally_specified_distribution/adaptive_rejection_sampling/libs/adaptive_rejection_sampling/src/example/test_data.cpp | 14
sandbox/conditionally_specified_distribution/adaptive_rejection_sampling/libs/adaptive_rejection_sampling/src/example/test_data.h | 4
sandbox/conditionally_specified_distribution/adaptive_rejection_sampling/libs/adaptive_rejection_sampling/src/example/test_gaussian.cpp | 56 +-
sandbox/conditionally_specified_distribution/adaptive_rejection_sampling/libs/adaptive_rejection_sampling/src/main.cpp | 47 +
10 files changed, 673 insertions(+), 462 deletions(-)
Modified: sandbox/conditionally_specified_distribution/adaptive_rejection_sampling/boost/adaptive_rejection_sampling/approximation.hpp
==============================================================================
--- sandbox/conditionally_specified_distribution/adaptive_rejection_sampling/boost/adaptive_rejection_sampling/approximation.hpp (original)
+++ sandbox/conditionally_specified_distribution/adaptive_rejection_sampling/boost/adaptive_rejection_sampling/approximation.hpp 2009-03-26 18:18:07 EDT (Thu, 26 Mar 2009)
@@ -36,6 +36,10 @@
// Modifier:
// initialize(x1,x2)
// update(x)
+ //
+ // In case this is used to sample from {Y~(Y|x[i]):i=1,...,n},
+ // consider using reset_distribution_function(const dist_f_t& f)
+ // rather than creating a new object to save allocation
template<
typename DistFun,
template<typename,typename> class Cont = std::vector,
@@ -77,7 +81,6 @@
max_recursion_(max_recursion){
this->initialize(x0_init,x1_init);
}
-
approximation(const approximation& that)
:super_t(that),
dist_f_(that.dist_f_),
@@ -93,6 +96,11 @@
return *this;
}
public:
+ void reset_distribution_function(const dist_f_t& f){
+ dist_f_ = f;
+ //next : initialize
+ }
+
void initialize(
typename super_t::value_type x1,typename super_t::value_type x2){
try{
Modified: sandbox/conditionally_specified_distribution/adaptive_rejection_sampling/boost/adaptive_rejection_sampling/approximation_impl.hpp
==============================================================================
--- sandbox/conditionally_specified_distribution/adaptive_rejection_sampling/boost/adaptive_rejection_sampling/approximation_impl.hpp (original)
+++ sandbox/conditionally_specified_distribution/adaptive_rejection_sampling/boost/adaptive_rejection_sampling/approximation_impl.hpp 2009-03-26 18:18:07 EDT (Thu, 26 Mar 2009)
@@ -7,20 +7,27 @@
#ifndef BOOST_ADAPTIVE_REJECTION_SAMPLING_APPROXIMATION_IMPL_HPP_ER_2009
#define BOOST_ADAPTIVE_REJECTION_SAMPLING_APPROXIMATION_IMPL_HPP_ER_2009
#include <vector>
+#include <list>
#include <ext/algorithm> // is_sorted
#include <algorithm>
#include <functional>
#include <iterator>
#include <cmath>
#include <boost/mpl/identity.hpp>
+#include <boost/mpl/bool.hpp>
#include <boost/assert.hpp>
+#include <boost/static_assert.hpp>
#include <boost/bind.hpp>
#include <boost/function.hpp>
#include <boost/range.hpp>
+#include <boost/math/special_functions/fpclassify.hpp>
#include <boost/math/tools/precision.hpp>
+#include <boost/math/special_functions/next.hpp>
#include <boost/math_limit/miscellanea.hpp>
#include <boost/utility/container_to_string.hpp>
#include <boost/adaptive_rejection_sampling/point.hpp>
+//TODO : make use of
+//<boost/math/special_functions/fpclassify.hpp>
namespace boost{
namespace adaptive_rejection_sampling{
@@ -49,265 +56,328 @@
// Modifier:
// void initialize_impl(const point_t& a,const point_t& b)
// update_impl(const point_t& p)
- template<
- typename RealType = double,
- template<typename,typename> class Cont = std::vector,
- template<typename> class Alloc = std::allocator
- >
- class approximation_impl{
- public:
- typedef RealType value_type;
- typedef math_limit::limits<value_type> limits_t;
- typedef point<value_type> point_t;
- typedef Alloc<point_t> alloc_t;
- typedef Cont<point_t,alloc_t> points_t;
- typedef std::vector<value_type> upper_t;
- typedef typename range_iterator<points_t>::type iter_t;
- typedef typename range_iterator<upper_t>::type upper_iter_t;
- typedef typename range_iterator<const points_t>::type
+ //
+ template<
+ typename RealType = double,
+ template<typename,typename> class Cont = std::list,
+ template<typename> class Alloc = std::allocator
+ >
+ class approximation_impl{
+ public:
+ typedef RealType value_type;
+ typedef point<value_type> point_t;
+ typedef Alloc<point_t> alloc_t;
+ typedef Cont<point_t,alloc_t> points_t;
+ typedef std::vector<value_type> upper_t;
+ typedef typename range_iterator<points_t>::type iter_t;
+ typedef typename range_iterator<upper_t>::type upper_iter_t;
+ typedef typename range_iterator<const points_t>::type
citer_t;
- typedef typename range_iterator<const upper_t>::type
+ typedef typename range_iterator<const upper_t>::type
cupper_iter_t;
- typedef typename range_size<points_t>::type size_type;
- struct result_of_range_points : mpl::identity<const points_t&>{};
- struct result_of_range_upper : mpl::identity<const upper_t&>{};
-
- bool is_valid_leftmost_point(value_type dy){
- return (!limits_t::is_negative_infinity(min())) || (dy>0.0);
- }
- bool is_valid_rightmost_point(value_type dy){
- return (!limits_t::is_infinity(max())) || (dy<0.0);
- }
-
- approximation_impl(
- value_type min,
- value_type max
- ):min_(min),max_(max),subtract_from_y_(0.0)
- {
- //Warning not initialized.
- }
-
- approximation_impl(
- value_type min,
- value_type max,
- const point_t& a,
- const point_t& b
- ):min_(min),max_(max_),subtract_from_y_(0.0)
- {
- initialize_impl(a,b);
- }
- //TODO protected
- void initialize_impl(
- const point_t& a,
- const point_t& b
- ){
- BOOST_ASSERT(a.x()<b.x());
- BOOST_ASSERT(a.dy()>b.dy());
- BOOST_ASSERT(is_valid_leftmost_point(a.dy()));
- BOOST_ASSERT(is_valid_rightmost_point(b.dy()));
- points_.clear();
- points_.push_back(a);
- points_.push_back(b);
- upper_.clear();
- upper_.push_back(min());
- value_type ab = tangent_intersection(a,b);
- upper_.push_back(ab);
- upper_.push_back(max());
- assert_upper_points();
- assert_points_sorted();
- assert_range_dy_sorted();
- assert_upper_sorted();
- update_unnormalized_cdf();
-
- }
- approximation_impl(const approximation_impl& that)
- :min_(that.min_),
- max_(that.max_),
- points_(that.points_),
- upper_(that.upper_),
- range_unnormalized_cdf_(
- that.range_unnormalized_cdf_){}
-
- approximation_impl&
- operator=(const approximation_impl& that){
- if(&that!=this){
- min_ = that.min_;
- max_ = that.max_;
- points_ = that.points_;
- upper_ = that.upper_;
- range_unnormalized_cdf_
+ typedef typename range_size<points_t>::type size_type;
+ struct result_of_range_points : mpl::identity<const points_t&>{};
+ struct result_of_range_upper : mpl::identity<const upper_t&>{};
+
+ bool is_valid_leftmost_point(value_type dy){
+ return ( (min()>lowest()) || (dy>0.0) );
+ }
+ bool is_valid_rightmost_point(value_type dy){
+ return ( (max()<highest()) || (dy<0.0) );
+ }
+
+ approximation_impl(
+ value_type min,
+ value_type max
+ ):min_(min),max_(max)
+ {
+ //Warning not initialized.
+ }
+
+ approximation_impl(
+ value_type min,
+ value_type max,
+ const point_t& a,
+ const point_t& b
+ ):min_(min),max_(max_)
+ {
+ initialize_impl(a,b);
+ }
+
+ // Optional, before initialize
+ void reserve(size_type n){
+ //Cont == vector only
+ points_.reserve(n);
+ upper_.reserve(n);
+ range_unnormalized_cdf_.reserve(n);
+ }
+
+ size_type size_range_point()const{
+ return size(points_);
+ }
+
+ //TODO protected
+ void initialize_impl(
+ const point_t& a,
+ const point_t& b
+ ){
+ BOOST_ASSERT(a.x()<b.x());
+ BOOST_ASSERT(a.dy()>b.dy());
+ BOOST_ASSERT(is_valid_leftmost_point(a.dy()));
+ BOOST_ASSERT(is_valid_rightmost_point(b.dy()));
+ points_.clear();
+ points_.push_back(a);
+ points_.push_back(b);
+ upper_.clear();
+ upper_.push_back(min());
+ value_type ab = tangent_intersection(a,b);
+ upper_.push_back(ab);
+ upper_.push_back(max());
+ assert_upper_points();
+ assert_points_sorted();
+ assert_range_dy_sorted();
+ assert_upper_sorted();
+ update_unnormalized_cdf();
+
+ }
+ approximation_impl(const approximation_impl& that)
+ :min_(that.min_),
+ max_(that.max_),
+ points_(that.points_),
+ upper_(that.upper_),
+ range_unnormalized_cdf_(that.range_unnormalized_cdf_),
+ shift_(that.shift_){}
+
+ approximation_impl&
+ operator=(const approximation_impl& that){
+ if(&that!=this){
+ min_ = that.min_;
+ max_ = that.max_;
+ points_ = that.points_;
+ upper_ = that.upper_;
+ range_unnormalized_cdf_
= that.range_unnormalized_cdf_;
+ shift_ = that.shift_;
//unnormalized_cdf_ = that.unnormalized_cdf_;
- }
- return *this;
- }
-
- value_type min()const{return min_;}
- value_type max()const{return max_;}
- value_type
- total_unnormalized_cdf()const{
- BOOST_ASSERT(size(points_)>1);
- return range_unnormalized_cdf_.back();
- //warning : this is an approximation, whose magnitude
- // is compounded by subtract_y (Ctrl+F), but as n-->inf,
- // sutract_y should --->0
-
- }
- value_type
- lower(value_type x)const{
- BOOST_ASSERT(size(points_)>1);
- value_type res = 0.0;
- const point_t p(x,0.0,0.0);
- citer_t i = upper_bound(
- begin(points_),
- end(points_),
- p,
- less<value_type>()
- );
- if(i!=end(points_)){
- if(begin(points_)!=i){
- value_type x_h = i -> x();
- value_type y_h = i -> y();
- --i;
- value_type x_l = i -> x();
- value_type y_l = i -> y();
- res = y_l + (x-x_l) * (y_h-y_l) / (x_h-x_l);
- }else{
- res = negative_infinity();
- }
- }else{
- res = negative_infinity();
- }
- return res;
- }
- value_type
- upper(value_type x)const{
- BOOST_ASSERT(size(points_)>1);
+ }
+ return *this;
+ }
+
+ value_type min()const{return min_;}
+ value_type max()const{return max_;}
+ value_type
+ total_unnormalized_cdf()const{
+ BOOST_ASSERT(size(points_)>1);
+ return range_unnormalized_cdf_.back();
+ // Warning : This is a shifted normalizing constant i.e. in
+ // general it will not approach the true normalizing const
+ }
+ value_type
+ lower(value_type x)const{
+ BOOST_ASSERT(size(points_)>1);
+ try{
+ const point_t p(x,0.0,0.0);
+ citer_t i = upper_bound(
+ begin(points_),
+ end(points_),
+ p,
+ less<value_type>()
+ );
+ if(i!=end(points_)){
+ if(begin(points_)!=i){
+ citer_t prior = i; --prior;
+ return linear_approximation(*prior,*i,x);
+ }else{
+ //TODO That, or lower_bound?
+ return lowest();
+ }
+ }else{
+ return lowest();
+ }
+ }catch(std::exception& e){
+ std::string str = "lower(%1%)";
+ format f(str); f%x;
+ throw std::runtime_error(f.str());
+ }
+ }
+ value_type
+ upper(value_type x)const{
+ BOOST_ASSERT(size(points_)>1);
// Loops
// a b c
// / \ / \ / \
// z ab bc z (tangent intersections)
// x is_in [ab,bc], we need tangent at b
- cupper_iter_t u = upper_bound(
- begin(upper_),
- end(upper_),
- x
- );
- typedef typename range_difference<upper_t>::type u_diff_t;
-
- u_diff_t u_diff = std::distance(begin(upper_),u);
- BOOST_ASSERT(u_diff>0);
- citer_t i = begin(points_);
- std::advance(i,u_diff-1);
- BOOST_ASSERT((i->x())<*u);
- --u; BOOST_ASSERT(*u<=x);
- return tangent(*i,x);
- }
- value_type inverse_cdf(value_type u)const{
- BOOST_ASSERT(size(points_)>1);
-
- // a---b---c---d
- // /\ / \ / \ / \
- // z- ab bc cd z+
- // | | | |
-
- const upper_t& rcdf
- = range_unnormalized_cdf_;
- cupper_iter_t i_cdf = begin(rcdf);
- cupper_iter_t i_u = begin(upper_);
- citer_t i_p = begin(points_);
- value_type thresh = u * total_unnormalized_cdf();
- while(*i_cdf< thresh){
- ++i_cdf;
- ++i_u;
- ++i_p;
- BOOST_ASSERT(
+ cupper_iter_t u = upper_bound(
+ begin(upper_),
+ end(upper_),
+ x
+ );
+ typedef typename range_difference<upper_t>::type u_diff_t;
+
+ u_diff_t u_diff = std::distance(begin(upper_),u);
+ BOOST_ASSERT(u_diff>0);
+ citer_t i = begin(points_);
+ std::advance(i,u_diff-1);
+ BOOST_ASSERT((i->x())<*u);
+ --u; BOOST_ASSERT(*u<=x);
+ return tangent_impl(*i,x);
+ }
+ value_type inverse_cdf(value_type u)const{
+ BOOST_ASSERT(size(points_)>1);
+// x x0 x1 x
+// /\ / \ / \ / \
+// z- z z0 z1 z+
+// | | | |
+// c0 c1 c
+// ^
+// c*u
+// where c1 = (exp(t1-shift)-exp(t0-shift)) / dy1 > 0.
+// Let a = exp(t0 - shift) + dy1 (c u - c0)
+// The quantile is
+// q = x1 + ( log ( a ) - (y1-shift) ) / dy1
+// We have: (c u - c0) < c1 hence |dy1 (cu - c0)| < |dy1| c1
+// if dy1>0, 0< exp(t0-shift) < a < exp(t1-shift)
+// if dy1<0, 0< exp(t1-shift) < a < exp(t0-shift)
+// Hence we want shift such that
+// exp(t1-shift) < inf
+// exp(t0-shift) < inf
+//
+// From the above, we have:
+// dy1>0 : x1 + (1/dy) * (t0-y1) < q < x1 + (1/dy) * (t1-y1)
+// dy1>0 : x1 + (1/dy) * (t0-y1) < q < x1 + (1/dy) * (t1-y1)
+
+ value_type thresh, x, y, dy, t, cumul_cdf, a, res;
+ try{
+
+ const upper_t& rcdf = range_unnormalized_cdf_;
+ cupper_iter_t i_cdf = begin(rcdf);
+ cupper_iter_t i_u = begin(upper_);
+ citer_t i_p = begin(points_);
+ thresh = u * total_unnormalized_cdf();
+ while(*i_cdf< thresh){
+ ++i_cdf;
+ ++i_u;
+ ++i_p;
+ BOOST_ASSERT(
i_cdf!=end(range_unnormalized_cdf_)
- );
- }
- value_type x = i_p->x();
- value_type y = i_p->y() - subtract_from_y_;
- value_type dy = i_p->dy();
- //TODO distinguish based on is_negative_infinity(*i_u)
- value_type t = tangent(*i_p,*i_u) - subtract_from_y_;
- value_type cumul_cdf = 0.0;
- if(i_cdf!=begin(rcdf)){
- --i_cdf;
- cumul_cdf = *i_cdf;
- }
- value_type res
- = math_limit::safer_exp(t)+ dy * (thresh-cumul_cdf);
- if(!(res>0.0)){
- std::string str = "adaptive_rejection_sampling::";
- str+= "approximation_impl::";
- str+= "inverse_cdf : ";
- str+= "!res=%1%>0.0";
- str+= "x=%2% subtract_from_y_=%3% t=%4%";
- format f(str);
- f%res%x%subtract_from_y_;
- f%t;
+ );
+ }
+ x = i_p->x();
+ y = i_p->y();
+ dy = i_p->dy();
+ t = tangent_impl(*i_p,*i_u);
+
+ cumul_cdf = 0.0;
+ if(i_cdf!=begin(rcdf)){
+ --i_cdf;
+ cumul_cdf = *i_cdf;
+ }
+
+ //value_type e = exp_shift(t-log(fabs(dy)));
+ //a = e;
+ //if(dy<0.0){
+ // a -= (thresh-cumul_cdf);
+ //}else{
+ // a += (thresh-cumul_cdf);
+ //}
+ //a *= fabs(dy);
+
+ value_type e = exp_shift(t);
+ a = e + dy * (thresh-cumul_cdf);
+
+ if(!(a>0.0)){
+ format f("!a=%1%>0.0 \n"); f%a;
+ throw std::runtime_error(f.str());
+ }else{
+ if(!(a<highest())){
+ format f("!(a=%1%<highest()=%2%)");
+ f%1%highest();
throw std::runtime_error(f.str());
- }else{ //safer_log instead?
- if(!(res<limits_t::infinity())){
- std::string str = "inverse_cdf";
- str+= "!(res=%1%<limits_t::max_value()=%2%)";
- str+= "x=%3%";
- format f(str);
- f%res%limits_t::infinity()%x;
+ }else{
+ res = x + ( log(a) - shift(y) ) / dy;
+
+ if((math::isnan)(res)){
+ std::string str = "!isnan( x + ( log(a) - y ) / dy)";
+ str+= "x = %1%, a = %2%, y = %3%, dy = %4%";
+ format f(str); f%x%a%y%dy;
throw std::runtime_error(f.str());
- }else{
- res = x + ( log(res) - y ) / dy;
}
}
- return res;
- }
-
- // For testing only.
- typename result_of_range_points::type
- range_point()const{ return points_; }
- typename result_of_range_upper::type
- range_upper()const{ return upper_; }
- //TODO protected
- void update_impl(value_type x,value_type y,value_type dy){
- return update_impl(point_t(x,y,dy));
}
- //TODO protected
- void update_impl(const point_t& p){
- BOOST_ASSERT(size(points_)>1);
-
+ }catch(std::exception& e){
+ std::string str = "adaptive_rejection_sampling::";
+ str+= "approximation_impl::";
+ str+= "inverse_cdf(%1%) : ";
+ str+= "thresh = %2%, x = %3%, y = %4%, dy = %5%, t = %6%";
+ str+= "cumul_cdf = %7%, res = %8%";
+ format f(str); f%u%thresh%x%y%dy%t%cumul_cdf%res;
+ str+= e.what();
+ str+= as_string();
+ }
+ return res;
+ }
+
+ // For testing only.
+ typename result_of_range_points::type
+ range_point()const{ return points_; }
+ typename result_of_range_upper::type
+ range_upper()const{ return upper_; }
+ //TODO protected
+ void update_impl(value_type x,value_type y,value_type dy){
+ return update_impl(point_t(x,y,dy));
+ }
+ //TODO protected
+ void update_impl(const point_t& p){
+ BOOST_ASSERT(size(points_)>1);
+ try{
// lower_bound returns furthermost iterator i in [first, last) such that,
// for every j in [first, i), *j < value or (comp(*j, value) is true)
// upper_bound returns the furthermost iterator i in [first, last) such that
//, for every iterator j in [first, i), value < *j is false.
// or comp(value, *j) is false.
- iter_t i = upper_bound(
+ iter_t i = upper_bound(
begin(points_),
end(points_),
p,
less<value_type>()
- );
- bool is_not_largest = (i!=end(points_));
- i = points_.insert(i,p);
- bool is_not_smallest = (i!=begin(points_));
- value_type z_l,z_r = 0.0; //left
-
- if(is_not_smallest){
- iter_t l = i; --l;
- z_l = tangent_intersection(*l,*i);
- }
- if(is_not_largest){
+ );
+ bool is_not_largest = (i!=end(points_));
+ i = points_.insert(i,p);
+ bool is_not_lowest = (i!=begin(points_));
+ value_type z_l,z_r = 0.0;
+
+ if(is_not_lowest){
+ iter_t l = i; --l;
+ z_l = tangent_intersection(*l,*i);
+ if((math::isnan)(z_l)){
+ std::string str = "(math::isnan)(z_l=%1%) with ";
+ str+="l : "; str+=l->as_string();
+ str+="i : "; str+=i->as_string();
+ format f(str); f%z_l;
+ throw std::runtime_error(f.str());
+ }
+ }
+ if(is_not_largest){
iter_t r = i; ++r;
z_r = tangent_intersection(*i,*r);
- }
- // z[0] = m, z[1],...., z[n-1]=M, end
- if(!is_not_smallest){
+ if((math::isnan)(z_l)){
+ std::string str = "(math::isnan)(z_r=%1%) with ";
+ str+="i : "; str+=i->as_string();
+ str+="r : "; str+=r->as_string();
+ format f(str); f%z_r;
+ throw std::runtime_error(f.str());
+ }
+ }
+ // z[0] = m, z[1],...., z[n-1]=M, (end)
+ if(!is_not_lowest){
upper_iter_t j = begin(upper_); ++j; //*j = z[1]
upper_.insert(j,z_r);
- }else{
- if(is_not_largest){
+ }else{
+ if(is_not_largest){
// a b c
// / \ / \
// z- ac z+
@@ -321,149 +391,259 @@
// a b c
// / \ / \ / \
// z- ab bc z+
- }else{
+ }else{
upper_iter_t j = begin(upper_);
std::advance(j,size(upper_)-1); //*j = z+
upper_.insert(j,z_l);
- }
- }
-
- assert_upper_points();
- assert_points_sorted();
- assert_range_dy_sorted();
- assert_upper_sorted();
- update_unnormalized_cdf();
}
+ }
+ assert_upper_points();
+ assert_points_sorted();
+ assert_range_dy_sorted();
+ assert_upper_sorted();
+ update_unnormalized_cdf();
+ }catch(std::exception& e){
+ std::string str = "adaptive_rejection_sampling::";
+ str+= "approximation_impl::";
+ str+= "update_impl : ";
+ str+= e.what();
+ str+= " ";
+ str+= this->as_string();
+ throw std::runtime_error(str);
+ }
- std::string as_string()const{
+ }
+
+ std::string as_string()const{
typedef utility::container_to_string to_str_t;
std::string str;
str+= "[ upper_ = %1%\n";
str+= " range_unnormalized_cdf_ = %2%\n";
str+= " {x} = %3%\n";
str+= " {y} = %4%\n";
- str+= " {dy} = %5% ]";
+ str+= " {dy} = %5%";
+ str+= " shift = %6%]";
format f(str);
f%to_str_t::get_indexed(upper_);
f%to_str_t::get_indexed(range_unnormalized_cdf_);
f%to_str_t::get_indexed(range_x());
f%to_str_t::get_indexed(range_y());
f%to_str_t::get_indexed(range_dy());
+ f%shift_;
return f.str();
- }
- protected:
- value_type
- negative_infinity()const{
- return limits_t::negative_infinity();
- }
- value_type infinity()const{
- return limits_t::infinity();
- }
-
- private:
-
- void update_unnormalized_cdf(){
- upper_t upper_ys;
- value_type max_t = negative_infinity();
- //min_nni_t = infinity(); // no use, for now.
- subtract_from_y_ = 0.0;
- {
- citer_t i = begin(points_);
- citer_t e = end(points_);
- cupper_iter_t u = begin(upper_);
- value_type t = tangent(*i,*u);
- upper_ys.push_back(t);
- ++u;
- while(i < e){
- BOOST_ASSERT(u!=end(upper_));
- t = tangent(*i,*u);
- upper_ys.push_back(t);
- max_t = (t<max_t)?max_t : t;
- //min_nni_t
- // = ( (t<min_nni_t) && (!is_negative_infinity(t)) )
- // ?t : min_nn_t;
- ++i, ++u;
- }
- BOOST_ASSERT(max_t > negative_infinity());
- }
+ }
+ protected:
+ value_type log_max_value()const{
+ typedef math_limit::exp_support<value_type> exp_sup_t;
+ return exp_sup_t::highest();
+ }
+ value_type log_min_value()const{
+ typedef math_limit::exp_support<value_type> exp_sup_t;
+ return exp_sup_t::lowest();
+ }
+ value_type highest()const{
+ return math::tools::max_value<value_type>();
+ }
+ value_type lowest()const{
+ return (-highest());
+ }
+
+ private:
+
+ value_type tangent_impl(const point_t& p0,value_type x)const{
+ value_type t = tangent(p0,x);
+ if(math::isnan(t)){
+ std::string str = "at p0 : %1% isnan(tangent(p0,%2%))";
+ format f(str); f%p0.as_string()%x;
+ std::runtime_error(f.str());
+ }
+ return t;
+ }
+ // Warning: use only in conjunction with unnormalized cdf
+ // i.e. not for lower(...), upper(..)
+ value_type shift(value_type y)const{
+ return (y - shift_);
+ }
+
+ // If need arises later : tangent_shift
+ // return (p.y()-shift_)+((x-p.x()) * p.dy()-shift_);
+
+ value_type exp_shift(value_type y)const{
+ value_type s = shift(y);
+ value_type e = exp(s);
+ if(math::isinf(e)){
+ std::string str = "isinf(exp_shift(y)), y = %1%, ";
+ str += "shift(y) = %2%.";
+ format f(str); f%y%s;
+ throw std::runtime_error(f.str());
+ }else{
+ return e;
+ }
+ }
+
+ void update_unnormalized_cdf(){
+ // Since we work with the unnormalized cdf we are free
+ // to shift it by an arbitrary amount:
+ // s0 := t0 - shift
+ // s1 := t1 - shift
+ // delta = (exp(s1)-exp(s0))/dy1
+ // Property of log-concave distribution :
+ // (exp(s1)-exp(s0))/dy1 > 0.0
+ //
+ // m = min{shift :
+ // exp(s0) < inf (see inverse_cdf)
+ // exp(s1) < inf (see inverse_cdf)
+ // delta < inf
+ // sum{delta[i]:i=1,...,n} < inf
+ // }
+ //
+ // M = max{shift: exp(max(s0,s1)) > 0.0 }
+ //
+ // We may not always have m<M for all i, so shift>m takes precedence
- { // narrows support to [log_min_value,log_max_value]
- value_type max_allowed = limits_t::log_max_value();
- subtract_from_y_ = max_t - max_allowed;
- subtract_from_y_
- = (subtract_from_y_>0.0)? subtract_from_y_ : 0.0;
- upper_iter_t b = begin(upper_ys);
- upper_iter_t e = end(upper_ys);
- if(subtract_from_y_>0.0){
- for(upper_iter_t i = b; i<e; i++){
- if(
- *i>subtract_from_y_+limits_t::log_min_value()
- ){
- (*i) -= subtract_from_y_;
- }
- }
- }
+ bool do_restart = true;
+ value_type t0, t1, e0, e1, delta, l1, w0, w1, dy;
+ value_type min_shift = lowest();
+ shift_ = 0.0;
+ value_type max_shift = highest();
+ value_type unnormalized_cdf;
+ size_type count = 0;
+ // (1) incompatible with (2)-(4) (TODO verify), hence at most 3
+ size_type max_count = size(points_) * 3;
+ try{
+ while(do_restart){
+ if(++count > max_count){
+ std::string str = " exceed max_count = %1%";
+ format f(str); f%(value_type)(max_count);
+ throw std::runtime_error(str);
}
range_unnormalized_cdf_.clear();
- value_type unnormalized_cdf = 0.0;
citer_t i = begin(points_);
citer_t e = end(points_);
- cupper_iter_t i_t = begin(upper_ys);
- value_type t_prior = *i_t; //z-
- value_type et_prior = math_limit::safer_exp(t_prior);
- ++i_t;
- value_type t, et, dy, delta = 0.0;
+ cupper_iter_t u0 = begin(upper_);
+ cupper_iter_t u1 = u0; ++u1;
+ unnormalized_cdf = 0.0;
while(i < e){
- BOOST_ASSERT(i_t<end(upper_ys));
- t = *i_t;
- et = math_limit::safer_exp(t);
+
+ // t0[i] = t1[i-1], but probably safer this way :
+ t0 = tangent_impl(*i,*u0);
+ t1 = tangent_impl(*i,*u1);
+ // if i==begin and *u0 = min(),
+ // t0 < 0 (must) and -inf allowed
+ // if ++i==end and *u1 = max(),
+ // t1 < 0 (must) and -inf allowed
+
dy = i->dy();
- // sometimes et-et_prior = 0 even if (t-t_prior) !=0
- // due to rounding
- delta = (t - t_prior);
+ delta = (t1 - t0);
if(
- ((dy>0.0) && (delta >0.0))
- || ((dy<0.0) && (delta<0.0))
+ ( (delta>0.0) && (dy<0.0) )
+ ||( (delta<0.0) && (dy>0.0) )
)
{
- //Warning. Do not un-comment as leads to failure
- //if(delta<limits_t::min_value() * dy){
- // delta = limits_t::min_value();
- //}else{
- delta = (et-et_prior)/dy;
- //}
- }else{
-
- //math_limit::safer_exp(t); //debugging only
- typedef utility::container_to_string to_str_t;
- std::string str = "adaptive_rejection_sampling::";
- str+= "approximation_impl";
- str+= "!(delta=%1%>0.0)\n";
- str+= "subtract_from_y_=%2%\n, t_prior = %3%";
- str+= " t = %4%, dy = %5%, et_prior = %6%";
- str+= " et = %7%\n";
- str+= " upper_ys = %8%\n";
- str+= this->as_string();
- format f(str);
- f%delta%subtract_from_y_;
- f%t_prior%t%dy%et_prior%et;
- f%to_str_t::get_indexed(upper_ys);
+ std::string str = "!((t1-t0)*dy=%1%>0.0)\n";
+ format f(str); f%delta;
throw std::runtime_error(f.str());
+ }else{
+ // Allow possibility that delta == 0
}
- t_prior = t;
- et_prior = et;
- unnormalized_cdf += delta;
- range_unnormalized_cdf_.push_back(unnormalized_cdf);
- ++i; ++i_t;
- } //while
- assert_range_unnormalized_cdf_points();
- assert_range_unnormalized_cdf_sorted();
- }
-
+ //l1 = fabs(dy);
+ //l1 = log(l1);
+ //w0 = t0 - l1;
+ //w1 = t1 - l1;
+ {
+ value_type tmp = ( (t1<t0)?t0:t1 ) - log_min_value();
+ max_shift = (tmp<max_shift)? tmp : max_shift;
+ }
+ if((max_shift > min_shift) && (shift_>max_shift)){
+ shift_ = boost::math::float_prior(max_shift);
+ do_restart = true; // (1)
+ }else{
+ {
+ value_type tmp = ((t1>t0)? t1 : t0 );
+ tmp -= log_max_value();
+ min_shift = (tmp>min_shift)? tmp : min_shift;
+ }
+ if(shift_<min_shift){
+ shift_ = boost::math::float_next(min_shift);
+ do_restart = true; // (2)
+ }else{
+// e0 = exp_shift(w0);
+// e1 = exp_shift(w1);
+// delta = e1 - e0;
+// if(dy<0.0){ delta = (-delta); }
+ e0 = exp_shift(t0);
+ e1 = exp_shift(t1);
+ delta = (e1 - e0);
+
+ value_type c = math_limit::scale_to_max::divides(
+ fabs(delta),
+ fabs(dy)
+ );
+ if(c>1.0){
+ min_shift = shift_ + log(c);
+ shift_ = boost::math::float_next(min_shift);
+ do_restart = true; // (3)
+ }else{
+ delta /= dy;
+ if(math::isnan(delta)){
+ std::string str;
+ str = "math::isnan(delta)";
+ throw std::runtime_error(str);
+ }
+ value_type c
+ = math_limit::scale_to_max::plus(
+ unnormalized_cdf,
+ delta
+ );
+ if(c>1.0){
+ min_shift = shift_ + log(c);
+ shift_
+ = boost::math::float_next(min_shift);
+ do_restart = true; // (4)
+ }else{
+ unnormalized_cdf+=delta;
+ if(math::isinf(unnormalized_cdf)){
+ std::string str;
+ str = "math::isinf(unnormalized_cdf)";
+ throw std::runtime_error(str);
+ }
+ range_unnormalized_cdf_.push_back(
+ unnormalized_cdf);
+ ++i; ++u0; ++u1;
+ do_restart = false;
+ }
+ }
+ }
+ }
+ if(do_restart){
+ i = e; //ensures exit
+ }
+ }//while
+ } //while
+ assert_range_unnormalized_cdf_points();
+ assert_range_unnormalized_cdf_sorted();
+ if(!(range_unnormalized_cdf_.back()>0.0)){
+ std::string str = "!(range_unnormalized_cdf_.back()>0.0)";
+ throw std::runtime_error(str);
+ }
+ }catch(std::exception& e){
+ std::string str = "";
+ str += "update_unnormalized_cdf ";
+ str += e.what();
+ str += " t0 = %1%, t1 = %2%, e0 = %3%, e1 =%4%, ";
+ str += "l1 = %5%, w0 = %6%, w1 = %7%, ";
+ str += "delta = %8%, dy = %9%, ";
+ str += "unnormalized_cdf = %10%";
+ format f(str); f%t0%t1%e0%e1%l1%w0%w1%delta%dy;
+ f%unnormalized_cdf;
+ str = f.str();
+ str+= (this->as_string());
+ throw std::runtime_error(str);
+ }
+ }
/// --> For testing only
- const upper_t&
- range_dy()const{
+ const upper_t& range_dy()const{
static upper_t result;
result.clear();
transform(
@@ -476,9 +656,8 @@
)
);
return result;
- }
- const upper_t&
- range_x()const{
+ }
+ const upper_t& range_x()const{
static upper_t result;
result.clear();
transform(
@@ -491,9 +670,8 @@
)
);
return result;
- }
- const upper_t&
- range_y()const{
+ }
+ const upper_t& range_y()const{
static upper_t result;
result.clear();
transform(
@@ -506,14 +684,16 @@
)
);
return result;
- }
- /// <--- For Testing only
- void assert_upper_points()const{
+ }
+ /// <--- For Testing only
+
+
+ void assert_upper_points()const{
BOOST_ASSERT(
size(upper_) - (size(points_)+1) ==0
);
- }
- void assert_points_sorted()const{
+ }
+ void assert_points_sorted()const{
BOOST_ASSERT(
is_sorted(
begin(points_),
@@ -521,14 +701,14 @@
less<value_type>()
)
);
- }
- void assert_upper_sorted()const{
+ }
+ void assert_upper_sorted()const{
BOOST_ASSERT(
is_sorted(begin(upper_),end(upper_))
);
- }
+ }
- void assert_range_dy_sorted()const{
+ void assert_range_dy_sorted()const{
BOOST_ASSERT(
is_sorted(
begin(range_dy()),
@@ -536,29 +716,27 @@
std::greater<value_type>()
)
);
- }
-
+ }
- void assert_range_unnormalized_cdf_points()const{
+ void assert_range_unnormalized_cdf_points()const{
BOOST_ASSERT(size(range_unnormalized_cdf_)==size(points_));
- }
- void assert_range_unnormalized_cdf_sorted()const{
+ }
+ void assert_range_unnormalized_cdf_sorted()const{
BOOST_ASSERT(
is_sorted(
begin(range_unnormalized_cdf_),
end(range_unnormalized_cdf_)
)
);
- }
+ }
value_type min_;
value_type max_;
points_t points_;
upper_t upper_;
upper_t range_unnormalized_cdf_;
- //value_type unnormalized_cdf_;
- value_type subtract_from_y_;
- };
+ value_type shift_;
+ };
}//adaptive rejection sampling
}//boost
Modified: sandbox/conditionally_specified_distribution/adaptive_rejection_sampling/boost/adaptive_rejection_sampling/example/standard_gaussian.hpp
==============================================================================
--- sandbox/conditionally_specified_distribution/adaptive_rejection_sampling/boost/adaptive_rejection_sampling/example/standard_gaussian.hpp (original)
+++ sandbox/conditionally_specified_distribution/adaptive_rejection_sampling/boost/adaptive_rejection_sampling/example/standard_gaussian.hpp 2009-03-26 18:18:07 EDT (Thu, 26 Mar 2009)
@@ -6,7 +6,7 @@
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
#ifndef BOOST_ADAPTIVE_REJECTION_SAMPLING_STANDARD_GAUSSIAN_HPP_ER_2009
#define BOOST_ADAPTIVE_REJECTION_SAMPLING_STANDARD_GAUSSIAN_HPP_ER_2009
-#include <boost/math_limit/miscellanea.hpp>
+#include <boost/math/tools/precision.hpp>
namespace boost{
namespace adaptive_rejection_sampling{
@@ -14,7 +14,6 @@
/// \brief Example modeling LogConcaveEvaluator
template<typename RealType>
class standard_gaussian_evaluator{
- typedef math_limit::limits<RealType> limits_t;
public:
typedef RealType value_type;
typedef standard_gaussian_evaluator base_type;
@@ -25,10 +24,10 @@
return -x;
}
value_type min()const{
- return limits_t::negative_infinity();
+ return (-max());
}
value_type max()const{
- return limits_t::infinity();
+ return math::tools::max_value<RealType>();
}
};
Modified: sandbox/conditionally_specified_distribution/adaptive_rejection_sampling/boost/adaptive_rejection_sampling/point.hpp
==============================================================================
--- sandbox/conditionally_specified_distribution/adaptive_rejection_sampling/boost/adaptive_rejection_sampling/point.hpp (original)
+++ sandbox/conditionally_specified_distribution/adaptive_rejection_sampling/boost/adaptive_rejection_sampling/point.hpp 2009-03-26 18:18:07 EDT (Thu, 26 Mar 2009)
@@ -6,6 +6,7 @@
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
#ifndef BOOST_ADAPTIVE_REJECTION_SAMPLING_POINT_HPP_ER_2009
#define BOOST_ADAPTIVE_REJECTION_SAMPLING_POINT_HPP_ER_2009
+#include <string>
#include <stdexcept>
#include <boost/format.hpp>
#include <boost/math_limit/miscellanea.hpp>
@@ -23,6 +24,12 @@
value_type y()const{return y_;}
value_type dy()const{return dy_;}
+ std::string as_string()const{
+ std::string str ="x = %1%, y = %2%, dy = %3%";
+ format f(str);
+ f%x()%y()%dy();
+ return f.str();
+ }
private:
value_type x_;
value_type y_;
@@ -31,49 +38,46 @@
template<typename RealType>
RealType
- tangent_at_negative_infinity(const point<RealType>& p){
- typedef math_limit::limits<RealType> limits_t;
-
- RealType result = 0.0;
- if(p.dy()>0.0){
- result = limits_t::negative_infinity();
- }else{
- result = limits_t::infinity();
- }
- return result;
- }
- template<typename RealType>
- RealType
- tangent_at_infinity(const point<RealType>& p){
- typedef math_limit::limits<RealType> limits_t;
- RealType result = 0.0;
- if(p.dy()>0.0){
- result = limits_t::infinity();
- }else{
- result = limits_t::negative_infinity();
+ linear_approximation(
+ const point<RealType>& p0,
+ const point<RealType>& p1,
+ RealType x
+ )
+ { try{
+ if(x<=p1.x()){
+ }else{
+ std::string str = "x<=p1.x()";
+ throw std::runtime_error(str);
+ }
+ if(p0.x()<=x){}else{
+ std::string str = "p0.x()<=x";
+ throw std::runtime_error(str);
+ }
+ if(p0.x()<p1.x()){}else{
+ std::string str = "p0.x()<p1.x()";
+ throw std::runtime_error(str);
+ }
+ }catch(std::exception& e){
+ std::string str = "linear_approximation(p0,p1,x) : ";
+ str+= e.what();
+ str+= " p0 : "; str+= p0.as_string();
+ str+= " p1 : "; str+= p1.as_string();
+ str+= " x = %1%.";
+ format f(str); f%x;
+ throw std::runtime_error(f.str());
}
+ RealType result = (x-p0.x()) * (p1.y()-p0.y());
+ result /= (p1.x()-p0.x());
+ result += p0.y();
return result;
}
- template<typename RealType>
- RealType tangent_at_finite(const point<RealType>& p,RealType x){
- return p.y()+(x-p.x()) * p.dy();
- }
+
template<typename RealType>
RealType tangent(const point<RealType>& p,RealType x){
- typedef math_limit::limits<RealType> limits_t;
- if(limits_t::is_infinity(x)){
- return tangent_at_infinity(p);
- }else{
- if(limits_t::is_negative_infinity(x)){
- return tangent_at_negative_infinity(p);
- }else{
- return tangent_at_finite(p,x);
- }
- }
+ return p.y()+(x-p.x()) * p.dy();
}
-
template<typename RealType>
RealType
tangent_intersection(
Modified: sandbox/conditionally_specified_distribution/adaptive_rejection_sampling/boost/adaptive_rejection_sampling/sampler.hpp
==============================================================================
--- sandbox/conditionally_specified_distribution/adaptive_rejection_sampling/boost/adaptive_rejection_sampling/sampler.hpp (original)
+++ sandbox/conditionally_specified_distribution/adaptive_rejection_sampling/boost/adaptive_rejection_sampling/sampler.hpp 2009-03-26 18:18:07 EDT (Thu, 26 Mar 2009)
@@ -7,6 +7,7 @@
#ifndef BOOST_ADAPTIVE_REJECTION_SAMPLING_SAMPLER_HPP_ER_2009
#define BOOST_ADAPTIVE_REJECTION_SAMPLING_SAMPLER_HPP_ER_2009
#include <vector>
+#include <cmath>
#include <stdexcept>
#include <boost/assert.hpp>
#include <stdexcept>
@@ -95,16 +96,13 @@
u = gen_excl_endpoints(e);
//std::cout << "u1=" << std::setprecision(10) << u;
x = approx_t::inverse_cdf(u);
- //std::cout << "x=" << std::setprecision(10) << x;
thresh = exp_lower_minus_upper(x);
u = gen_excl_endpoints(e);
- //std::cout << "u2=" << std::setprecision(10) << u;
reject = !( u <= thresh);
if( reject ){
thresh = exp_y_minus_upper(x);
u = gen_excl_endpoints(e);
- //std::cout << "u3=" << std::setprecision(10) << u;
reject = !(u<=thresh);
if(reject){
approx_t::update(x);
@@ -121,15 +119,16 @@
result_type postponed_;
result_type
exp_lower_minus_upper(result_type x)const{
- result_type a = approx_t::lower(x); // - inf?!
- a -= approx_t::upper(x);
- return math_limit::safer_exp(a);
+ result_type l = approx_t::lower(x);
+ result_type h = approx_t::upper(x);
+ result_type e = exp(l-h); //safer_exp
+ return e;
}
result_type
exp_y_minus_upper(result_type x)const{
result_type a = approx_t::get_y(x);
a -= approx_t::upper(x);
- return math_limit::safer_exp(a);
+ return exp(a); //safer_exp
}
void postpone_update(result_type x){
is_postponed_ = true;
Modified: sandbox/conditionally_specified_distribution/adaptive_rejection_sampling/libs/adaptive_rejection_sampling/src/doc/readme.txt
==============================================================================
--- sandbox/conditionally_specified_distribution/adaptive_rejection_sampling/libs/adaptive_rejection_sampling/src/doc/readme.txt (original)
+++ sandbox/conditionally_specified_distribution/adaptive_rejection_sampling/libs/adaptive_rejection_sampling/src/doc/readme.txt 2009-03-26 18:18:07 EDT (Thu, 26 Mar 2009)
@@ -85,4 +85,10 @@
At i = 181, normalizing constant = 2.52269 vs true value 2.50663
<-
+///////////
+/ History /
+///////////
+2009/03/18
+Fixed some stability issues pertaining to shift_ and removed the reliance
+on a truncated range e.g. exp(-inf) is not only fine but desirable.
Modified: sandbox/conditionally_specified_distribution/adaptive_rejection_sampling/libs/adaptive_rejection_sampling/src/example/test_data.cpp
==============================================================================
--- sandbox/conditionally_specified_distribution/adaptive_rejection_sampling/libs/adaptive_rejection_sampling/src/example/test_data.cpp (original)
+++ sandbox/conditionally_specified_distribution/adaptive_rejection_sampling/libs/adaptive_rejection_sampling/src/example/test_data.cpp 2009-03-26 18:18:07 EDT (Thu, 26 Mar 2009)
@@ -8,8 +8,11 @@
namespace libs{namespace adaptive_rejection_sampling{
const double limitingNc=2.50663;//sqrt(2*pi)
- const double initAr[rowsCount][2]={{-0.5,1.5},{-0.001,0.001}};
- //TODO {-40,40} or above will generate overflow.
+ const double initAr[rowsCount][2]={
+ {-0.5,1.5},
+ {-0.001,0.001},
+ {-40.0,40.0}
+ };
const double unifsAr[]={
//SeedRandom[0]; Table[Random[], {i, 1, 500}]
@@ -155,7 +158,7 @@
-1.09034, -1.45428, -0.504806, 0.00461931, -1.4878, -0.495555, -0.121771,
-1.2577, 0.204893, 1.33436, 0.934877, -0.475671, 0.486964, -0.694213,
-0.309207, 0.651196, -0.497478, -0.174707, 1.07433, -0.450561
- }/*,
+ },//here
{
0.529408, -0.572962, 0.513134, 1.09102, -0.446931, -2.31367, -1.07024,
-1.65772, -0.991474, -1.41893, 0.316418, 0.101749, 2.20439, -0.789688,
@@ -186,9 +189,9 @@
-1.94589, -0.699194, -1.02964, -0.457675, 0.161575, 0.447034, 0.597158,
1.50555, 1.48887, -0.832003, -0.372902, 0.0313541
}
-*/
+//here
};
-
+/*
const double dbarsAr_1[rowsCount][colsCount]={
{
0.141745, 1.33381, 0.426831, 1.26225, -1.89854, -0.897882, 0.838643,
@@ -253,6 +256,7 @@
}
};
+ */
}//adaptive_rejection_sampling
}//libs
Modified: sandbox/conditionally_specified_distribution/adaptive_rejection_sampling/libs/adaptive_rejection_sampling/src/example/test_data.h
==============================================================================
--- sandbox/conditionally_specified_distribution/adaptive_rejection_sampling/libs/adaptive_rejection_sampling/src/example/test_data.h (original)
+++ sandbox/conditionally_specified_distribution/adaptive_rejection_sampling/libs/adaptive_rejection_sampling/src/example/test_data.h 2009-03-26 18:18:07 EDT (Thu, 26 Mar 2009)
@@ -16,7 +16,7 @@
namespace libs{namespace adaptive_rejection_sampling{
extern const double limitingNc;//sqrt(2*pi)
- const unsigned rowsCount=2;
+ const unsigned rowsCount=3;//2
const unsigned colsCount=200;
const unsigned unifsCount=500;
const unsigned show_count=1; //show every show_count iteration
@@ -37,11 +37,9 @@
Uniform_sampler_mathematica(){
unifs.reserve(unifsCount);
for(size_t i=0; i<unifsCount; ++i){
- //std::cout << "unifsAr[i]=" << unifsAr[i] << std::endl;
unifs.push_back(unifsAr[i]);
};
it = unifs.begin();
- //std::cout << "unifs.size()=" << unifs.size() << std::endl;
};
unsigned long it_distance()const{return (it-unifs.begin());};
double operator()(){
Modified: sandbox/conditionally_specified_distribution/adaptive_rejection_sampling/libs/adaptive_rejection_sampling/src/example/test_gaussian.cpp
==============================================================================
--- sandbox/conditionally_specified_distribution/adaptive_rejection_sampling/libs/adaptive_rejection_sampling/src/example/test_gaussian.cpp (original)
+++ sandbox/conditionally_specified_distribution/adaptive_rejection_sampling/libs/adaptive_rejection_sampling/src/example/test_gaussian.cpp 2009-03-26 18:18:07 EDT (Thu, 26 Mar 2009)
@@ -6,6 +6,7 @@
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
#include <iostream>
#include <list>
+#include <string>
#include <vector>
#include <boost/assert.hpp>
#include <boost/static_assert.hpp>
@@ -18,20 +19,19 @@
namespace adaptive_rejection_sampling{
void example_standard_gaussian(){
std::cout<<"-> test_gaussian"<<std::endl;
- const double tolerance = 1e-4; //1e-5 causes fail
- const unsigned report_every = 20;
+ const double tolerance = 1e-4;
std::string str;
- str+= "We draw from the standard-gaussian using the algorithm.";
- str+= " Results are compared with a similar implementation in";
- str+= "Mathematica. A failure will be generated if the difference";
- str+= " is greater than %1%";
+ str+= "n = %1% samples from N(0,1) using adaptive rejection";
+ str+= "sampling vs a similar implementation in Mathematica.";
+ str+= " Tolerance is set at %2%.";
str+= " The uniform draws were Mathematica generated.";
- boost::format f(str); f % tolerance;
+ boost::format f(str); f % colsCount % tolerance;
std::cout << f.str() << std::endl;
Uniform_sampler_mathematica usm;
+ typedef double value_type;
typedef boost::adaptive_rejection_sampling
- ::standard_gaussian_evaluator<double>
+ ::standard_gaussian_evaluator<value_type>
dist_fun_type;
typedef boost::adaptive_rejection_sampling
::sampler<dist_fun_type,std::vector> dbars_t;
@@ -47,33 +47,33 @@
double drawVal,ncVal;
BOOST_STATIC_ASSERT(colsCount<=unifsCount);
usm.reset();
- dbars.initialize(initAr[j][0], initAr[j][1]);
std::cout<<"Initialized with (x1,x2)=("<<
- initAr[j][0]<<","<<initAr[j][1]<<")"<<std::endl;
+ initAr[j][0]<<","<<initAr[j][1]<<") ";
+ dbars.initialize(initAr[j][0], initAr[j][1]);
for(unsigned int i=0; i<colsCount; i++){
- //std::cout << "usm.it_distance()="
- //<< usm.it_distance() << std::endl;
+ //std::cout << "i=" << i << std::endl;
drawVal = dbars(usm);
ncVal = dbars.total_unnormalized_cdf();
- BOOST_ASSERT(drawVal < dbarsAr[j][i]+tolerance);
- BOOST_ASSERT(drawVal > dbarsAr[j][i]-tolerance);
- if(!(i%report_every)){
- std::cout
- << "At i = "<< i+1
- << ", normalizing constant = "
- << ncVal
- << " vs true value "
- << limitingNc
- << std::endl;
- }
+ if(
+ (drawVal > dbarsAr[j][i]+tolerance)
+ ||
+ (drawVal < dbarsAr[j][i]-tolerance)
+ )
+ {
+ std::string str = "(drawVal = %1%)";
+ str+= "- (dbarsAr[j][i] = %2% ) < %3%";
+ boost::format f(str); f%drawVal%dbarsAr[j][i];
+ f % tolerance;
+ throw std::runtime_error(f.str());
+ }
};
+ std::cout << ": OK." << std::endl;
};
}catch(std::exception& e){
- std::cerr << e.what() << std::endl;
- std::cout << "This throw is meant to illustrate the"
- << "!(delta>0.0) pitfall described in"
- << " approximation_impl::update_unnormalized_cdf_impl"
- << std::endl;
+ std::string str = e.what();
+ str+= " ";
+ str+= dbars.as_string();
+ std::cerr << str << std::endl;
};
std::cout<<"<-"<<std::endl;
Modified: sandbox/conditionally_specified_distribution/adaptive_rejection_sampling/libs/adaptive_rejection_sampling/src/main.cpp
==============================================================================
--- sandbox/conditionally_specified_distribution/adaptive_rejection_sampling/libs/adaptive_rejection_sampling/src/main.cpp (original)
+++ sandbox/conditionally_specified_distribution/adaptive_rejection_sampling/libs/adaptive_rejection_sampling/src/main.cpp 2009-03-26 18:18:07 EDT (Thu, 26 Mar 2009)
@@ -5,6 +5,7 @@
// Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
#include <iostream>
+#include <cmath>
#include <boost/math_limit/miscellanea.hpp>
#include <libs/adaptive_rejection_sampling/src/example/test_gaussian.h>
@@ -12,23 +13,37 @@
int main()
{
+ using namespace boost;
- typedef boost::math_limit::limits<double> tmp_t;
- std::cout
- << "limits::negative_infinity()"
- << tmp_t::negative_infinity() << std::endl;
- std::cout
- << "limits::infinity()"
- << tmp_t::infinity() << std::endl;
- std::cout
- << "limits::min_value()"
- << tmp_t::min_value() << std::endl;
- std::cout
- << "limits::log_min_value()"
- << tmp_t::log_min_value() << std::endl;
- std::cout
- << "limits::log_max_value()"
- << tmp_t::log_max_value() << std::endl;
+ typedef double d_t;
+ typedef long double ld_t;
+ typedef boost::math_limit::exp_support<d_t> exp_dsup_t;
+ typedef boost::math_limit::exp_support<ld_t> exp_ldsup_t;
+ d_t d_es_l = exp_dsup_t::lowest();
+ d_t d_es_h = exp_dsup_t::highest();
+ ld_t ld_es_l = exp_ldsup_t::lowest();
+ ld_t ld_es_h = exp_ldsup_t::highest();
+ d_t d_h = math::tools::max_value<d_t>();
+ d_t d_l = - d_h;
+ ld_t ld_h = math::tools::max_value<ld_t>();
+ ld_t ld_l = - ld_h;
+
+ d_t d_e_es_l = math_limit::truncated_exp(d_es_l);
+ d_t d_e_es_h = math_limit::truncated_exp(d_es_h);
+ ld_t ld_e_es_l = math_limit::truncated_exp(ld_es_l);
+ ld_t ld_e_es_h = math_limit::truncated_exp(ld_es_h);
+ d_t d_e_l = math_limit::truncated_exp(d_l);
+ d_t d_e_h = math_limit::truncated_exp(d_h);
+ ld_t ld_e_l = math_limit::truncated_exp(ld_l);
+ ld_t ld_e_h = math_limit::truncated_exp(ld_h);
+
+ const d_t dtol = 1e-10;
+ const ld_t ldtol = 1e-10;
+
+ BOOST_ASSERT(fabs(d_e_es_l-d_e_l)<dtol);
+ BOOST_ASSERT(fabs(d_e_es_h-d_e_h)<dtol);
+ BOOST_ASSERT(fabs(ld_e_es_l-ld_e_l)<ldtol);
+ BOOST_ASSERT(fabs(ld_e_es_h-ld_e_h)<ldtol);
libs::adaptive_rejection_sampling::example_standard_gaussian();
return 0;
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