|
|
Boost-Commit : |
Subject: [Boost-commit] svn:boost r56901 - sandbox/statistics/importance_sampling/boost/statistics/detail/importance_sampling/weights
From: erwann.rogard_at_[hidden]
Date: 2009-10-16 01:13:07
Author: e_r
Date: 2009-10-16 01:13:06 EDT (Fri, 16 Oct 2009)
New Revision: 56901
URL: http://svn.boost.org/trac/boost/changeset/56901
Log:
a
Added:
sandbox/statistics/importance_sampling/boost/statistics/detail/importance_sampling/weights/effective_sample_size.hpp (contents, props changed)
sandbox/statistics/importance_sampling/boost/statistics/detail/importance_sampling/weights/find_scale_to_finite_sum.hpp (contents, props changed)
sandbox/statistics/importance_sampling/boost/statistics/detail/importance_sampling/weights/maximal_finite_sums.hpp (contents, props changed)
sandbox/statistics/importance_sampling/boost/statistics/detail/importance_sampling/weights/prepare_weights.hpp (contents, props changed)
sandbox/statistics/importance_sampling/boost/statistics/detail/importance_sampling/weights/scale_to_finite_sum.hpp (contents, props changed)
Added: sandbox/statistics/importance_sampling/boost/statistics/detail/importance_sampling/weights/effective_sample_size.hpp
==============================================================================
--- (empty file)
+++ sandbox/statistics/importance_sampling/boost/statistics/detail/importance_sampling/weights/effective_sample_size.hpp 2009-10-16 01:13:06 EDT (Fri, 16 Oct 2009)
@@ -0,0 +1,77 @@
+///////////////////////////////////////////////////////////////////////////////
+// importance_sampling::effective_sample_size.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_STATISTICS_DETAIL_IMPORTANCE_SAMPLING_WEIGHTS_EFFECTIVE_SAMPLE_SIZE_HPP_ER_2009
+#define BOOST_STATISTICS_DETAIL_IMPORTANCE_SAMPLING_WEIGHTS_EFFECTIVE_SAMPLE_SIZE_HPP_ER_2009
+#include <algorithm>
+#include <numeric>
+#include <iterator>
+#include <boost/lambda/lambda.hpp>
+#include <boost/lambda/bind.hpp>
+#include <boost/ref.hpp>
+#include <boost/iterator/iterator_traits.hpp>
+#include <boost/functional/mean_var_accumulator.hpp>
+
+namespace boost{
+namespace statistics{
+namespace detail{
+namespace importance_sampling{
+
+ // Finds the number of iid observations whose standard error equates that
+ // of the importance sample mean. The number is reported in proportion to
+ // the input sample size (m). See sources [3] or [4].
+ //
+ // n = m / (1+Var(w)) where the w's are standardized weights
+ // i.e. w[j]<- uw[j]/mean(uw[j]:j=1,...,m), solves A = B, where
+ // A = V((1/n) sum{y[i],i=1,...,n}), y~p(y) (iid)
+ // B = V(sum{y[j]uw[j],i=1,...,m}/sum{uw[j]:j=1,...,m}), y~q, uw propto p/q
+ template<typename InIt>
+ typename iterator_value<InIt>::type
+ percentage_effective_sample_size(
+ InIt b_w, // un-normalized weights
+ InIt e_w
+ );
+
+ // Implementation //
+
+ template<typename InIt>
+ typename iterator_value<InIt>::type
+ percentage_effective_sample_size(
+ InIt b_w, // un-normalized weights
+ InIt e_w
+ ){
+ typedef typename iterator_value<InIt>::type val_;
+ typedef typename functional::mean_var_accumulator<val_>::type acc_;
+
+ // Var(w/c) = Var(w) / c^2
+
+ static val_ one = static_cast<val_>(1);
+
+ acc_ acc;
+ std::for_each(
+ b_w,
+ e_w,
+ lambda::bind<void>(boost::ref(acc),lambda::_1)
+ );
+ val_ v = static_cast<val_>(
+ boost::accumulators::variance(acc)
+ );
+ val_ c = static_cast<val_>(
+ boost::accumulators::mean(acc)
+ );
+
+ v /= (c*c);
+
+ return one / (one + v);
+ }
+
+}// importance_weights
+}// detail
+}// statistics
+}// boost
+
+#endif
\ No newline at end of file
Added: sandbox/statistics/importance_sampling/boost/statistics/detail/importance_sampling/weights/find_scale_to_finite_sum.hpp
==============================================================================
--- (empty file)
+++ sandbox/statistics/importance_sampling/boost/statistics/detail/importance_sampling/weights/find_scale_to_finite_sum.hpp 2009-10-16 01:13:06 EDT (Fri, 16 Oct 2009)
@@ -0,0 +1,143 @@
+///////////////////////////////////////////////////////////////////////////////
+// importance_sampling::find_scale_to_finite_sum.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_STATISTICS_DETAIL_IMPORTANCE_SAMPLING_WEIGHTS_FIND_SCALE_TO_FINITE_SUM_HPP_ER_2009
+#define BOOST_STATISTICS_DETAIL_IMPORTANCE_SAMPLING_WEIGHTS_FIND_SCALE_TO_FINITE_SUM_HPP_ER_2009
+#include <cmath>
+#include <numeric>
+#include <stdexcept>
+#include <boost/lambda/lambda.hpp>
+#include <boost/range.hpp>
+#include <boost/format.hpp>
+#include <boost/math/special_functions/fpclassify.hpp>
+#include <boost/math/tools/precision.hpp>
+#include <boost/numeric/conversion/bounds.hpp>
+#include <boost/statistics/detail/importance_sampling/weights/maximal_finite_sums.hpp>
+
+namespace boost{
+namespace statistics{
+namespace detail{
+namespace importance_sampling{
+
+ // Finds c such that the sum over sum{*i/c, i in [b,e)} < inf using the
+ // bisection method.
+ //
+ // Warning: The value c is not insensitive to permutations of [b,e), due
+ // to non-associativity in the fp system
+ template<typename InIt>
+ typename iterator_value<InIt>::type
+ find_scale_to_finite_sum(
+ InIt b,InIt e,
+ typename iterator_value<InIt>::type low_init,
+ typename iterator_value<InIt>::type high_init
+ );
+
+ // This version may be faster than that above with low = 1, high = highest
+ template<typename InIt>
+ typename iterator_value<InIt>::type
+ find_scale_to_finite_sum(InIt b,InIt e);
+
+ // Implementation //
+
+ template<typename InIt>
+ typename iterator_value<InIt>::type
+ find_scale_to_finite_sum(
+ InIt b,InIt e,
+ typename iterator_value<InIt>::type low_init,
+ typename iterator_value<InIt>::type high_init
+ ){
+ typedef typename iterator_value<InIt>::type val_;
+ static val_ zero = static_cast<val_>(0);
+ static val_ two = static_cast<val_>(2);
+ static val_ eps = math::tools::epsilon<val_>();
+
+ BOOST_ASSERT(low_init < high_init);
+ BOOST_ASSERT(low_init > zero);
+ BOOST_ASSERT(high_init > zero);
+ val_ low = low_init;
+ val_ high = high_init;
+ val_ mid = (low + high)/two;
+ val_ delta, acc;
+
+ do{
+ delta = high - low;
+ acc =
+ std::accumulate(b,e,zero, lambda::_1 + ( lambda::_2 / mid) );
+ if(boost::math::isinf(acc)){
+ low = mid;
+ }else{
+ high = mid;
+ }
+ mid = (low+high)/two;
+ }while(
+ delta - (high - low)>eps
+ );
+ static
+ const char* str = "%3% = find_scale_to_finite_sum(b,e,%1%,%2%) = inf";
+ if(
+ boost::math::isinf(
+ std::accumulate(b,e,zero, lambda::_1 + ( lambda::_2 / high) )
+ )
+ ){
+ format f(str); f%low_init%high_init%high;
+ throw std::runtime_error(f.str());
+ }
+ return high;
+ }
+
+ template<typename InIt>
+ typename iterator_value<InIt>::type
+ find_scale_to_finite_sum(InIt b,InIt e){
+ typedef typename iterator_value<InIt>::type val_;
+ typedef numeric::bounds<val_> bounds_;
+ typedef std::vector<val_> vec_;
+ static val_ zero = static_cast<val_>(0);
+ static val_ one = static_cast<val_>(1);
+ static val_ two = static_cast<val_>(2);
+ static val_ highest = bounds_::highest();
+ vec_ vec;
+ maximal_finite_sums(b,e,std::back_inserter(vec));
+ val_ low = one;
+ val_ high = highest;
+ val_ mid = find_scale_to_finite_sum(
+ boost::begin(vec),
+ boost::end(vec),
+ low,high
+ );
+ low = mid;
+ high = mid;
+ while(
+ !boost::math::isinf(
+ std::accumulate(b,e,zero, lambda::_1 + ( lambda::_2 / low) )
+ )
+ ){
+ low /= two;
+ if(low<one){
+ low = one;
+ break;
+ }
+ }
+ while(
+ boost::math::isinf(
+ std::accumulate(b,e,zero, lambda::_1 + ( lambda::_2 / high) )
+ )
+ ){
+ high *= two;
+ if(boost::math::isinf(high)){
+ high = highest;
+ break;
+ }
+ }
+ return find_scale_to_finite_sum(b,e,low,high);
+ }
+
+}// importance_weights
+}// detail
+}// statistics
+}// boost
+
+#endif
\ No newline at end of file
Added: sandbox/statistics/importance_sampling/boost/statistics/detail/importance_sampling/weights/maximal_finite_sums.hpp
==============================================================================
--- (empty file)
+++ sandbox/statistics/importance_sampling/boost/statistics/detail/importance_sampling/weights/maximal_finite_sums.hpp 2009-10-16 01:13:06 EDT (Fri, 16 Oct 2009)
@@ -0,0 +1,46 @@
+///////////////////////////////////////////////////////////////////////////////
+// importance_sampling::maximal_finite_sums.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_STATISTICS_DETAIL_IMPORTANCE_SAMPLING_WEIGHTS_MAXIMAL_FINITE_SUMS_HPP_ER_2009
+#define BOOST_STATISTICS_DETAIL_IMPORTANCE_SAMPLING_WEIGHTS_MAXIMAL_FINITE_SUMS_HPP_ER_2009
+#include <boost/iterator/iterator_traits.hpp>
+#include <boost/math/special_functions/fpclassify.hpp>
+
+namespace boost{
+namespace statistics{
+namespace detail{
+namespace importance_sampling{
+
+ // Breaks down the summation of *i, i in [b,e) into maximal finite amounts,
+ // each of which are copied to the output iterator i
+ //
+ // Note: Primarily for use by find_scale_to_finite_sum
+ template<typename InIt,typename OutIt>
+ OutIt maximal_finite_sums(InIt b, InIt e, OutIt i)
+ {
+ typedef typename iterator_value<InIt>::type value_type;
+ value_type sum = static_cast<value_type>(0);
+ while(b!=e){
+ value_type d = *b;
+ if( boost::math::isinf( sum + d ) ){
+ *i = sum;
+ sum = d;
+ ++i;
+ }else{
+ sum += d;
+ }
+ ++b;
+ }
+ return i;
+ };
+
+}// importance_weights
+}// detail
+}// statistics
+}// boost
+
+#endif
\ No newline at end of file
Added: sandbox/statistics/importance_sampling/boost/statistics/detail/importance_sampling/weights/prepare_weights.hpp
==============================================================================
--- (empty file)
+++ sandbox/statistics/importance_sampling/boost/statistics/detail/importance_sampling/weights/prepare_weights.hpp 2009-10-16 01:13:06 EDT (Fri, 16 Oct 2009)
@@ -0,0 +1,160 @@
+///////////////////////////////////////////////////////////////////////////////
+// importance_sampling::prepare_weights.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_STATISTICS_DETAIL_IMPORTANCE_SAMPLING_WEIGHTS_PREPARE_WEIGHTS_HPP_ER_2009
+#define BOOST_STATISTICS_DETAIL_IMPORTANCE_SAMPLING_WEIGHTS_PREPARE_WEIGHTS_HPP_ER_2009
+#include <iterator>
+#include <functional>
+#include <boost/format.hpp>
+#include <boost/lambda/lambda.hpp>
+#include <boost/math/tools/precision.hpp>
+#include <boost/statistics/detail/importance_sampling/weights/apply_exp_offset.hpp>
+#include <boost/statistics/detail/importance_sampling/weights/scale_to_finite_sum.hpp>
+#include <boost/statistics/detail/importance_sampling/weights/effective_sample_size.hpp>
+
+namespace boost{
+namespace statistics{
+namespace detail{
+namespace importance_sampling{
+
+ // Warning: read side effects carefully.
+ template<typename T>
+ class prepare_weights{
+ public:
+ typedef T value_type;
+ typedef std::size_t size_type;
+
+ prepare_weights();
+ prepare_weights(value_type max_log);
+ // Default copy/assign
+
+ // [ Input ]
+ // max_log controls precision hence raising it should decr pc_lt_eps
+ // but also incr risk that cum_sum isinf.
+ value_type max_log;
+
+ // [ Output ]
+ value_type offset; // lw <- lw + offset, max{lw}+offset = max_log
+ value_type scaling_factor;// w <- w/c such that sum{w/c}<inf
+ value_type pc_ess; // pc effective sample size
+ value_type pc_lt_eps; // pc w<eps
+
+ // [ Side effect ]
+ // 1) w <- exp(lw+offset)
+ // 2) if needed, w <- w/c such that sum{w} < inf
+ template<typename ItW>
+ void operator()(
+ ItW b_w, // log( unnormalized weights )
+ ItW e_w
+ );
+
+ public:
+ static value_type zero;
+ static value_type eps;
+ static value_type default_max_log;
+ static const char* header;
+
+ };
+
+ template<typename T>
+ std::ostream& operator<<(std::ostream& out,
+ const prepare_weights<T>& that);
+
+ // Implementation
+
+ template<typename T>
+ std::ostream& operator<<(
+ std::ostream& out,
+ const prepare_weights<T>& that
+ ){
+ out <<
+ (
+ format("(%1%,%2%,%3%,%4%)")
+ % that.offset
+ % that.scaling_factor
+ % that.pc_ess
+ % that.pc_lt_eps
+ ).str();
+ return out;
+ }
+
+ template<typename T>
+ const char* prepare_weights<T>::header
+ = "(offset,scaling_factor,pc_ess,pc_lt_eps)";
+
+ template<typename T>
+ typename prepare_weights<T>::value_type
+ prepare_weights<T>::eps = math::tools::epsilon<value_type>();
+
+ template<typename T>
+ typename prepare_weights<T>::value_type
+ prepare_weights<T>::default_max_log = static_cast<value_type>(0);
+
+ template<typename T>
+ typename prepare_weights<T>::value_type
+ prepare_weights<T>::zero = static_cast<value_type>(0);
+
+ template<typename T>
+ prepare_weights<T>::prepare_weights()
+ :max_log(default_max_log),
+ offset(zero),scaling_factor(zero),pc_ess(zero),pc_lt_eps(zero){}
+
+ template<typename T>
+ prepare_weights<T>::prepare_weights(value_type ml)
+ :max_log(ml),
+ offset(zero),scaling_factor(zero),pc_ess(zero),pc_lt_eps(zero){}
+
+ template<typename T>
+ template<typename ItW>
+ void
+ prepare_weights<T>::operator()(
+ ItW b_w,
+ ItW e_w
+ ){
+ offset = apply_exp_offset(
+ b_w,
+ e_w,
+ max_log
+ );
+
+ // if max_log is small enough (which costs precision), this does not
+ // nothing i.e. scaling_factor = 1
+ scaling_factor = scale_to_finite_sum(
+ b_w,
+ e_w
+ );
+
+ ItW i_lt_eps = std::lower_bound(
+ b_w,
+ e_w,
+ eps,
+ ( boost::lambda::_1 >= boost::lambda::_2 )
+ );
+
+ value_type n_gt_eps
+ = static_cast<value_type>( std::distance(b_w,i_lt_eps) );
+ value_type n_lt_eps
+ = static_cast<value_type>( std::distance(i_lt_eps,e_w) );
+
+
+ // Increasing max_log should decrease this number
+ pc_lt_eps = n_lt_eps / ( n_lt_eps + n_gt_eps ) ;
+
+ // Beware that pc_lt_eps >0 may distort ess
+ pc_ess = percentage_effective_sample_size(
+ b_w,
+ e_w
+ );
+
+ }
+
+}// importance_weights
+}// detail
+}// statistics
+}// boost
+
+#endif
\ No newline at end of file
Added: sandbox/statistics/importance_sampling/boost/statistics/detail/importance_sampling/weights/scale_to_finite_sum.hpp
==============================================================================
--- (empty file)
+++ sandbox/statistics/importance_sampling/boost/statistics/detail/importance_sampling/weights/scale_to_finite_sum.hpp 2009-10-16 01:13:06 EDT (Fri, 16 Oct 2009)
@@ -0,0 +1,76 @@
+///////////////////////////////////////////////////////////////////////////////
+// importance_sampling::scale_to_finite_sum.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_STATISTICS_DETAIL_IMPORTANCE_SAMPLING_WEIGHTS_SCALE_TO_FINITE_SUM_HPP_ER_2009
+#define BOOST_STATISTICS_DETAIL_IMPORTANCE_SAMPLING_WEIGHTS_SCALE_TO_FINITE_SUM_HPP_ER_2009
+#include <numeric>
+#include <boost/lambda/lambda.hpp>
+#include <boost/range.hpp>
+#include <boost/statistics/detail/importance_sampling/weights/find_scale_to_finite_sum.hpp>
+
+namespace boost{
+namespace statistics{
+namespace detail{
+namespace importance_sampling{
+
+ // Scales each element [b,e) by the smallest factor, c, such that the sum
+ // is finite.
+ //
+ // [ Warning ] c is not insensitive to permutations of [b,e), due
+ // to non-associativity in the fp system
+ template<typename InIt>
+ typename iterator_value<InIt>::type
+ scale_to_finite_sum(InIt b,InIt e,
+ typename iterator_value<InIt>::type low,
+ typename iterator_value<InIt>::type high
+ );
+
+ template<typename InIt>
+ typename iterator_value<InIt>::type
+ scale_to_finite_sum(InIt b,InIt e);
+
+ // Implementation //
+
+ template<typename InIt>
+ typename iterator_value<InIt>::type
+ scale_to_finite_sum(InIt b,InIt e,
+ typename iterator_value<InIt>::type low,
+ typename iterator_value<InIt>::type high
+ ){
+ typedef typename iterator_value<InIt>::type val_;
+ val_ c = find_scale_finite_sum(
+ b,e,low,high
+ );
+ std::transform(
+ b,
+ e,
+ b,
+ boost::lambda::_1 / c
+ );
+ return c;
+ }
+
+ template<typename InIt>
+ typename iterator_value<InIt>::type
+ scale_to_finite_sum(InIt b,InIt e){
+ typedef typename iterator_value<InIt>::type val_;
+ val_ c = find_scale_to_finite_sum(b,e);
+ std::transform(
+ b,
+ e,
+ b,
+ boost::lambda::_1 / c
+ );
+ return c;
+ }
+
+}// importance_weights
+}// detail
+}// statistics
+}// boost
+
+#endif
\ No newline at end of file
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