Subject: [Boost-commit] svn:boost r52002 - in sandbox/conditionally_specified_distribution: conditionally_specified_distribution/boost/conditionally_specified_distribution/function conditionally_specified_distribution/boost/conditionally_specified_distribution/parameter conditionally_specified_distribution/boost/conditionally_specified_distribution/parameter/dependent conditionally_specified_distribution/boost/conditionally_specified_distribution/parameter/dependent/detail conditionally_specified_distribution/boost/conditionally_specified_distribution/sampler/adaptive_rejection conditionally_specified_distribution/boost/conditionally_specified_distribution/support conditionally_specified_distribution/boost/conditionally_specified_distribution/support/detail conditionally_specified_distribution/libs/conditionally_specified_distribution/src conditionally_specified_distribution/libs/conditionally_specified_distribution/src/doc conditionally_specified_distribution/libs/conditionally_specified_distribution/src/example/sampler math_limit/boost/math_limit
From: erwann.rogard_at_[hidden]
Date: 2009-03-26 18:24:28
Author: e_r
Date: 2009-03-26 18:24:23 EDT (Thu, 26 Mar 2009)
New Revision: 52002
URL: http://svn.boost.org/trac/boost/changeset/52002
Log:
Fixed bug support(), added as_string() and misc in libs
Added:
sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/support/pair_independent.hpp (contents, props changed)
Text files modified:
sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/function/normal.hpp | 1
sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/parameter/dependent/detail/regression_coefficient.hpp | 28 ++++++-
sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/parameter/dependent/normal_given_normal.hpp | 20 +++++
sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/parameter/normal.hpp | 18 ++++
sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/parameter/pair_independent.hpp | 16 ++++
sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/sampler/adaptive_rejection/sampler.hpp | 57 ++++++++++++----
sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/support/detail/negative_max_to_max.hpp | 8 -
sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/support/include.hpp | 5 +
sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/support/normal.hpp | 9 ++
sandbox/conditionally_specified_distribution/conditionally_specified_distribution/libs/conditionally_specified_distribution/src/doc/readme.txt | 31 +++++++-
sandbox/conditionally_specified_distribution/conditionally_specified_distribution/libs/conditionally_specified_distribution/src/example/sampler/gibbs_ars_bivariate_normal.cpp | 64 +++++++++++++----
sandbox/conditionally_specified_distribution/conditionally_specified_distribution/libs/conditionally_specified_distribution/src/main.cpp | 2
sandbox/conditionally_specified_distribution/math_limit/boost/math_limit/miscellanea.hpp | 139 +++++++++++++++++++++++++++++----------
13 files changed, 307 insertions(+), 91 deletions(-)
Modified: sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/function/normal.hpp
==============================================================================
--- sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/function/normal.hpp (original)
+++ sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/function/normal.hpp 2009-03-26 18:24:23 EDT (Thu, 26 Mar 2009)
@@ -9,6 +9,7 @@
#define BOOST_CONDITIONALLY_SPECIFIED_DISTRIBUTION_FUNCTION_NORMAL_HPP_ER_2009
#include <cmath>
#include <boost/math/special_functions/fpclassify.hpp>
+#include <boost/math/constants/constants.hpp>
namespace boost{
namespace conditionally_specified_distribution{
namespace function{
Modified: sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/parameter/dependent/detail/regression_coefficient.hpp
==============================================================================
--- sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/parameter/dependent/detail/regression_coefficient.hpp (original)
+++ sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/parameter/dependent/detail/regression_coefficient.hpp 2009-03-26 18:24:23 EDT (Thu, 26 Mar 2009)
@@ -8,8 +8,10 @@
#ifndef BOOST_CONDITIONALLY_SPECIFIED_DISTRIBUTION_PARAMETER_DEPENDENT_REGRESSION_COEFFICIENT_HPP_ER_2009
#define BOOST_CONDITIONALLY_SPECIFIED_DISTRIBUTION_PARAMETER_DEPENDENT_REGRESSION_COEFFICIENT_HPP_ER_2009
#include <vector>
+#include <string>
#include <boost/range.hpp>
#include <boost/bind.hpp>
+#include <boost/format.hpp>
#include <boost/type_traits.hpp>
#include <boost/mpl/assert.hpp>
#include <boost/mpl/identity.hpp>
@@ -18,6 +20,7 @@
#include <boost/utility/dont_care.hpp>
#include <boost/utility/remove_qualifier.hpp>
#include <boost/utility/assert_is_base_of.hpp>
+#include <boost/utility/container_to_string.hpp>
#include <boost/shared_features/feature/keyword.hpp>
#include <boost/shared_features/parameter.hpp>
#include <boost/shared_features/depends_on.hpp>
@@ -86,17 +89,14 @@
>::type
>::type value_type;
typedef std::vector<value_type> container_type;
-
public:
struct result_of_range_skip_position_covariate_dot_coefficient
: mpl::identity<const container_type&>
{};
-
// outside interface
struct function_value : mpl::identity<value_type>{};
struct function_argument : mpl::identity<value_type>{};
};
-
template<typename Id, typename Dataset,typename Coefficients>
class regression_coefficient
: public regression_coefficient_base<
@@ -115,8 +115,6 @@
:position_(pos){
//Warning : state_ default initializes
}
-
-
template<typename Args>
regression_coefficient(const Args& args)
:position_(args[shared_features::kwd<Id>::position|0]){
@@ -188,6 +186,8 @@
return this->matching_covariate(u);
}
+ //TODO range_matching_covariate
+
// {<x[-i],b[j]> : i =1,...,n}
typename conditionally_specified_distribution::result_of
::range_skip_position_covariate_dot_coefficient<
@@ -222,6 +222,24 @@
the_op
);
}
+
+ template<typename Args>
+ std::string
+ as_string(const Args& args)const{
+ typedef utility::container_to_string to_str_t;
+ std::string str = "[position = %1%\n";
+ format f(str);
+ f%(this->position(args));
+ str = f.str();
+ str += "range_coefficient =";
+ str += to_str_t::get_indexed(this->range_coefficient(args));
+ str += "\nrange_skip_position_covariate_dot_coefficient = ";
+ str += to_str_t::get_indexed(
+ this->range_skip_position_covariate_dot_coefficient(args)
+ );
+ str +="]";
+ return str;
+ }
private:
void alert_if(mpl::bool_<true>,int i)const{
BOOST_ASSERT(i>-1);
Modified: sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/parameter/dependent/normal_given_normal.hpp
==============================================================================
--- sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/parameter/dependent/normal_given_normal.hpp (original)
+++ sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/parameter/dependent/normal_given_normal.hpp 2009-03-26 18:24:23 EDT (Thu, 26 Mar 2009)
@@ -8,6 +8,8 @@
#ifndef BOOST_CONDITIONALLY_SPECIFIED_DISTRIBUTION_PARAMETER_NORMAL_GIVEN_NORMAL_HPP_ER_2009
#define BOOST_CONDITIONALLY_SPECIFIED_DISTRIBUTION_PARAMETER_NORMAL_GIVEN_NORMAL_HPP_ER_2009
#include <cmath>
+#include <string>
+#include <boost/format.hpp>
#include <boost/mpl/apply.hpp>
#include <boost/mpl/vector.hpp>
#include <boost/mpl/identity.hpp>
@@ -124,10 +126,28 @@
typename super_t::result_of_beta::type
beta(utility::dont_care)const{ return this->beta(); }
+ template<typename Args>
+ std::string
+ as_string(const Args& args)const{
+ return this->as_string_impl(
+ args[::boost::shared_features::kwd_set]);
+ }
+
private:
typename super_t::value_type mu_;
typename super_t::value_type sigma_;
typename super_t::value_type beta_;
+
+ template<typename Set>
+ std::string
+ as_string_impl(const Set& set)const{
+ std::string str = "parameter::normal_given_normal (y|x)";
+ str+= "mu = %1%, sigma = %2%";
+ str+= "x= %3% ";
+ format f(str);
+ f%(this->mu())%(this->sigma())%((this->get_x(set))());
+ return f.str();
+ }
template<typename Set>
typename super_t::value_type sigma(const Set& set){
typename super_t::value_type sigma_y = get_par_y(set).sigma();
Modified: sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/parameter/normal.hpp
==============================================================================
--- sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/parameter/normal.hpp (original)
+++ sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/parameter/normal.hpp 2009-03-26 18:24:23 EDT (Thu, 26 Mar 2009)
@@ -7,9 +7,12 @@
////////////////////////////////////////////////////////////////////////////
#ifndef BOOST_CONDITIONALLY_SPECIFIED_DISTRIBUTION_PARAMETER_NORMAL_HPP_ER_2009
#define BOOST_CONDITIONALLY_SPECIFIED_DISTRIBUTION_PARAMETER_NORMAL_HPP_ER_2009
+#include <string>
#include <boost/mpl/apply.hpp>
#include <boost/parameter/keyword.hpp>
+#include <boost/format.hpp>
#include <boost/utility/assert_is_base_of.hpp>
+#include <boost/utility/dont_care.hpp>
#include <boost/conditionally_specified_distribution/keyword/parameter.hpp>
#include <boost/conditionally_specified_distribution/result_of/include.hpp>
#include <boost/conditionally_specified_distribution/crtp/include.hpp>
@@ -76,8 +79,7 @@
return *this;
}
- template<typename Args>
- void update(const Args& args){}
+ void update(utility::dont_care){}
typename result_of::mu<super_t>::type
mu() const {
@@ -97,6 +99,18 @@
typename result_of::sigma<super_t>::type
sigma(const Args& args) const { return this->sigma(); }
+ std::string as_string()const{
+ std::string str = "parameter::impl::normal : ";
+ str += " mu = %1%";
+ str += " sigma = %2%";
+ format f(str); f%mu_%sigma_;
+ return f.str();
+ }
+
+ std::string as_string(utility::dont_care)const{
+ return this->as_string();
+ }
+
private:
typename super_t::value_type mu_;
typename super_t::value_type sigma_;
Modified: sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/parameter/pair_independent.hpp
==============================================================================
--- sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/parameter/pair_independent.hpp (original)
+++ sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/parameter/pair_independent.hpp 2009-03-26 18:24:23 EDT (Thu, 26 Mar 2009)
@@ -7,6 +7,7 @@
////////////////////////////////////////////////////////////////////////////
#ifndef BOOST_CONDITIONALLY_SPECIFIED_DISTRIBUTION_PAIR_INDEPENDENT_HPP_ER_2009
#define BOOST_CONDITIONALLY_SPECIFIED_DISTRIBUTION_PAIR_INDEPENDENT_HPP_ER_2009
+#include <string>
#include <boost/mpl/apply.hpp>
#include <boost/mpl/vector.hpp>
#include <boost/parameter/keyword.hpp>
@@ -85,7 +86,7 @@
return *this;
}
- template<typename Set,typename Args>
+ template<typename Args>
void update(const Args& args){
this->update_impl(args[shared_features::kwd_set],args);
}
@@ -101,6 +102,19 @@
second(const Args& args) const {
return second_impl(args[shared_features::kwd_set],args);
}
+
+ template<typename Args>
+ std::string as_string(const Args& args)const{
+ std::string str = "parameter::pair_independent :";
+ str+= "\n first : ";
+ str+= first_impl(
+ args[shared_features::kwd_set],args).as_string(args);
+ str+= "\n second : ";
+ str+= second_impl(
+ args[shared_features::kwd_set],args).as_string(args);
+ return str;
+ }
+
private:
template<typename Set,typename Args>
typename result_of::first<super_t(const Args&)>::type
Modified: sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/sampler/adaptive_rejection/sampler.hpp
==============================================================================
--- sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/sampler/adaptive_rejection/sampler.hpp (original)
+++ sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/sampler/adaptive_rejection/sampler.hpp 2009-03-26 18:24:23 EDT (Thu, 26 Mar 2009)
@@ -35,6 +35,19 @@
namespace impl{
+ template<template<typename,typename> class Cont>
+ struct apply_reserve{
+ template<typename Impl,typename T>
+ void operator()(Impl& impl, T n){}
+ };
+ template<>
+ struct apply_reserve<std::vector>{
+ template<typename Impl,typename T>
+ void operator()(Impl& impl, T n){
+ impl.reserve(n);
+ }
+ };
+
template<typename Par>
class adaptive_rejection_sampler_base{
public:
@@ -67,10 +80,13 @@
typedef kwd<> kwd;
std::size_t static default_max_init_recursion(){
- std::size_t static result = 100;
+ static std::size_t result = 100;
+ return result;
+ }
+ std::size_t static default_reserve(){
+ static std::size_t result = 10;
return result;
}
-
public:
template<typename Args>
@@ -78,11 +94,13 @@
: init_(args),
max_init_recursion_(
args[kwd::ars_max_init_recursion|default_max_init_recursion()]
- ){}
+ ),
+ reserve_(default_reserve()){}
adaptive_rejection_sampler(const adaptive_rejection_sampler& that)
: init_(that.init_),
- max_init_recursion_(that.max_init_recursion_){}
+ max_init_recursion_(that.max_init_recursion_),
+ reserve_(that.reserve_){}
//TODO check
adaptive_rejection_sampler&
@@ -90,6 +108,7 @@
if(&that!=this){
init_ = that.init_;
max_init_recursion_ = that.max_init_recursion_;
+ reserve_ = that.reserve_;
}
return *this;
}
@@ -127,11 +146,6 @@
typedef typename
Set::template result_of_extract<State>::type extr_state_t;
typedef function::detail::adapter<extr_par_t,Args> adapter_t;
-// typedef ::boost::random::adaptive_rejection_sampler::simulator<
-// adapter_t,
-// Container,
-// Allocator
-// > impl_t;
typedef ::boost::adaptive_rejection_sampling::sampler<
adapter_t,
Container,
@@ -139,12 +153,24 @@
> impl_t;
const extr_par_t& par = set.template extract<Par>();
adapter_t adapter (par,args);
+ //TODO (*)
+ //static impl_t impl(
+ // adapter,
+ // this->max_init_recursion_
+ //);
+
impl_t impl(
adapter,
this->max_init_recursion_
);
+
interval_t interval = (this->init_)();
- // Have had px!=0 fails here with gibbs rte::regression_coeff
+ // TODO (*)
+ // This causes problem due I think to dependence adapter
+ // on Args (wrap around an abstract base class?)
+ // impl.reset_distribution_function(adapter);
+ typedef apply_reserve<Container> apply_res_t;
+ apply_res_t()(impl,reserve_);
impl.initialize(interval.first,interval.second);
gen_t gen(urng,unif01_t());
@@ -153,17 +179,20 @@
extr_state_t& extr_state = set.template extract<State>();
extr_state.set(args,x);
(this->init_).update(impl);
+ if(reserve_<impl.size_range_point()){
+ reserve_ = impl.size_range_point();
+ }
}
mutable init_t init_;
typename super_t::size_type max_init_recursion_;
+ mutable typename super_t::size_type reserve_;
};
}//impl
namespace feature{
- // Warning : pair_quantile_averaged_over_iterations is non_markov
- // Consequential theretically? Dunno, but from a practical
- // standpoint, less likely to lead to stability issues
- // such as exp(-inf).
+ // Warning : pair_quantile_averaged_over_iterations is likely
+ // to require fewer iterations than pair_quantile.
+ // Non-markovian : consequential for convergence?
template<
typename Par,
Modified: sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/support/detail/negative_max_to_max.hpp
==============================================================================
--- sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/support/detail/negative_max_to_max.hpp (original)
+++ sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/support/detail/negative_max_to_max.hpp 2009-03-26 18:24:23 EDT (Thu, 26 Mar 2009)
@@ -21,13 +21,11 @@
support_negative_max_to_max(
const D& d
){
- typedef typename function_support<D>::type fs;
+ typedef typename function_support<D>::type sup_t;
+ typedef typename function_argument<D>::type arg_t;
using boost::math::tools::max_value;
- return fs(
- -max_value<typename fs::first_type>(),
- max_value<typename fs::second_type>()
- );
+ return sup_t(-max_value<arg_t>(),max_value<arg_t>());
};
}//detail
Modified: sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/support/include.hpp
==============================================================================
--- sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/support/include.hpp (original)
+++ sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/support/include.hpp 2009-03-26 18:24:23 EDT (Thu, 26 Mar 2009)
@@ -8,7 +8,10 @@
#ifndef BOOST_CONDITIONALLY_SPECIFIED_DISTRIBUTION_SUPPORT_INCLUDE_HPP_ER_2009
#define BOOST_CONDITIONALLY_SPECIFIED_DISTRIBUTION_SUPPORT_INCLUDE_HPP_ER_2009
+#include <boost/conditionally_specified_distribution/support/pair_independent.hpp>
#include <boost/conditionally_specified_distribution/support/normal.hpp>
-#include <boost/conditionally_specified_distribution/support/fwd.hpp>
+
+// Bad idea:
+//#include <boost/conditionally_specified_distribution/support/fwd.hpp>
#endif
Modified: sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/support/normal.hpp
==============================================================================
--- sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/support/normal.hpp (original)
+++ sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/support/normal.hpp 2009-03-26 18:24:23 EDT (Thu, 26 Mar 2009)
@@ -7,6 +7,7 @@
////////////////////////////////////////////////////////////////////////////
#ifndef BOOST_CONDITIONALLY_SPECIFIED_DISTRIBUTION_SUPPORT_NORMAL_HPP_ER_2009
#define BOOST_CONDITIONALLY_SPECIFIED_DISTRIBUTION_SUPPORT_NORMAL_HPP_ER_2009
+#include <boost/utility/dont_care.hpp>
#include <boost/conditionally_specified_distribution/result_of/function_support.hpp>
#include <boost/conditionally_specified_distribution/crtp/normal.hpp>
#include <boost/conditionally_specified_distribution/support/detail/negative_max_to_max.hpp>
@@ -24,6 +25,14 @@
return detail::support_negative_max_to_max(d_);
}
+template <typename D>
+inline
+typename function_support<D>::type
+support(const crtp::normal<D>& d,utility::dont_care)
+{
+ return support(d);
+}
+
}
}
}
Added: sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/support/pair_independent.hpp
==============================================================================
--- (empty file)
+++ sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/support/pair_independent.hpp 2009-03-26 18:24:23 EDT (Thu, 26 Mar 2009)
@@ -0,0 +1,34 @@
+/////////////////////////////////////////////////////////////////////////////
+// pair_independent.hpp
+//
+// Copyright 2009 Erwann Rogard. Distributed under the Boost
+// Software License, Version 1.0. (See accompanying file
+// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+////////////////////////////////////////////////////////////////////////////
+#ifndef BOOST_CONDITIONALLY_SPECIFIED_DISTRIBUTION_SUPPORT_PAIR_INDEPENDENT_HPP_ER_2009
+#define BOOST_CONDITIONALLY_SPECIFIED_DISTRIBUTION_SUPPORT_PAIR_INDEPENDENT_HPP_ER_2009
+#include <boost/conditionally_specified_distribution/result_of/function_support.hpp>
+#include <boost/conditionally_specified_distribution/crtp/pair_independent.hpp>
+namespace boost{
+namespace conditionally_specified_distribution{
+namespace function{
+
+template <typename D,typename Args>
+inline
+typename function_support<D>::type
+support(const crtp::pair_independent<D>& d,const Args& args)
+{
+ const D& d_ = static_cast<const D&>(d);
+
+ // TODO :
+ // pdf(x) = pdf_1(x)pdf_2(x) = 0 if either pdf_1(x) = 0 or pdf_2(y) =0
+ // so support (x,y) is intersection(suppor_1,support_2)
+
+ // For now, assume same support:
+ return support(d_.first(args));
+}
+
+}
+}
+}
+#endif // BOOST_CONDITIONALLY_SPECIFIED_DISTRIBUTION_SUPPORT_PAIR_INDEPENDENT_HPP_ER_2009
Modified: sandbox/conditionally_specified_distribution/conditionally_specified_distribution/libs/conditionally_specified_distribution/src/doc/readme.txt
==============================================================================
--- sandbox/conditionally_specified_distribution/conditionally_specified_distribution/libs/conditionally_specified_distribution/src/doc/readme.txt (original)
+++ sandbox/conditionally_specified_distribution/conditionally_specified_distribution/libs/conditionally_specified_distribution/src/doc/readme.txt 2009-03-26 18:24:23 EDT (Thu, 26 Mar 2009)
@@ -229,7 +229,6 @@
////////////////////////
/ Output from main.cpp /
////////////////////////
-
->example_parameter_normal_given_normal
mu_y|x=0.1
beta_y|x=0.75
@@ -276,12 +275,32 @@
1000000 samples from (x,y) jointly normal, mu_x = -0.1 sigma_x = 1
rho_xy = 0.9 mu_y = 0.1 sigma_y = 1.5 using the decomposition (y|x, x|y)
and an ars sampler for each component.
+i=0
+i=100000
+i=200000
+i=300000
+i=400000
+i=500000
+i=600000
+i=700000
+i=800000
+i=900000
sample statistics :
- mu_x : -0.0976351
- sigma_x : 0.972748
- mu_y : 0.103237
- sigma_y : 1.45947
- rho_xy :0.896274
+ mu_x : -0.0984778
+ sigma_x : 0.974992
+ mu_y : 0.102567
+ sigma_y : 1.46175
+ rho_xy :0.896784
<-
->example_function_pair_independent
OK <-
+
+
+///////////
+/ History /
+///////////
+
+On 2009/03/18, added:
+
+parameter::dependent::detail::regression_coefficient::as_string(args)
+
Modified: sandbox/conditionally_specified_distribution/conditionally_specified_distribution/libs/conditionally_specified_distribution/src/example/sampler/gibbs_ars_bivariate_normal.cpp
==============================================================================
--- sandbox/conditionally_specified_distribution/conditionally_specified_distribution/libs/conditionally_specified_distribution/src/example/sampler/gibbs_ars_bivariate_normal.cpp (original)
+++ sandbox/conditionally_specified_distribution/conditionally_specified_distribution/libs/conditionally_specified_distribution/src/example/sampler/gibbs_ars_bivariate_normal.cpp 2009-03-26 18:24:23 EDT (Thu, 26 Mar 2009)
@@ -8,8 +8,10 @@
#include <numeric>
#include <vector>
#include <cmath>
+#include <string>
#include <boost/assert.hpp>
#include <boost/bind.hpp>
+#include <boost/format.hpp>
#include <boost/random/mersenne_twister.hpp>
#include <boost/type_traits.hpp>
#include <boost/mpl/size_t.hpp>
@@ -24,6 +26,7 @@
#include <boost/shared_features/feature/keyword.hpp>
#include <boost/shared_features/set.hpp>
+#include <boost/shared_features/parameter.hpp>
#include <boost/shared_features/mpl_features.hpp>
#include <boost/shared_features/feature/scalar.hpp>
#include <boost/shared_features/functor/sample.hpp>
@@ -40,7 +43,7 @@
std::cout << "->example_sampler_gibbs_ars_bivariate_normal" << std::endl;
using namespace boost;
- namespace shared = shared_features;
+ namespace shared = boost::shared_features;
namespace csd = conditionally_specified_distribution;
namespace par = csd::parameter;
namespace sampl = csd::sampler;
@@ -186,24 +189,51 @@
std::cout << the_format.str() << std::endl;
- for(unsigned i=0; i<n; ++i){
- //if(i%1e4==0){
- // std::cout << "i=" << i << std::endl;
- //}
- set.visit_if<is_sampler_x>(
- shared::functor::sample(),
- (sampl::kwd<>::random_number_generator = urng)
+ try{
+ for(unsigned i=0; i<n; ++i){
+ unsigned mod = 1e5;
+ if(i%mod==0){
+ std::cout << "i=" << i << std::endl;
+ }
+ try{
+ set.visit_if<is_sampler_x>(
+ shared::functor::sample(),
+ (sampl::kwd<>::random_number_generator = urng)
+ );
+ }catch(std::exception& e){
+ std::string str = "x|y : ";
+ str+= e.what();
+ throw std::runtime_error(str);
+ }
+ try{
+ draws_x.push_back(
+ set.extract<state_x>()()
+ );
+ set.visit_if<is_sampler_y>(
+ shared::functor::sample(),
+ (sampl::kwd<>::random_number_generator = urng)
+ );
+ }catch(std::exception& e){
+ std::string str = "y|x : ";
+ str+= e.what();
+ throw std::runtime_error(str);
+ }
+ draws_y.push_back(
+ set.extract<state_y>()()
+ );
+ }
+ }catch(std::exception& e){
+ std::string str = e.what();
+ str += "y|x=";
+ str += set.extract<par_y_given_x>().as_string(
+ (shared::kwd_set = set)
);
- draws_x.push_back(
- set.extract<state_x>()()
- );
- set.visit_if<is_sampler_y>(
- shared::functor::sample(),
- (sampl::kwd<>::random_number_generator = urng)
- );
- draws_y.push_back(
- set.extract<state_y>()()
+ str += "x|y=";
+ str += set.extract<par_x_given_y>().as_string(
+ (shared::kwd_set = set)
);
+
+ std::cout << str << std::endl;
}
std::cout << "sample statistics : " << std::endl;
Modified: sandbox/conditionally_specified_distribution/conditionally_specified_distribution/libs/conditionally_specified_distribution/src/main.cpp
==============================================================================
--- sandbox/conditionally_specified_distribution/conditionally_specified_distribution/libs/conditionally_specified_distribution/src/main.cpp (original)
+++ sandbox/conditionally_specified_distribution/conditionally_specified_distribution/libs/conditionally_specified_distribution/src/main.cpp 2009-03-26 18:24:23 EDT (Thu, 26 Mar 2009)
@@ -19,7 +19,6 @@
{
example_parameter_normal_given_normal();
-
//normal (from simple to complex):
example_sampler_exact_normal_given_normal();
example_sampler_exact_bivariate_normal();
@@ -29,6 +28,5 @@
example_function_pair_independent();
-
return 0;
}
Modified: sandbox/conditionally_specified_distribution/math_limit/boost/math_limit/miscellanea.hpp
==============================================================================
--- sandbox/conditionally_specified_distribution/math_limit/boost/math_limit/miscellanea.hpp (original)
+++ sandbox/conditionally_specified_distribution/math_limit/boost/math_limit/miscellanea.hpp 2009-03-26 18:24:23 EDT (Thu, 26 Mar 2009)
@@ -6,56 +6,119 @@
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
#ifndef BOOST_MATH_LIMITS_MISCELLANEA_HPP_ER_2009
#define BOOST_MATH_LIMITS_MISCELLANEA_HPP_ER_2009
+#include <iostream>
+#include <cmath>
+#include <stdexcept>
+#include <boost/assert.hpp>
+#include <boost/format.hpp>
+#include <boost/mpl/identity.hpp>
+#include <boost/static_assert.hpp>
+#include <boost/math/special_functions/next.hpp>
#include <boost/math/tools/precision.hpp>
+#include <boost/numeric/conversion/converter.hpp>
+#include <boost/math/special_functions/fpclassify.hpp>
namespace boost{
namespace math_limit{
+ //Discussed in
+ //comp.lang.c++/browse_thread/thread/976701405e6ad0cf?hl=en#
template<typename RealType>
- struct limits{
- static RealType log_max_value(){
- return math::tools::log_max_value<RealType>();
+ struct exp_support{
+ static RealType highest(){
+ struct check{
+ // Warning: fails with RealType = long double with
+ // i686-apple-darwin9-gcc-4.0.1 (GCC) 4.0.1
+ check(RealType x){
+ BOOST_ASSERT(!math::isinf(exp(x)));
+ }
+ };
+ static RealType result
+ = math::tools::log_max_value<RealType>();
+ static check const the_check( result );
+ return result;
}
- static RealType log_min_value(){
- return math::tools::log_min_value<RealType>();
- }
- static RealType negative_infinity(){
- return -infinity();
- }
- static RealType min_value(){
- //beware math::tools::min_value<RealType>() is the
- // min positive value (for floats);
- return math::tools::min_value<RealType>();
- }
- static RealType infinity(){
- return math::tools::max_value<RealType>();
- }
- static RealType is_negative_infinity(RealType x){
- return !(x>negative_infinity());
- }
- static RealType is_infinity(RealType x){
- return !(x<infinity());
+
+ static RealType lowest(){
+ struct check{
+ //TODO is this what needs to be checked?
+ check(RealType x){
+ BOOST_ASSERT(!math::isnan(exp(x)));
+ }
+ };
+ static RealType result
+ = math::tools::log_min_value<RealType>();
+ static check const the_check( result );
+ return result;
}
};
- template<typename RealType>
- RealType
- safer_exp(RealType a){
- typedef math_limit::limits<RealType> limits_t;
- RealType e = 0.0;
- RealType lmx = limits_t::log_max_value();
- RealType lmi = limits_t::log_min_value();
- if(a < lmx){
- if(a > lmi){
- e = exp(a); // -inf < a < inf
+
+ // Warning : if possible avoid truncated expression
+ // Better to dispatch on inf
+ template<typename RealType>
+ RealType
+ truncated_exp(RealType a){
+ typedef math_limit::exp_support<RealType> sup_t;
+ RealType e = 0.0;
+ RealType h = sup_t::highest();
+ RealType l = sup_t::lowest();
+ if(a < h){
+ if(a > l){
+ e = exp(a); // -inf < a < inf
+ }else{
+ e = 0.0;
+ }
}else{
- e = 0.0;
+ e = math::tools::max_value<RealType>(); // a = inf
}
- }else{
- e = limits_t::infinity(); // a = inf
+ BOOST_ASSERT(!std::isnan(e));
+ return e;
}
- BOOST_ASSERT(!std::isnan(e));
- return e;
- }
+
+
+ struct scale_to_max{
+ //news://news.gmane.org:119/gq3755$vb9$1@ger.gmane.org
+ template<class T>// (x/c)+(y/c) <= max
+ static T plus(
+ const T& x,
+ const T& y
+ ) {
+ double max = (std::numeric_limits<T>::max)();
+ if(x > max - y) {
+ double c = x/max + y/max;
+ while(!boost::math::isinf(x/c + y/c)) {
+ c = boost::math::float_prior(c);
+ }
+ while(boost::math::isinf(x/c + y/c)) {
+ c = boost::math::float_next(c);
+ }
+ return(c);
+ } else {
+ return(static_cast<T>(1));
+ }
+ }
+ template<class T>// (x/c)+(y/c) <= max
+ static T divides(
+ const T& x,
+ const T& y
+ ) {
+ double max = (std::numeric_limits<T>::max)();
+ if(x / max > y) {
+ double c = (x/max)/y;
+ while(!boost::math::isinf((x/c)/y)) {
+ c = boost::math::float_prior(c);
+ }
+ while(boost::math::isinf((x/c)/y)) {
+ c = boost::math::float_next(c);
+ }
+ return(c);
+ } else {
+ return(static_cast<T>(1));
+ }
+ }
+
+ };
+
}
}