Subject: [Boost-commit] svn:boost r55847 - in sandbox/statistics/kernel/boost: . statistics statistics/kernel statistics/kernel/bandwidth statistics/kernel/functional statistics/kernel/functional/detail statistics/kernel/joint statistics/kernel/meta statistics/kernel/scalar
From: erwann.rogard_at_[hidden]
Date: 2009-08-28 18:41:49
Author: e_r
Date: 2009-08-28 18:41:48 EDT (Fri, 28 Aug 2009)
New Revision: 55847
URL: http://svn.boost.org/trac/boost/changeset/55847
Log:
replacing kernel/boost
Added:
sandbox/statistics/kernel/boost/
sandbox/statistics/kernel/boost/statistics/
sandbox/statistics/kernel/boost/statistics/kernel/
sandbox/statistics/kernel/boost/statistics/kernel/bandwidth/
sandbox/statistics/kernel/boost/statistics/kernel/bandwidth/include.hpp (contents, props changed)
sandbox/statistics/kernel/boost/statistics/kernel/bandwidth/normal_distribution.hpp (contents, props changed)
sandbox/statistics/kernel/boost/statistics/kernel/functional/
sandbox/statistics/kernel/boost/statistics/kernel/functional/benchmark_nw.hpp (contents, props changed)
sandbox/statistics/kernel/boost/statistics/kernel/functional/benchmark_rp.hpp (contents, props changed)
sandbox/statistics/kernel/boost/statistics/kernel/functional/detail/
sandbox/statistics/kernel/boost/statistics/kernel/functional/detail/mean_accumulator.hpp (contents, props changed)
sandbox/statistics/kernel/boost/statistics/kernel/functional/detail/return_if.hpp (contents, props changed)
sandbox/statistics/kernel/boost/statistics/kernel/functional/estimator.hpp (contents, props changed)
sandbox/statistics/kernel/boost/statistics/kernel/functional/include.hpp (contents, props changed)
sandbox/statistics/kernel/boost/statistics/kernel/functional/nw_visitor.hpp (contents, props changed)
sandbox/statistics/kernel/boost/statistics/kernel/functional/nw_visitor_tuple.hpp (contents, props changed)
sandbox/statistics/kernel/boost/statistics/kernel/functional/rp_visitor.hpp (contents, props changed)
sandbox/statistics/kernel/boost/statistics/kernel/include.hpp (contents, props changed)
sandbox/statistics/kernel/boost/statistics/kernel/joint/
sandbox/statistics/kernel/boost/statistics/kernel/joint/crtp.hpp (contents, props changed)
sandbox/statistics/kernel/boost/statistics/kernel/joint/include.hpp (contents, props changed)
sandbox/statistics/kernel/boost/statistics/kernel/joint/kernel_mono.hpp (contents, props changed)
sandbox/statistics/kernel/boost/statistics/kernel/meta/
sandbox/statistics/kernel/boost/statistics/kernel/scalar/
sandbox/statistics/kernel/boost/statistics/kernel/scalar/crtp.hpp (contents, props changed)
sandbox/statistics/kernel/boost/statistics/kernel/scalar/gaussian.hpp (contents, props changed)
sandbox/statistics/kernel/boost/statistics/kernel/scalar/include.hpp (contents, props changed)
Added: sandbox/statistics/kernel/boost/statistics/kernel/bandwidth/include.hpp
==============================================================================
--- (empty file)
+++ sandbox/statistics/kernel/boost/statistics/kernel/bandwidth/include.hpp 2009-08-28 18:41:48 EDT (Fri, 28 Aug 2009)
@@ -0,0 +1,13 @@
+///////////////////////////////////////////////////////////////////////////////
+// statistics::kernel::bandwidth::include.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_KERNEL_BANDWIDTH_INCLUDE_H_ER_2009
+#define BOOST_STATISTICS_KERNEL_BANDWIDTH_INCLUDE_H_ER_2009
+
+#include <boost/statistics/kernel/bandwidth/normal_distribution.hpp>
+
+#endif
Added: sandbox/statistics/kernel/boost/statistics/kernel/bandwidth/normal_distribution.hpp
==============================================================================
--- (empty file)
+++ sandbox/statistics/kernel/boost/statistics/kernel/bandwidth/normal_distribution.hpp 2009-08-28 18:41:48 EDT (Fri, 28 Aug 2009)
@@ -0,0 +1,49 @@
+//////////////////////////////////////////////////////////////////////////////
+// boost::statistics::kernel::bandwidth::normal_distribution.hpp //
+// //
+// (C) Copyright 2009 Erwann Rogard //
+// Use, modification and distribution are subject to 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_KERNEL_BANDWIDTH_NORMAL_DISRTIBUTION_HPP_ER_2009
+#define BOOST_STATISTICS_KERNEL_BANDWIDTH_NORMAL_DISRTIBUTION_HPP_ER_2009
+#include <cmath>
+
+namespace boost{
+namespace statistics{
+namespace kernel{
+
+
+// Finds the optimal bandwidth for density estimation (Rosenblatt-Parzen)
+// assuming
+// 1) The data is a normal vector each coordinate having stddev sigma
+// 2) The kernel is Gaussian.
+
+template<unsigned M, typename T>
+struct normal_distribution{
+ typedef T value_type;
+ static value_type rp_bandwidth(std::size_t n); // Assumes sigma = 1
+ static value_type rp_bandwidth(value_type sigma,unsigned n);
+ // For M == 1, bandwdith = sigma (3n/4)^(-1/5) = 1.06 sigma n^(-1/5)
+};
+template<unsigned M, typename T>
+typename normal_distribution<M,T>::value_type
+normal_distribution<M,T>::rp_bandwidth(std::size_t n){
+ static const value_type e = -static_cast<T>(1)/static_cast<T>(M+4);
+ static const value_type r = static_cast<T>(2+M)/ static_cast<T>(4);
+ value_type a = static_cast<T>(n) * r;
+ return std::pow( a, e );
+}
+
+template<unsigned M, typename T>
+typename normal_distribution<M,T>::value_type
+normal_distribution<M,T>::rp_bandwidth(value_type sigma,unsigned n){
+ return sigma * rp_bandwidth(n);
+}
+
+}// kernel
+}// statistics
+}// boost
+
+#endif
Added: sandbox/statistics/kernel/boost/statistics/kernel/functional/benchmark_nw.hpp
==============================================================================
--- (empty file)
+++ sandbox/statistics/kernel/boost/statistics/kernel/functional/benchmark_nw.hpp 2009-08-28 18:41:48 EDT (Fri, 28 Aug 2009)
@@ -0,0 +1,252 @@
+///////////////////////////////////////////////////////////////////////////////
+// statistics::kernel::functional::benchmark_nw.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_KERNEL_FUNCTIONAL_BENCHMARK_NW_H_ER_2009
+#define BOOST_STATISTICS_KERNEL_FUNCTIONAL_BENCHMARK_NW_H_ER_2009
+#include <cmath>
+#include <ostream>
+#include <string>
+#include <boost/format.hpp>
+#include <boost/timer.hpp>
+#include <boost/foreach.hpp>
+#include <boost/call_traits.hpp>
+#include <boost/mpl/not.hpp>
+#include <boost/range.hpp>
+#include <boost/type_traits.hpp>
+#include <boost/binary_op/data/tuple_range.hpp>
+#include <boost/statistics/kernel/functional/detail/mean_accumulator.hpp>
+#include <boost/statistics/kernel/functional/nw_visitor_tuple.hpp>
+#include <boost/statistics/kernel/functional/estimator.hpp>
+
+namespace boost{
+namespace statistics{
+namespace kernel{
+
+// Computes the error in estimating the conditional mean with a
+// particular kernel. see benchmark_rp.hpp
+template<typename Sx,typename Sp,typename Sy,typename K, typename Ae =
+ typename detail::mean_accumulator<typename K::result_type>::type
+>
+class benchmark_nw{
+public:
+ typedef typename K::result_type value_type;
+ typedef typename statistics::kernel::detail::mean_accumulator<value_type>::type acc_;
+ private:
+ typedef is_reference<Sx> is_ref_sx_;
+ typedef is_reference<Sp> is_ref_sp_;
+ typedef is_reference<Sp> is_ref_sy_;
+ typedef typename remove_reference<Sx>::type const_sx_;
+ typedef typename remove_reference<Sp>::type const_sp_;
+ typedef typename remove_reference<Sy>::type const_sy_;
+ public:
+ // Construct
+ benchmark_nw();
+ benchmark_nw(
+ typename call_traits<Sx>::param_type,
+ typename call_traits<Sp>::param_type,
+ typename call_traits<Sy>::param_type,
+ K k
+ );
+ benchmark_nw(
+ typename call_traits<Sx>::param_type,
+ typename call_traits<Sp>::param_type,
+ typename call_traits<Sy>::param_type
+ );
+ benchmark_nw(const benchmark_nw&);
+ benchmark_nw& operator=(const benchmark_nw&);
+
+ void set_kernel(K k);
+
+ // Update
+ template<typename Sx1,typename Sy1>
+ void operator()(const Sx1& in_x,const Sy1& in_y);
+
+ // Access
+ const K& kernel()const;
+ const Ae& error_accumulator_p()const;
+ const Ae& error_accumulator_y()const;
+ value_type error_p()const;
+ value_type error_y()const;
+ unsigned n()const; // in-sample size
+ unsigned m()const; // out-sample size
+ value_type time()const;
+
+ static std::string header;
+
+ private:
+ K k_;
+ typename call_traits<Sx>::value_type sx_;
+ typename call_traits<Sp>::value_type sp_;
+ typename call_traits<Sy>::value_type sy_;
+ Ae error_acc_p_;
+ Ae error_acc_y_;
+ value_type time_;
+ unsigned n_;
+};
+
+template<typename Sx,typename Sp, typename Sy,typename K,typename Ae>
+std::string benchmark_nw<Sx,Sp,Sy,K,Ae>::header = "(n, avg_t,e_p,e_y)";
+
+template<typename Sx,typename Sp, typename Sy,typename K,typename Ae>
+std::ostream&
+operator<<(std::ostream& out,const benchmark_nw<Sx,Sp,Sy,K,Ae>& b){
+ static const char* str = "(%1%, %2%, %3%, %4%)";
+ typedef benchmark_nw<Sx,Sp,Sy,K,Ae> bench_;
+ typedef typename bench_::value_type value_type;
+ value_type m_ = static_cast<value_type>(b.m());
+ value_type n_ = static_cast<value_type>(b.n());
+ value_type avg_t = (b.time()) / (m_); //Should increase linearly with n
+ format f(str);
+ f%(b.n())%(avg_t)%(b.error_p())%(b.error_y());
+ out << f.str();
+ return out;
+}
+
+// Construction
+template<typename Sx,typename Sp, typename Sy,typename K,typename Ae>
+benchmark_nw<Sx,Sp,Sy,K,Ae>::benchmark_nw(){
+ BOOST_MPL_ASSERT((mpl::not_<is_ref_sx_>));
+ BOOST_MPL_ASSERT((mpl::not_<is_ref_sp_>));
+ BOOST_MPL_ASSERT((mpl::not_<is_ref_sy_>));
+}
+
+template<typename Sx,typename Sp, typename Sy,typename K,typename Ae>
+benchmark_nw<Sx,Sp,Sy,K,Ae>::benchmark_nw(
+ typename call_traits<Sx>::param_type sx,
+ typename call_traits<Sp>::param_type sp,
+ typename call_traits<Sy>::param_type sy,
+ K k
+):k_(k),sx_(sx),sp_(sp),sy_(sy),
+error_acc_p_(),error_acc_y_(),time_(static_cast<value_type>(0)){}
+
+template<typename Sx,typename Sp, typename Sy,typename K,typename Ae>
+benchmark_nw<Sx,Sp,Sy,K,Ae>::benchmark_nw(
+ typename call_traits<Sx>::param_type sx,
+ typename call_traits<Sp>::param_type sp,
+ typename call_traits<Sy>::param_type sy
+):k_(K()),sx_(sx),sp_(sp),sy_(sy),
+error_acc_p_(),error_acc_y_(),time_(static_cast<value_type>(0)){}
+
+template<typename Sx,typename Sp, typename Sy,typename K,typename Ae>
+benchmark_nw<Sx,Sp,Sy,K,Ae>::benchmark_nw(const benchmark_nw& that)
+:k_(that.k_),sx_(that.sx_),sp_(that.sp_),sy_(that.sy_),
+error_acc_p_(that.error_acc_p_),error_acc_y_(that.error_acc_y_),
+time_(that.time_){}
+
+template<typename Sx,typename Sp, typename Sy,typename K,typename Ae>
+benchmark_nw<Sx,Sp,Sy,K,Ae>&
+benchmark_nw<Sx,Sp,Sy,K,Ae>::operator=(const benchmark_nw& that){
+ if(&that!=this){
+ BOOST_MPL_ASSERT((mpl::not_<is_ref_sx_>));
+ BOOST_MPL_ASSERT((mpl::not_<is_ref_sp_>));
+ BOOST_MPL_ASSERT((mpl::not_<is_ref_sy_>));
+ k_ = (that.k_);
+ sx_ = (that.sx_);
+ sp_ = (that.sp_);
+ sy_ = (that.sy_);
+ error_acc_p_ = (that.error_acc_p_);
+ error_acc_p_ = (that.error_acc_y_);
+ time_ = (that.time_);
+ }
+ return (*this);
+}
+
+template<typename Sx,typename Sp, typename Sy,typename K,typename Ae>
+void benchmark_nw<Sx,Sp,Sy,K,Ae>::set_kernel(K k){
+ k_ = k;
+}
+
+// Update
+template<typename Sx,typename Sp, typename Sy,typename K,typename Ae>
+template<typename Sx1,typename Sy1>
+void benchmark_nw<Sx,Sp,Sy,K,Ae>::operator()(const Sx1& sx1,const Sy1& sy1){
+ typedef const Sx1& call_x1_;
+ typedef const Sy1& call_y1_;
+ typedef binary_op::tuple_range<call_x1_,call_y1_> factory_;
+ typedef typename factory_::type range_tuple_;
+ typedef estimator<range_tuple_,nw_visitor_tuple,K,acc_> estimator_type;
+
+ timer t;
+ range_tuple_ range_tuple = factory_::make(sx1,sy1);
+ estimator_type estimator(this->kernel());
+ train(estimator,range_tuple);
+
+ typedef typename range_iterator<const_sx_>::type iter_sx_;
+ typedef typename range_iterator<const_sp_>::type iter_sp_;
+ typedef typename range_iterator<const_sy_>::type iter_sy_;
+ typedef typename range_iterator<range_tuple_>::type iter_data_;
+
+ iter_sx_ i_x = boost::begin(sx_);
+ iter_sx_ e_x = boost::end(sx_);
+ iter_sp_ i_p = boost::begin(sp_);
+ iter_sy_ i_y = boost::begin(sy_);
+ value_type a, b;
+ while(i_x != e_x){
+ BOOST_AUTO(obj,estimator(*i_x));
+ a = obj.rp_estimate();
+ b = (*i_p);
+ error_acc_p_(fabs(a-b));
+ a = obj.nw_estimate();
+ b = (*i_y);
+ error_acc_y_(fabs(a-b));
+ ++i_x;
+ ++i_p;
+ ++i_y;
+ }
+ n_ = boost::size(range_tuple);
+ time_ = t.elapsed();
+}
+
+// Access
+template<typename Sx,typename Sp, typename Sy,typename K,typename Ae>
+const K& benchmark_nw<Sx,Sp,Sy,K,Ae>::kernel()const{
+ return (this->k_);
+}
+
+template<typename Sx,typename Sp, typename Sy,typename K,typename Ae>
+const Ae& benchmark_nw<Sx,Sp,Sy,K,Ae>::error_accumulator_p()const{
+ return (this->error_acc_p_);
+}
+
+template<typename Sx,typename Sp, typename Sy,typename K,typename Ae>
+const Ae& benchmark_nw<Sx,Sp,Sy,K,Ae>::error_accumulator_y()const{
+ return (this->error_acc_y_);
+}
+
+template<typename Sx,typename Sp, typename Sy,typename K,typename Ae>
+typename benchmark_nw<Sx,Sp,Sy,K,Ae>::value_type
+benchmark_nw<Sx,Sp,Sy,K,Ae>::error_p()const{
+ return accumulators::mean(this->error_acc_p_);
+}
+
+template<typename Sx,typename Sp, typename Sy,typename K,typename Ae>
+typename benchmark_nw<Sx,Sp,Sy,K,Ae>::value_type
+benchmark_nw<Sx,Sp,Sy,K,Ae>::error_y()const{
+ return accumulators::mean(this->error_acc_y_);
+}
+
+template<typename Sx,typename Sp, typename Sy,typename K,typename Ae>
+unsigned benchmark_nw<Sx,Sp,Sy,K,Ae>::n()const{
+ return (this->n_);
+}
+
+template<typename Sx,typename Sp, typename Sy,typename K,typename Ae>
+unsigned benchmark_nw<Sx,Sp,Sy,K,Ae>::m()const{
+ return accumulators::count(this->error_acc_p_);
+}
+
+template<typename Sx,typename Sp, typename Sy,typename K,typename Ae>
+typename benchmark_nw<Sx,Sp,Sy,K,Ae>::value_type
+benchmark_nw<Sx,Sp,Sy,K,Ae>::time()const{
+ return (this->time_);
+}
+
+}// kernel
+}// statistics
+}// boost
+
+#endif
Added: sandbox/statistics/kernel/boost/statistics/kernel/functional/benchmark_rp.hpp
==============================================================================
--- (empty file)
+++ sandbox/statistics/kernel/boost/statistics/kernel/functional/benchmark_rp.hpp 2009-08-28 18:41:48 EDT (Fri, 28 Aug 2009)
@@ -0,0 +1,220 @@
+///////////////////////////////////////////////////////////////////////////////
+// statistics::kernel::functional::benchmark_rp.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_KERNEL_FUNCTIONAL_BENCHMARK_RP_H_ER_2009
+#define BOOST_STATISTICS_KERNEL_FUNCTIONAL_BENCHMARK_RP_H_ER_2009
+#include <cmath>
+#include <ostream>
+#include <string>
+#include <boost/format.hpp>
+#include <boost/timer.hpp>
+#include <boost/foreach.hpp>
+#include <boost/call_traits.hpp>
+#include <boost/mpl/not.hpp>
+#include <boost/range.hpp>
+#include <boost/type_traits.hpp>
+#include <boost/statistics/kernel/functional/detail/mean_accumulator.hpp>
+#include <boost/statistics/kernel/functional/rp_visitor.hpp>
+#include <boost/statistics/kernel/functional/estimator.hpp>
+
+namespace boost{
+namespace statistics{
+namespace kernel{
+
+// Computes the error in estimating the density with a particular kernel.
+//
+// Ae specifies the accumulator used to compute the error
+//
+// Let sx of type Sx denote a sequence of values and sp of type Sx
+// their correponding density under a given probability distribution. Let
+// sx1 of type Sx1 denote a (training) sample from the same distribution.
+//
+// Usage:
+// typedef benchmark_rp<Sx,Sp,K,Ae> bench_;
+// bench_ bench(sx,sp);
+// bench.set_kernel(bandwidth);
+// bench(sx1);
+// bench.error();
+template<typename Sx,typename Sp,typename K, typename Ae =
+ typename statistics::kernel::detail::mean_accumulator<typename K::result_type>::type
+>
+class benchmark_rp{
+ public:
+ typedef typename K::result_type value_type;
+ typedef typename statistics::kernel::detail::mean_accumulator<value_type>::type acc_;
+ private:
+ typedef is_reference<Sx> is_ref_sx_;
+ typedef is_reference<Sp> is_ref_sp_;
+ typedef typename remove_reference<Sp>::type const_sp_;
+ public:
+ // Construct
+ benchmark_rp();
+ benchmark_rp(
+ typename call_traits<Sx>::param_type,
+ typename call_traits<Sp>::param_type,
+ K k
+ );
+ benchmark_rp(
+ typename call_traits<Sx>::param_type,
+ typename call_traits<Sp>::param_type
+ );
+ benchmark_rp(const benchmark_rp&);
+ benchmark_rp& operator=(const benchmark_rp&);
+
+ void set_kernel(K k);
+
+ // Update
+ template<typename Sx1> void operator()(const Sx1& in_sample);
+
+ // Access
+ const K& kernel()const;
+ const Ae& error_accumulator()const;
+ value_type error()const;
+ unsigned n()const; // in-sample size
+ unsigned m()const; // out-sample size
+ value_type time()const;
+
+ static std::string header;
+
+ private:
+ K k_;
+ typename call_traits<Sx>::value_type sx_;
+ typename call_traits<Sp>::value_type sp_;
+ Ae error_acc_;
+ value_type time_;
+ unsigned n_;
+};
+
+template<typename Sx,typename Sp, typename K,typename Ae>
+std::string benchmark_rp<Sx,Sp,K,Ae>::header = "(n, avg_t,e)";
+
+template<typename Sx,typename Sp, typename K,typename Ae>
+std::ostream& operator<<(std::ostream& out,const benchmark_rp<Sx,Sp,K,Ae>& b){
+ static const char* str = "(%1%, %2%, %3%)";
+ typedef benchmark_rp<Sx,Sp,K,Ae> bench_;
+ typedef typename bench_::value_type value_type;
+ value_type m_ = static_cast<value_type>(b.m());
+ value_type n_ = static_cast<value_type>(b.n());
+ value_type avg_t = (b.time()) / (m_); //Should increase linearly with n
+ format f(str);
+ f%(b.n())%(avg_t)%(b.error());
+ out << f.str();
+ return out;
+}
+
+
+// Construction
+template<typename Sx,typename Sp, typename K,typename Ae>
+benchmark_rp<Sx,Sp,K,Ae>::benchmark_rp(){
+ BOOST_MPL_ASSERT((mpl::not_<is_ref_sx_>));
+ BOOST_MPL_ASSERT((mpl::not_<is_ref_sp_>));
+}
+
+template<typename Sx,typename Sp, typename K,typename Ae>
+benchmark_rp<Sx,Sp,K,Ae>::benchmark_rp(
+ typename call_traits<Sx>::param_type sx,
+ typename call_traits<Sp>::param_type sp,
+ K k
+):k_(k),sx_(sx),sp_(sp),
+error_acc_(),time_(static_cast<value_type>(0)){
+ BOOST_ASSERT(size(sx) == size(sp));
+}
+
+template<typename Sx,typename Sp, typename K,typename Ae>
+benchmark_rp<Sx,Sp,K,Ae>::benchmark_rp(
+ typename call_traits<Sx>::param_type sx,
+ typename call_traits<Sp>::param_type sp
+):k_(K()),sx_(sx),sp_(sp),
+error_acc_(),time_(static_cast<value_type>(0)){
+ BOOST_ASSERT(size(sx) == size(sp));
+}
+
+template<typename Sx,typename Sp, typename K,typename Ae>
+benchmark_rp<Sx,Sp,K,Ae>::benchmark_rp(const benchmark_rp& that)
+:k_(that.k_),sx_(that.sx_),sp_(that.sp_),
+error_acc_(that.error_acc_),time_(that.time_){}
+
+template<typename Sx,typename Sp, typename K,typename Ae>
+benchmark_rp<Sx,Sp,K,Ae>&
+benchmark_rp<Sx,Sp,K,Ae>::operator=(const benchmark_rp& that){
+ if(&that!=this){
+ BOOST_MPL_ASSERT((mpl::not_<is_ref_sx_>));
+ BOOST_MPL_ASSERT((mpl::not_<is_ref_sp_>));
+ k_ = (that.k_);
+ sx_ = (that.sx_);
+ sp_ = (that.sp_);
+ error_acc_ = (that.error_acc_);
+ time_ = (that.time_);
+ }
+ return (*this);
+}
+
+template<typename Sx,typename Sp, typename K,typename Ae>
+void benchmark_rp<Sx,Sp,K,Ae>::set_kernel(K k){
+ k_ = k;
+}
+
+// Update
+template<typename Sx,typename Sp, typename K,typename Ae>
+template<typename Sx1>
+void benchmark_rp<Sx,Sp,K,Ae>::operator()(const Sx1& data){
+ typedef const Sx1& call_;
+ typedef estimator<call_,rp_visitor,K,acc_> estimator_type;
+ timer t;
+ estimator_type estimator(data,this->kernel());
+ typedef typename range_value<Sx1>::type unit_;
+ typedef typename range_iterator<const_sp_>::type iter_y_;
+
+ iter_y_ i = boost::begin(sp_);
+ BOOST_FOREACH(const unit_& u, sx_){
+ value_type a = estimator(u).estimate();
+ value_type b = (*i);
+ error_acc_(fabs(a-b));
+ ++i;
+ }
+ n_ = boost::size(data);
+ time_ = t.elapsed();
+}
+
+// Access
+template<typename Sx,typename Sp, typename K,typename Ae>
+const K& benchmark_rp<Sx,Sp,K,Ae>::kernel()const{
+ return (this->k_);
+}
+
+template<typename Sx,typename Sp, typename K,typename Ae>
+const Ae& benchmark_rp<Sx,Sp,K,Ae>::error_accumulator()const{
+ return (this->error_acc_);
+}
+
+template<typename Sx,typename Sp, typename K,typename Ae>
+typename benchmark_rp<Sx,Sp,K,Ae>::value_type
+benchmark_rp<Sx,Sp,K,Ae>::error()const{
+ return accumulators::mean(this->error_acc_);
+}
+
+template<typename Sx,typename Sp, typename K,typename Ae>
+unsigned benchmark_rp<Sx,Sp,K,Ae>::n()const{
+ return (this->n_);
+}
+
+template<typename Sx,typename Sp, typename K,typename Ae>
+unsigned benchmark_rp<Sx,Sp,K,Ae>::m()const{
+ return accumulators::count(this->error_acc_);
+}
+
+template<typename Sx,typename Sp, typename K,typename Ae>
+typename benchmark_rp<Sx,Sp,K,Ae>::value_type
+benchmark_rp<Sx,Sp,K,Ae>::time()const{
+ return (this->time_);
+}
+
+}// kernel
+}// statistics
+}// boost
+
+#endif
Added: sandbox/statistics/kernel/boost/statistics/kernel/functional/detail/mean_accumulator.hpp
==============================================================================
--- (empty file)
+++ sandbox/statistics/kernel/boost/statistics/kernel/functional/detail/mean_accumulator.hpp 2009-08-28 18:41:48 EDT (Fri, 28 Aug 2009)
@@ -0,0 +1,33 @@
+///////////////////////////////////////////////////////////////////////////////
+// statistics::kernel::functional::detail::mean_accumulator.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_KERNEL_FUNCTIONAL_DETAIL_MEAN_ACCUMULATOR_H_ER_2009
+#define BOOST_STATISTICS_KERNEL_FUNCTIONAL_DETAIL_MEAN_ACCUMULATOR_H_ER_2009
+#include <boost/accumulators/accumulators.hpp>
+#include <boost/accumulators/statistics/stats.hpp>
+#include <boost/accumulators/statistics/mean.hpp>
+
+namespace boost{
+namespace statistics{
+namespace kernel{
+namespace detail{
+
+ // The Rosenblatt-Parzen estimator is a density estimator
+ template<typename T>
+ struct mean_accumulator{
+ typedef accumulators::stats<
+ accumulators::tag::mean
+ > stat_;
+ typedef accumulators::accumulator_set<T,stat_> type;
+ };
+
+}// detail
+}// kernel
+}// statistics
+}// boost
+
+#endif
\ No newline at end of file
Added: sandbox/statistics/kernel/boost/statistics/kernel/functional/detail/return_if.hpp
==============================================================================
--- (empty file)
+++ sandbox/statistics/kernel/boost/statistics/kernel/functional/detail/return_if.hpp 2009-08-28 18:41:48 EDT (Fri, 28 Aug 2009)
@@ -0,0 +1,37 @@
+///////////////////////////////////////////////////////////////////////////////
+// statistics::kernel::functional::detail::return_if.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_KERNEL_FUNCTIONAL_DETAIL_RETURN_IF_H_ER_2009
+#define BOOST_STATISTICS_KERNEL_FUNCTIONAL_DETAIL_RETURN_IF_H_ER_2009
+#include <boost/utility/enable_if.hpp>
+#include <boost/type_traits/is_scalar.hpp>
+
+namespace boost{
+namespace statistics{
+namespace kernel{
+namespace detail{
+
+// TODO no longer needed
+
+template<typename X1,typename R>
+struct return_if_scalar : boost::enable_if<
+ boost::is_scalar<X1>,
+ R
+>{};
+
+template<typename X1,typename R>
+struct return_if_multi : boost::disable_if<
+ boost::is_scalar<X1>,
+ R
+>{};
+
+}
+}
+}
+}
+
+#endif
\ No newline at end of file
Added: sandbox/statistics/kernel/boost/statistics/kernel/functional/estimator.hpp
==============================================================================
--- (empty file)
+++ sandbox/statistics/kernel/boost/statistics/kernel/functional/estimator.hpp 2009-08-28 18:41:48 EDT (Fri, 28 Aug 2009)
@@ -0,0 +1,229 @@
+///////////////////////////////////////////////////////////////////////////////
+// statistics::kernel::functional::rp_visitor.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_KERNEL_FUNCTIONAL_DENSITY_VISITOR_HPP_ER_2009
+#define BOOST_STATISTICS_KERNEL_FUNCTIONAL_DENSITY_VISITOR_HPP_ER_2009
+#include <algorithm>
+#include <boost/call_traits.hpp>
+#include <boost/call_traits.hpp>
+#include <boost/type_traits/is_reference.hpp>
+#include <boost/mpl/nested_type.hpp>
+#include <boost/statistics/estimator_concept/meta/result_of_estimate.hpp>
+#include <boost/statistics/kernel/functional/detail/mean_accumulator.hpp>
+
+namespace boost{
+namespace statistics{
+namespace kernel{
+
+ // R: range of training_data
+ // V: visitor
+ // K: kernel
+ // A: accumulator
+ //
+ // V R
+ // rp_visitor {x[i]:i=1,...,n}
+ // nw_visitor_tuple {(x,y)[i]:i=1,...,n}
+ //
+ // See statistics/estimator_concept
+ // Models TrainableEstimator
+ // If V == nw_visitor_tuple, models TrainableRegressionEstimator
+ template<
+ typename R,
+ template<typename,typename,typename> class V,
+ typename K,
+ typename A = typename
+ detail::mean_accumulator<typename K::result_type>::type
+ >
+ class estimator{
+
+ public:
+ typedef R training_data_type;
+ typedef K kernel_type;
+ typedef A accumulator_type;
+
+ template<typename X>
+ struct result{
+ typedef V<K,const X&,A> arg_;
+ typedef typename mpl::nested_type<arg_>::type type;
+ };
+
+ // Constructor
+ estimator();
+ estimator(K k);
+ estimator(K k,const A&);
+ estimator(const estimator&);
+ estimator& operator=(const estimator&);
+
+ // Evaluate
+ template<typename X>
+ typename result<X>::type operator()(const X& x)const;
+
+ // Access
+ const R& train_data()const;
+ const K& kernel()const;
+ const A& accumulator()const;
+
+ void set_train_data(const R& train_data);
+
+ private:
+ R train_data_;
+ K k_;
+ A a_; // Initialized accumulator
+ };
+
+}//kernel
+
+ template<
+ typename R,
+ template<typename,typename,typename> class V,
+ typename K,
+ typename A
+ >
+ void train(kernel::estimator<R,V,K,A>& e, const R& train_data){
+ e.set_train_data(train_data);
+ }
+
+ template<
+ typename R,
+ template<typename,typename,typename> class V,
+ typename K, typename A, typename X
+ >
+ struct result_of_estimate< kernel::estimator<R,V,K,A>,X > {
+ typedef kernel::estimator<R,V,K,A> est_;
+ typedef typename est_::template result<X> meta_vis_;
+ typedef typename meta_vis_::type vis_;
+ typedef typename vis_::result_type type;
+ };
+
+ template<
+ typename R,
+ template<typename,typename,typename> class V,
+ typename K, typename A, typename X
+ >
+ typename result_of_estimate< kernel::estimator<R,V,K,A>, X >::type
+ estimate(const kernel::estimator<R,V,K,A>& e, const X& x){
+ return e(x).estimate();
+ }
+
+// Implementation //
+namespace kernel{
+
+// Constructor
+
+template<
+ typename R,
+ template<typename,typename,typename> class V,
+ typename K,
+ typename A
+>
+estimator<R,V,K,A>::estimator():train_data_(),k_(),a_(){}
+
+template<
+ typename R,
+ template<typename,typename,typename> class V,
+ typename K,
+ typename A
+>
+estimator<R,V,K,A>::estimator(K k):k_(k),a_(){}
+
+template<
+ typename R,
+ template<typename,typename,typename> class V,
+ typename K,
+ typename A
+>
+estimator<R,V,K,A>::estimator(K k,const A& a)
+:train_data_(),k_(k),a_(a){}
+
+template<
+ typename R,
+ template<typename,typename,typename> class V,
+ typename K,
+ typename A
+>
+estimator<R,V,K,A>::estimator(const estimator& that)
+:train_data_(that.train_data_),k_(that.k_),a_(that.a_){}
+
+template<
+ typename R,
+ template<typename,typename,typename> class V,
+ typename K,
+ typename A
+>
+estimator<R,V,K,A>&
+estimator<R,V,K,A>::operator=(const estimator& that){
+ if(&that = this){
+ this->train_data_ = that.train_data_;
+ this->k_ = that.k_;
+ this->a_ = that.a_;
+ }
+ return *this;
+}
+
+template<
+ typename R,
+ template<typename,typename,typename> class V,
+ typename K,
+ typename A
+>
+void estimator<R,V,K,A>::set_train_data(const R& train_data){
+ this->train_data_ = train_data;
+}
+
+// Evaluate
+template<
+ typename R,
+ template<typename,typename,typename> class V,
+ typename K,
+ typename A
+>
+ template<typename X>
+typename estimator<R,V,K,A>::template result<X>::type
+estimator<R,V,K,A>::operator()(const X& x)const{
+ typedef typename estimator<R,V,K,A>::template result<X>::type result_;
+ return std::for_each(
+ boost::begin( this->train_data() ),
+ boost::end( this->train_data() ),
+ result_( this->kernel(), x, this->a_)
+ );
+}
+
+// Access
+template<
+ typename R,
+ template<typename,typename,typename> class V,
+ typename K,
+ typename A
+>
+const R&
+estimator<R,V,K,A>::train_data()const{ return this->train_data_; }
+
+template<
+ typename R,
+ template<typename,typename,typename> class V,
+ typename K,
+ typename A
+>
+const K&
+estimator<R,V,K,A>::kernel()const{ return this->k_; }
+
+template<
+ typename R,
+ template<typename,typename,typename> class V,
+ typename K,
+ typename A
+>
+const A&
+estimator<R,V,K,A>::accumulator()const{
+ return this->a_;
+}
+
+}// kernel
+}// statistics
+}// boost
+
+#endif
\ No newline at end of file
Added: sandbox/statistics/kernel/boost/statistics/kernel/functional/include.hpp
==============================================================================
--- (empty file)
+++ sandbox/statistics/kernel/boost/statistics/kernel/functional/include.hpp 2009-08-28 18:41:48 EDT (Fri, 28 Aug 2009)
@@ -0,0 +1,18 @@
+///////////////////////////////////////////////////////////////////////////////
+// statistics::kernel::functional::include.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_KERNEL_FUNCTIONAL_INCLUDE_H_ER_2009
+#define BOOST_STATISTICS_KERNEL_FUNCTIONAL_INCLUDE_H_ER_2009
+
+#include <boost/statistics/kernel/functional/benchmark_nw.hpp>
+#include <boost/statistics/kernel/functional/benchmark_rp.hpp>
+#include <boost/statistics/kernel/functional/estimator.hpp>
+#include <boost/statistics/kernel/functional/nw_visitor.hpp>
+#include <boost/statistics/kernel/functional/nw_visitor_tuple.hpp>
+#include <boost/statistics/kernel/functional/rp_visitor.hpp>
+
+#endif
Added: sandbox/statistics/kernel/boost/statistics/kernel/functional/nw_visitor.hpp
==============================================================================
--- (empty file)
+++ sandbox/statistics/kernel/boost/statistics/kernel/functional/nw_visitor.hpp 2009-08-28 18:41:48 EDT (Fri, 28 Aug 2009)
@@ -0,0 +1,142 @@
+///////////////////////////////////////////////////////////////////////////////
+// statistics::kernel::functional::nw_visitor.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_KERNEL_FUNCTIONAL_NW_VISITOR_HPP_ER_2009
+#define BOOST_STATISTICS_KERNEL_FUNCTIONAL_NW_VISITOR_HPP_ER_2009
+#include <boost/type_traits/is_reference.hpp>
+#include <boost/mpl/not.hpp>
+#include <boost/call_traits.hpp>
+#include <boost/statistics/kernel/functional/detail/mean_accumulator.hpp>
+#include <boost/statistics/kernel/functional/rp_visitor.hpp>
+
+namespace boost{
+namespace statistics{
+namespace kernel{
+
+// Let E[Y|X=x] denote the conditional mean of Y given X = x
+//
+// This class is to be used to visit the training data, to obtain an estimator,
+// by the Nadaraya-Watson method, of E[Y|X=x]
+template<
+ typename K,
+ typename X,
+ typename A = typename
+ statistics::kernel::detail::mean_accumulator<typename K::result_type>::type
+>
+class nw_visitor{
+ public:
+ typedef rp_visitor<K,X,A> rp_visitor_type;
+ typedef typename rp_visitor_type::result_type result_type;
+ typedef K kernel_type;
+ typedef A accumulator_type;
+
+ //Construct
+ nw_visitor();
+ nw_visitor(typename call_traits<X>::param_type);
+ nw_visitor(
+ K k, // passing radius should call implicit conversion
+ typename call_traits<X>::param_type x
+ );
+ nw_visitor(
+ K k,
+ typename call_traits<X>::param_type,
+ const accumulator_type&
+ );
+ nw_visitor(const nw_visitor&);
+ nw_visitor& operator=(const nw_visitor&);
+
+ // Update
+ template<typename X1,typename Y1> // Training data point
+ result_type operator()(const X1& x1,const Y1& y1);
+
+ // Access
+ result_type unnormalized_estimate()const;
+ result_type normalizing_constant()const;
+ result_type estimate()const;
+
+ const A& accumulator()const;
+ const rp_visitor_type& rp_visitor()const;
+
+ private:
+ rp_visitor_type rp_visitor_;
+ A a_;
+};
+
+//Construction
+template<typename K,typename X,typename A>
+nw_visitor<K,X,A>::nw_visitor():rp_visitor_(),a_(){}
+
+template<typename K,typename X,typename A>
+nw_visitor<K,X,A>::nw_visitor(K k,typename call_traits<X>::param_type x)
+:rp_visitor_(k,x),a_(){}
+
+template<typename K,typename X,typename A>
+nw_visitor<K,X,A>::nw_visitor(
+ K k,typename call_traits<X>::param_type x,const A& a
+):rp_visitor_(k,x,a),a_(a){}
+
+template<typename K,typename X,typename A>
+nw_visitor<K,X,A>::nw_visitor(const nw_visitor& that)
+:rp_visitor_(that.rp_visitor_),a_(that.a_){}
+
+template<typename K,typename X,typename A>
+typename nw_visitor<K,X,A>::nw_visitor&
+nw_visitor<K,X,A>::operator=(const nw_visitor& that){
+ if(&that!=this){
+ rp_visitor_ = that.rp_visitor_;
+ a_ = that.a_;
+ }
+ return *this;
+}
+
+// Update
+template<typename K,typename X,typename A>
+template<typename X1,typename Y1>
+typename nw_visitor<K,X,A>::result_type
+nw_visitor<K,X,A>::operator()(const X1& x1,const Y1& y){
+ result_type w = (this->rp_visitor_(x1));
+ result_type wy = w * y;
+ this->a_(wy);
+ return wy;
+}
+
+// Access
+template<typename K,typename X,typename A>
+typename nw_visitor<K,X,A>::result_type
+nw_visitor<K,X,A>::unnormalized_estimate()const{
+ return accumulators::mean(
+ this->accumulator()
+ );
+}
+
+template<typename K,typename X,typename A>
+typename nw_visitor<K,X,A>::result_type
+nw_visitor<K,X,A>::normalizing_constant()const{
+ return (this->rp_visitor_).estimate();
+}
+
+template<typename K,typename X,typename A>
+typename nw_visitor<K,X,A>::result_type
+nw_visitor<K,X,A>::estimate()const{
+ return (this->unnormalized_estimate()/this->normalizing_constant());
+}
+
+template<typename K,typename X,typename A>
+const A& nw_visitor<K,X,A>::accumulator()const{ return this->a_; }
+
+template<typename K,typename X,typename A>
+const rp_visitor<K,X,A>&
+nw_visitor<K,X,A>::rp_visitor()const{
+ return (this->rp_visitor_);
+}
+
+
+}// kernel
+}// statistics
+}// boost
+
+#endif
\ No newline at end of file
Added: sandbox/statistics/kernel/boost/statistics/kernel/functional/nw_visitor_tuple.hpp
==============================================================================
--- (empty file)
+++ sandbox/statistics/kernel/boost/statistics/kernel/functional/nw_visitor_tuple.hpp 2009-08-28 18:41:48 EDT (Fri, 28 Aug 2009)
@@ -0,0 +1,126 @@
+///////////////////////////////////////////////////////////////////////////////
+// statistics::kernel::functional::nw_visitor.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_KERNEL_FUNCTIONAL_TUPLE_NW_VISITOR_HPP_ER_2009
+#define BOOST_STATISTICS_KERNEL_FUNCTIONAL_TUPLE_NW_VISITOR_HPP_ER_2009
+#include <boost/type_traits.hpp>
+#include <boost/mpl/assert.hpp>
+#include <boost/call_traits.hpp>
+#include <boost/binary_op/functional/untupler.hpp>
+#include <boost/binary_op/data/tuple_range.hpp>
+#include <boost/statistics/kernel/functional/detail/mean_accumulator.hpp>
+#include <boost/statistics/kernel/functional/nw_visitor.hpp>
+
+namespace boost{
+namespace statistics{
+namespace kernel{
+
+// This functor (f) is like nw_visitor, but rather than f(x,y),
+// the appropriate call is f(t), where t=(x,y) is a tuple.
+template<typename K,typename X, typename A = typename
+ statistics::kernel::detail::mean_accumulator<typename K::result_type>::type >
+class nw_visitor_tuple
+ : public binary_op::untupler<nw_visitor<K,X,A>,0,1>{
+ typedef nw_visitor<K,X,A> binary_;
+ typedef binary_op::untupler<binary_,0,1> super_t;
+
+ public:
+ typedef binary_ nw_visitor_type;
+ typedef typename binary_::rp_visitor_type rp_visitor_type;
+ typedef typename super_t::result_type result_type;
+
+ // Construction
+ nw_visitor_tuple(
+ K k,
+ typename call_traits<X>::param_type
+ );
+ nw_visitor_tuple(
+ K k,
+ typename call_traits<X>::param_type,
+ const A&
+ );
+ nw_visitor_tuple(const nw_visitor_tuple&);
+ nw_visitor_tuple& operator=(const nw_visitor_tuple&);
+
+ // Access
+ rp_visitor_type rp_visitor()const;
+ nw_visitor_type nw_visitor()const;
+ result_type rp_estimate()const;
+ result_type nw_estimate()const;
+ result_type estimate()const; //same as nw_estimate
+
+ private:
+ nw_visitor_tuple();
+};
+
+// Construction
+template<typename K,typename X,typename A>
+nw_visitor_tuple<K,X,A>::nw_visitor_tuple(
+ K k,
+ typename call_traits<X>::param_type x
+):super_t(binary_(k,x)){}
+
+ template<typename K,typename X,typename A>
+nw_visitor_tuple<K,X,A>::nw_visitor_tuple(
+ K k,
+ typename call_traits<X>::param_type x,
+ const A& a
+):super_t(binary_(k,x,a)){}
+
+
+template<typename K,typename X,typename A>
+nw_visitor_tuple<K,X,A>::nw_visitor_tuple(const nw_visitor_tuple& that)
+:super_t(static_cast<const super_t&>(that)){}
+
+template<typename K,typename X,typename A>
+typename nw_visitor_tuple<K,X,A>::nw_visitor_tuple&
+nw_visitor_tuple<K,X,A>::nw_visitor_tuple::operator=(
+ const nw_visitor_tuple& that
+){
+ if(&that!=this){
+ super_t::operator=(that);
+ }
+ return *this;
+}
+
+// Access
+
+template<typename K,typename X,typename A>
+typename nw_visitor_tuple<K,X,A>::nw_visitor_type
+nw_visitor_tuple<K,X,A>::nw_visitor()const{
+ return this->base();
+}
+
+template<typename K,typename X,typename A>
+typename nw_visitor_tuple<K,X,A>::rp_visitor_type
+nw_visitor_tuple<K,X,A>::rp_visitor()const{
+ return (this->nw_visitor()).rp_visitor();
+}
+
+template<typename K,typename X,typename A>
+typename nw_visitor_tuple<K,X,A>::result_type
+nw_visitor_tuple<K,X,A>::rp_estimate()const{
+ return (this->rp_visitor()).estimate();
+}
+
+template<typename K,typename X,typename A>
+typename nw_visitor_tuple<K,X,A>::result_type
+nw_visitor_tuple<K,X,A>::nw_estimate()const{
+ return (this->nw_visitor()).estimate();
+}
+
+template<typename K,typename X,typename A>
+typename nw_visitor_tuple<K,X,A>::result_type
+nw_visitor_tuple<K,X,A>::estimate()const{
+ return (this->nw_estimate());
+}
+
+}// kernel
+}// statistics
+}// boost
+
+#endif
\ No newline at end of file
Added: sandbox/statistics/kernel/boost/statistics/kernel/functional/rp_visitor.hpp
==============================================================================
--- (empty file)
+++ sandbox/statistics/kernel/boost/statistics/kernel/functional/rp_visitor.hpp 2009-08-28 18:41:48 EDT (Fri, 28 Aug 2009)
@@ -0,0 +1,149 @@
+///////////////////////////////////////////////////////////////////////////////
+// statistics::kernel::functional::rp_visitor.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_KERNEL_FUNCTIONAL_RP_VISITOR_H_ER_2009
+#define BOOST_STATISTICS_KERNEL_FUNCTIONAL_RP_VISITOR_H_ER_2009
+#include <boost/type_traits/is_reference.hpp>
+#include <boost/mpl/not.hpp>
+#include <boost/call_traits.hpp>
+#include <boost/statistics/kernel/functional/detail/mean_accumulator.hpp>
+//#include <boost/statistics/kernel/functional/detail/return_if.hpp>
+//#include <boost/statistics/kernel/functional/detail/range_difference.hpp>
+
+namespace boost{
+namespace statistics{
+namespace kernel{
+
+// Let p(x) denote the density of X = x
+//
+// This class is to be used to visit the training data, to obtain an estimator,
+// by the Rosenblatt-Parzen method, of p(x)
+template<
+ typename K,
+ typename X,
+ typename A = typename
+ statistics::kernel::detail::mean_accumulator<typename K::result_type>::type
+>
+class rp_visitor : K{ //, addable<rp_visitor<K,X,A> >{ //
+ typedef is_reference<X> is_ref_;
+ public:
+ typedef K kernel_type;
+ typedef A accumulator_type;
+ typedef typename K::result_type result_type;
+
+ // Construct
+ rp_visitor();
+ rp_visitor(typename call_traits<X>::param_type);
+ rp_visitor(
+ K k, // passing radius calls implicit conversion
+ typename call_traits<X>::param_type x
+ );
+ rp_visitor(
+ K k,
+ typename call_traits<X>::param_type,
+ const accumulator_type&
+ );
+ rp_visitor(const rp_visitor&);
+ rp_visitor& operator=(const rp_visitor&);
+
+ // Update
+ // Passing the training data x1 updates the estimator
+ template<typename X1> result_type operator()(const X1& x1);
+
+ // Access
+ typename call_traits<X>::const_reference x()const;
+ const A& accumulator()const;
+ const result_type& normalizing_constant()const;
+
+ result_type estimate()const;
+
+ private:
+ typename call_traits<X>::value_type x_;
+ A acc_;
+};
+
+//Construction
+template<typename K,typename X,typename A>
+rp_visitor<K,X,A>::rp_visitor(){
+ BOOST_MPL_ASSERT((
+ mpl::not_<is_ref_>
+ ));
+}
+
+template<typename K,typename X,typename A>
+rp_visitor<K,X,A>::rp_visitor(K k,typename call_traits<X>::param_type x)
+:K(k),x_(x),acc_(){}
+
+template<typename K,typename X,typename A>
+rp_visitor<K,X,A>::rp_visitor(
+ K k,
+ typename call_traits<X>::param_type x,
+ const A& a
+):K(k),x_(x),acc_(a){}
+
+template<typename K,typename X,typename A>
+rp_visitor<K,X,A>::rp_visitor(const rp_visitor& that)
+:K(static_cast<const K&>(that)),x_(that.x_),acc_(that.acc_){}
+
+template<typename K,typename X,typename A>
+typename rp_visitor<K,X,A>::rp_visitor&
+rp_visitor<K,X,A>::operator=(const rp_visitor& that){
+ if(&that!=this){
+ BOOST_MPL_ASSERT((mpl::not_<is_ref_>));
+ K::operator=(static_cast<const K&>(*that));
+ x_ = that.x_;
+ acc_ = that.acc_;
+ }
+ return *this;
+}
+
+// Evaluate
+template<typename K,typename X,typename A>
+template<typename X1>
+typename rp_visitor<K,X,A>::result_type
+rp_visitor<K,X,A>::operator()(const X1& x1){
+ const K& kernel = static_cast<const K&>(*this);
+ result_type t = kernel(x(),x1);
+ this->acc_(t);
+ return t;
+}
+
+
+// Access
+template<typename K,typename X,typename A>
+typename rp_visitor<K,X,A>::result_type
+rp_visitor<K,X,A>::estimate()const{
+ return accumulators::mean(
+ this->accumulator()
+ );
+}
+
+template<typename K,typename X,typename A>
+const A&
+rp_visitor<K,X,A>::accumulator()const{
+ return this->acc_;
+}
+
+template<typename K,typename X,typename A>
+const typename rp_visitor<K,X,A>::result_type&
+rp_visitor<K,X,A>::normalizing_constant()const{
+ const K& k = static_cast<const K&>(*this);
+ return k.normalizing_constant();
+}
+
+template<typename K,typename X,typename A>
+typename call_traits<X>::const_reference
+rp_visitor<K,X,A>::x()const{
+ return this->x_;
+}
+
+
+}// kernel
+}// statistics
+}// boost
+
+#endif
\ No newline at end of file
Added: sandbox/statistics/kernel/boost/statistics/kernel/include.hpp
==============================================================================
--- (empty file)
+++ sandbox/statistics/kernel/boost/statistics/kernel/include.hpp 2009-08-28 18:41:48 EDT (Fri, 28 Aug 2009)
@@ -0,0 +1,16 @@
+///////////////////////////////////////////////////////////////////////////////
+// statistics::kernel::include.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_KERNEL_INCLUDE_H_ER_2009
+#define BOOST_STATISTICS_KERNEL_INCLUDE_H_ER_2009
+
+#include <boost/statistics/kernel/bandwidth/include.hpp>
+#include <boost/statistics/kernel/functional/include.hpp>
+#include <boost/statistics/kernel/joint/include.hpp>
+#include <boost/statistics/kernel/scalar/include.hpp>
+
+#endif
\ No newline at end of file
Added: sandbox/statistics/kernel/boost/statistics/kernel/joint/crtp.hpp
==============================================================================
--- (empty file)
+++ sandbox/statistics/kernel/boost/statistics/kernel/joint/crtp.hpp 2009-08-28 18:41:48 EDT (Fri, 28 Aug 2009)
@@ -0,0 +1,25 @@
+///////////////////////////////////////////////////////////////////////////////
+// statistics::kernel::joint::crtp.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_KERNEL_JOINT_CRTP_H_ER_2009
+#define BOOST_STATISTICS_KERNEL_JOINT_CRTP_H_ER_2009
+
+namespace boost{
+namespace statistics{
+namespace kernel{
+namespace joint{
+
+ // This has no purpose other signifying that the derived class is
+ // a kernel, and it is currently not in use, but this may change.
+ template<typename D> struct crtp{};
+
+}
+}
+}
+}
+
+#endif
Added: sandbox/statistics/kernel/boost/statistics/kernel/joint/include.hpp
==============================================================================
--- (empty file)
+++ sandbox/statistics/kernel/boost/statistics/kernel/joint/include.hpp 2009-08-28 18:41:48 EDT (Fri, 28 Aug 2009)
@@ -0,0 +1,14 @@
+///////////////////////////////////////////////////////////////////////////////
+// statistics::kernel::joint::include.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_KERNEL_JOINT_INCLUDE_H_ER_2009
+#define BOOST_STATISTICS_KERNEL_JOINT_INCLUDE_H_ER_2009
+
+#include <boost/statistics/kernel/joint/crtp.hpp>
+#include <boost/statistics/kernel/joint/kernel_mono.hpp>
+
+#endif
\ No newline at end of file
Added: sandbox/statistics/kernel/boost/statistics/kernel/joint/kernel_mono.hpp
==============================================================================
--- (empty file)
+++ sandbox/statistics/kernel/boost/statistics/kernel/joint/kernel_mono.hpp 2009-08-28 18:41:48 EDT (Fri, 28 Aug 2009)
@@ -0,0 +1,126 @@
+///////////////////////////////////////////////////////////////////////////////
+// statistics::kernel::joint::kernel_mono.hpp //
+// //
+// Copyright 2009 Erwann Rogard. Distributed under the Boost //
+// Software License, Version 1.0. (See accompanying file //
+// LICEMSE_1_0.txt or copy at http://www.boost.org/LICEMSE_1_0.txt) //
+///////////////////////////////////////////////////////////////////////////////
+#ifndef BOOST_STATISTICS_KERNEL_JOINT_KERNEL_MONO_H_ER_2009
+#define BOOST_STATISTICS_KERNEL_JOINT_KERNEL_MONO_H_ER_2009
+#include <cmath>
+#include <numeric>
+#include <boost/lambda/lambda.hpp>
+#include <boost/vector_space/data/lazy_difference.hpp>
+#include <boost/statistics/kernel/joint/crtp.hpp>
+
+namespace boost{
+namespace statistics{
+namespace kernel{
+namespace joint{
+
+// Overview: Multivariate kernel, common bandwidth across coordinates.
+//
+// Notation: Let x denote a vector of size M,
+//
+// Usage:
+// typedef kernel_mono<scalar::gaussian<T>,M> mult_;
+// typedef mult_::result_type result_;
+// mult_ mult(bandwidth);
+// result_ r = mult(x);
+template<typename K,unsigned M>
+class kernel_mono : K{
+ public:
+ typedef typename K::result_type result_type;
+
+ //Construction
+ kernel_mono();
+ template<typename K1> kernel_mono(K1 k1); //pass bandwith, or kernel
+
+ // Evaluate
+ template<typename X> result_type profile(const X& x)const;
+ template<typename X> result_type operator()(const X& x)const;
+ template<typename X,typename X1>
+ result_type operator()(const X& x,const X1& x1)const;
+
+ // Access
+ static unsigned dimension;
+ result_type radius()const;
+ result_type normalizing_constant()const;
+
+ private:
+ result_type normalizing_constant_; //no ambiguity because inherit K privately
+ result_type comp_nc()const;
+};
+
+// Construction
+template<typename K,unsigned M>
+kernel_mono<K,M>::kernel_mono():K(),normalizing_constant_(comp_nc()){}
+
+template<typename K,unsigned M>
+template<typename K1>
+kernel_mono<K,M>::kernel_mono(K1 k1):K(k1),normalizing_constant_(comp_nc()){}
+
+// Evaluate
+template<typename K,unsigned M>
+template<typename X>
+typename kernel_mono<K,M>::result_type
+kernel_mono<K,M>::profile(const X& x)const{
+ static result_type init = static_cast<result_type>(0);
+ const K& k = static_cast<const K&>(*this);
+ result_type norm = std::accumulate(
+ boost::begin(x),
+ boost::end(x),
+ init,
+ ( lambda::_1 + (lambda::_2 * lambda::_2 ) )
+ );
+ norm = std::sqrt(norm);
+ return k.profile(norm);
+}
+
+template<typename K,unsigned M>
+template<typename X>
+typename kernel_mono<K,M>::result_type
+kernel_mono<K,M>::operator()(const X& x)const{
+ return this->profile(x) / this->normalizing_constant();
+}
+
+template<typename K,unsigned M>
+template<typename X,typename X1>
+typename kernel_mono<K,M>::result_type
+kernel_mono<K,M>::operator()(const X& x,const X1& x1)const{
+ typedef vector_space::lazy_difference<X,X1> diff_;
+ typedef typename range_size<X>::type size_type;
+ BOOST_ASSERT(size(x) == static_cast<size_type>(size(x1)));
+ BOOST_ASSERT(size(x) == static_cast<size_type>(M));
+ diff_ diff(x,x1);
+ return (*this)(diff);
+}
+
+// Access
+template<typename K,unsigned M> unsigned kernel_mono<K,M>::dimension = M;
+
+template<typename K,unsigned M>
+typename kernel_mono<K,M>::result_type
+kernel_mono<K,M>::radius()const{
+ const K& k = static_cast<const K&>(*this);
+ return k.radius();
+}
+
+template<typename K,unsigned M>
+typename kernel_mono<K,M>::result_type
+kernel_mono<K,M>::normalizing_constant()const{ return normalizing_constant_; }
+
+template<typename K,unsigned M>
+typename kernel_mono<K,M>::result_type
+kernel_mono<K,M>::comp_nc()const{
+ const K& k = static_cast<const K&>(*this);
+ static result_type m = static_cast<result_type>(M);
+ result_type nc = k.normalizing_constant();
+ return std::pow(nc,m);
+}
+
+}// joint
+}// kernel
+}// statistics
+}// boost
+#endif
\ No newline at end of file
Added: sandbox/statistics/kernel/boost/statistics/kernel/scalar/crtp.hpp
==============================================================================
--- (empty file)
+++ sandbox/statistics/kernel/boost/statistics/kernel/scalar/crtp.hpp 2009-08-28 18:41:48 EDT (Fri, 28 Aug 2009)
@@ -0,0 +1,116 @@
+///////////////////////////////////////////////////////////////////////////////
+// statistics::kernel::scalar::crtp.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_KERNEL_SCALAR_CRTP_H_ER_2009
+#define BOOST_STATISTICS_KERNEL_SCALAR_CRTP_H_ER_2009
+#include <boost/format.hpp>
+
+namespace boost{
+namespace statistics{
+namespace kernel{
+namespace scalar{
+
+// Overview: provides common operations for a kernel whose implementation
+// is specified by a derived class, D.
+//
+// Let d note an instance of D, and T == D::result_type
+//
+// Requirements
+// crtp<D,T> is a base of D
+// D::result_type is defined
+// D::core_profile(const T&) returns an object of type T
+// D::core_nc returns an object of type T (normalzing c)
+// D d; Constructor Forwards to crtp<D,T>
+// D d(r) Constructor Forwards to crtp<D,T>
+template<typename D,typename T>
+class crtp{
+ public:
+ typedef T result_type;
+
+ // Construction
+ crtp();
+ crtp(const result_type& bandwidth);
+
+ // Evaluate
+ template<typename X> result_type profile(const X& x)const;
+ template<typename X> result_type operator()(const X& x)const;
+ template<typename X,typename X1>
+ result_type operator()(const X& x,const X1& x1)const;
+
+ // Access
+ result_type bandwidth()const;
+ result_type normalizing_constant()const;
+
+ private:
+ result_type bandwidth_;
+ result_type normalizing_constant_;
+ result_type comp_nc(result_type bandwidth);
+};
+
+template<typename D,typename T>
+std::ostream& operator<<(std::ostream& out,const crtp<D,T>& k){
+ const char* str = "(%1%,%2%)";
+ format f(str); f%k.bandwidth()%k.normalizing_constant();
+ out << f.str();
+ return out;
+}
+
+// Construction
+template<typename D,typename T> crtp<D,T>::crtp()
+:bandwidth_(static_cast<result_type>(1)),
+normalizing_constant_(this->comp_nc(bandwidth())){}
+
+template<typename D,typename T> crtp<D,T>::crtp(const result_type& bandwidth)
+:bandwidth_(bandwidth),
+normalizing_constant_(this->comp_nc(this->bandwidth())){}
+
+
+// Evaluate
+template<typename D,typename T>
+template<typename X>
+typename crtp<D,T>::result_type
+crtp<D,T>::profile(const X& x)const{
+ result_type u = x / this->bandwidth();
+ return D::core_profile(u);
+}
+
+template<typename D,typename T>
+template<typename X>
+typename crtp<D,T>::result_type
+crtp<D,T>::operator()(const X& x)const{
+ return ( this->profile(x) ) / ( this->normalizing_constant()) ;
+}
+
+template<typename D,typename T>
+template<typename X,typename X1>
+typename crtp<D,T>::result_type
+crtp<D,T>::operator()(const X& x,const X1& x1)const{
+ return (*this)(x-x1);
+}
+
+// Access
+template<typename D,typename T>
+typename crtp<D,T>::result_type crtp<D,T>::bandwidth()const{
+ return this->bandwidth_;
+}
+
+template<typename D,typename T>
+typename crtp<D,T>::result_type crtp<D,T>::normalizing_constant()const{
+ return this->normalizing_constant_;
+}
+
+// Implem
+template<typename D,typename T>
+typename crtp<D,T>::result_type crtp<D,T>::comp_nc(result_type bandwidth){
+ return D::core_nc * bandwidth;
+}
+
+}// scalar
+}// kernel
+}// statistics
+}// boost
+#endif
\ No newline at end of file
Added: sandbox/statistics/kernel/boost/statistics/kernel/scalar/gaussian.hpp
==============================================================================
--- (empty file)
+++ sandbox/statistics/kernel/boost/statistics/kernel/scalar/gaussian.hpp 2009-08-28 18:41:48 EDT (Fri, 28 Aug 2009)
@@ -0,0 +1,53 @@
+///////////////////////////////////////////////////////////////////////////////
+// statistics::kernel::scalar::gaussian.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_KERNEL_SCALAR_GAUSSIAN_H_ER_2009
+#define BOOST_STATISTICS_KERNEL_SCALAR_GAUSSIAN_H_ER_2009
+#include <cmath>
+#include <boost/statistics/kernel/scalar/crtp.hpp>
+#include <boost/math/constants/constants.hpp>
+
+namespace boost{
+namespace statistics{
+namespace kernel{
+namespace scalar{
+
+template<typename T>
+struct gaussian : scalar::crtp<gaussian<T>,T >{
+ typedef gaussian<T> this_type;
+ typedef scalar::crtp<this_type,T> crtp_;
+ public:
+ typedef typename crtp_::result_type result_type;
+
+ gaussian();
+ gaussian(const result_type& bandwidth);
+
+ static result_type core_profile(const result_type& x);
+ static result_type core_nc;
+};
+
+template<typename T> gaussian<T>::gaussian():crtp_(){}
+template<typename T> gaussian<T>::gaussian(const result_type& bandwidth)
+:crtp_(bandwidth){}
+
+template<typename T>
+typename gaussian<T>::result_type
+gaussian<T>::core_profile(const result_type& x){
+ static result_type two = static_cast<T>(2);
+ return exp(- x * x / two);
+}
+
+template<typename T>
+typename gaussian<T>::result_type
+gaussian<T>::core_nc = math::constants::root_two_pi<T>();
+
+}// scalar
+}// kernel
+}// statistics
+}// boost
+
+#endif
\ No newline at end of file
Added: sandbox/statistics/kernel/boost/statistics/kernel/scalar/include.hpp
==============================================================================
--- (empty file)
+++ sandbox/statistics/kernel/boost/statistics/kernel/scalar/include.hpp 2009-08-28 18:41:48 EDT (Fri, 28 Aug 2009)
@@ -0,0 +1,14 @@
+///////////////////////////////////////////////////////////////////////////////
+// statistics::kernel::scalar::include.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_KERNEL_SCALAR_INCLUDE_H_ER_2009
+#define BOOST_STATISTICS_KERNEL_SCALAR_INCLUDE_H_ER_2009
+
+#include <boost/statistics/kernel/scalar/crtp.hpp>
+#include <boost/statistics/kernel/scalar/gaussian.hpp>
+
+#endif
\ No newline at end of file