|
|
Boost-Commit : |
Subject: [Boost-commit] svn:boost r55848 - in sandbox/statistics/kernel/libs: . statistics statistics/kernel statistics/kernel/doc statistics/kernel/example statistics/kernel/src
From: erwann.rogard_at_[hidden]
Date: 2009-08-28 18:46:31
Author: e_r
Date: 2009-08-28 18:46:30 EDT (Fri, 28 Aug 2009)
New Revision: 55848
URL: http://svn.boost.org/trac/boost/changeset/55848
Log:
replacing libs/kernel
Added:
sandbox/statistics/kernel/libs/
sandbox/statistics/kernel/libs/statistics/
sandbox/statistics/kernel/libs/statistics/kernel/
sandbox/statistics/kernel/libs/statistics/kernel/doc/
sandbox/statistics/kernel/libs/statistics/kernel/doc/readme.txt (contents, props changed)
sandbox/statistics/kernel/libs/statistics/kernel/example/
sandbox/statistics/kernel/libs/statistics/kernel/example/benchmark_scalar.cpp (contents, props changed)
sandbox/statistics/kernel/libs/statistics/kernel/example/benchmark_scalar.h (contents, props changed)
sandbox/statistics/kernel/libs/statistics/kernel/example/kernel_mono_rp.cpp (contents, props changed)
sandbox/statistics/kernel/libs/statistics/kernel/example/kernel_mono_rp.h (contents, props changed)
sandbox/statistics/kernel/libs/statistics/kernel/example/scalar_nw.cpp (contents, props changed)
sandbox/statistics/kernel/libs/statistics/kernel/example/scalar_nw.h (contents, props changed)
sandbox/statistics/kernel/libs/statistics/kernel/example/scalar_rp.cpp (contents, props changed)
sandbox/statistics/kernel/libs/statistics/kernel/example/scalar_rp.h (contents, props changed)
sandbox/statistics/kernel/libs/statistics/kernel/src/
sandbox/statistics/kernel/libs/statistics/kernel/src/main.cpp (contents, props changed)
Added: sandbox/statistics/kernel/libs/statistics/kernel/doc/readme.txt
==============================================================================
--- (empty file)
+++ sandbox/statistics/kernel/libs/statistics/kernel/doc/readme.txt 2009-08-28 18:46:30 EDT (Fri, 28 Aug 2009)
@@ -0,0 +1,162 @@
+//////////////////////////////////////////////////////////////////////////////
+// statistics::kernel::doc::readme //
+// //
+// (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) //
+//////////////////////////////////////////////////////////////////////////////
+
+[ Contact ]
+
+erwann.rogard_at_[hidden]
+
+[ Overview ]
+
+These are C++ classes for kernel based density estimation and conditional mean
+estimation by the method of Rosenblatt-Parzen and Nadaraya-Watson, respectively.
+
+[ Compiler ]
+
+gcc version i686-apple-darwin9-gcc-4.0.1 (GCC) 4.0.1
+
+[ Dependencies ]
+
+boost_1_39_0
+/sandbox/statistics/vector_space/
+/sandbox/statistics/binary_op/
+/sandbox/statistics/arithmetic/
+/sandbox/statistics/scalar_dist/
+/sandbox/statistics/mpl/
+/sandbox/statistics/standard_distribution/
+/sandbox/statistics/kernel/
+
+[ Sources ]
+
+http://en.wikipedia.org/wiki/Kernel_(statistics)
+http://en.wikipedia.org/wiki/Kernel_density_estimation
+@book{citeulike:300228,
+ author = {Silverman, B. W.},
+ howpublished = {Hardcover},
+ isbn = {0412246201},
+ keywords = {density-estimation},
+ month = {April},
+ posted-at = {2008-01-13 17:43:25},
+ priority = {2},
+ publisher = {{Chapman \& Hall/CRC}},
+ title = {Density Estimation for Statistics and Data Analysis},
+ year = {1986}
+}
+
+[ History ]
+
+July 2009 : Current version
+March 2009: improved_fast_gauss_transform is now deprecated
+
+[Design]
+
+Let h denote bandwidth, d dimension, ||.|| norm.
+
+To perform density or conditional mean estimation, we are interested in
+operations that involve an abitrary bandwidth binary kernel, K[h](x0,x1),
+and/or possibly its kernel_monovariate extension K[d,h](x0,x1).
+
+Our design rests on the following relationships (see [Math representation])
+ - standardized kernel (bandwidth = 1) and kernels indexed by a bandwidth
+ - unary kernels, K(x), and binary kernels, K(x0,x1)
+ - univariate and kernel_monovariate kernels.
+
+The client need only specify the standardized kernel, such as
+scalar::gaussian, while a crtp mechanism, implemented in scalar::crtp,
+adds the remaining functionality.
+
+The class template statistics::kernel::kernel_mono, parameterized by an arbitrary unavariate
+kernel provides a kernel_monovariate analogue to the latter.
+
+The classes rp_visitor and nw_visitor are implementations for the Rosenblatt-
+Parzen and Nadaraya-Watson estimators for a given test value, x.
+
+Iterating either of the above estimators over a training sample is implemented
+by statistics::kernel::estimator, which, in turn, can be used in a std::for_each construct
+to iterate over test data.
+
+Bandwidth selection for density estimation is provided by
+bandwith::normal_distribution.
+
+[ Directory structure ]
+
+/bandwidth
+ Tools for bandwidth selection
+
+/functional
+
+ Tools for iterating over a training or test data
+
+/scalar
+
+ Univariate kernels
+
+/joint
+
+ Multivariate kernels
+
+[TODO]
+
+- Bandwidth selection by cross validation and other methods.
+- Provide a IFGT version of the above (see our deprecated
+improved_fast_gauss_transform package)
+- Kernel regression etc.
+
+[ Math representation ]
+
+Univariate
+ Unary:
+ Unnormalized: Ku(x) : [0,inf) -> [0,inf)
+ Normalized: K(x) : [0,inf) -> [0,1)
+ K(x) = Ku(x) / c such that integral K(x) dx = 1
+ Parameterized
+ by bandwidth K[h](x) = K(x/h) / h
+
+ e.g.
+ Ku(x) : exp(-x^2/2)
+ c : sqrt(2 pi)
+ K(x) : exp(-x^2/2) / sqrt(2 pi)
+ K[h](x) : exp(-(x/h)^2 /2) / ( sqrt(2 pi) * h)
+
+Binary:
+ K(x0,x1) = K(|x0-x1|)
+
+Multivariate:
+ Ku[d](x) : Ku(||x||)
+
+ e.g
+ Ku[d](x) : exp(-||x||^2) : Ku(||x||)
+ c[d] : (sqrt(2 pi))^d : c^d
+ K[d](x) : Ku[d](x)/c[d]
+
+ K[d,h](x) : exp(-||x/h||^2)/(c * h)^d
+ = Ku(||x||/h) / (c * h)^d
+
+[Output]
+
+Here's a copy of the output from libs/src/main.cpp:
+
+-> example_scalar_rp : 0.256009 0.390986 0.229151 0.261866 0.232843 0.228334
+0.185519 0.382394 0.388576 0.383203 <-
+-> example_scalar_nw : <-
+-> example_kernel_mono_rp : <-
+-> example_benchmark_scalar : m = 10000
+rp : (n, avg_t,e)nw : (n, avg_t,e_p,e_y)
+rp : (2, 4.0672e-06, 0.173297) nw : (2, 8.3822e-06, 0.173297, 1.2265)
+rp : (4, 2.859e-06, 0.138165) nw : (4, 5.86085e-06, 0.138165, 1.06847)
+rp : (8, 3.01747e-06, 0.104602) nw : (8, 6.289e-06, 0.104602, 0.831601)
+rp : (16, 3.96233e-06, 0.0978568) nw : (16, 8.16075e-06, 0.0978568, 0.702637)
+rp : (32, 5.94066e-06, 0.0818718) nw : (32, 1.20646e-05, 0.0818718, 0.595622)
+rp : (64, 9.40283e-06, 0.0719671) nw : (64, 1.95258e-05, 0.0719671, 0.531848)
+rp : (128, 1.57752e-05, 0.0634377) nw : (128, 3.26977e-05, 0.0634377, 0.477009)
+rp : (256, 2.72618e-05, 0.0565864) nw : (256, 5.64065e-05, 0.0565864, 0.432414)
+rp : (512, 4.81542e-05, 0.0511678) nw : (512, 0.00010056, 0.0511678, 0.398781)
+rp : (1024, 8.67954e-05, 0.0469132) nw : (1024, 0.000179972, 0.0469132, 0.371483)
+<-
+
+
Added: sandbox/statistics/kernel/libs/statistics/kernel/example/benchmark_scalar.cpp
==============================================================================
--- (empty file)
+++ sandbox/statistics/kernel/libs/statistics/kernel/example/benchmark_scalar.cpp 2009-08-28 18:46:30 EDT (Fri, 28 Aug 2009)
@@ -0,0 +1,274 @@
+///////////////////////////////////////////////////////////////////////////////
+// statistics::kernel::example::benchmark_scalar.cpp //
+// //
+// 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) //
+///////////////////////////////////////////////////////////////////////////////
+#include <cmath>
+#include <vector>
+#include <algorithm>
+#include <iterator>
+#include <boost/function.hpp>
+#include <boost/range.hpp>
+#include <boost/foreach.hpp>
+#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/normal_distribution.hpp>
+#include <boost/random/variate_generator.hpp>
+#include <boost/math/special_functions/fpclassify.hpp> //needed?
+#include <boost/math/tools/precision.hpp>
+
+// Order of the files matters!
+#include <boost/standard_distribution/distributions/normal.hpp>
+#include <boost/scalar_dist/fun_wrap/pdf.hpp>
+#include <boost/dist_random/distributions/normal.hpp>
+#include <boost/dist_random/random/generate_n.hpp>
+
+//#include <boost/scalar_dist/fun_wrap/pdf.hpp>
+//#include <boost/scalar_dist/meta/delegate.hpp>
+//#include <boost/scalar_dist/algorithm/transform.hpp>
+
+#include <boost/binary_op/data/tuple_range.hpp>
+
+#include <boost/statistics/kernel/include.hpp>
+
+//#include <boost/statistics/kernel/bandwidth/normal_distribution.hpp>
+//#include <boost/statistics/kernel/scalar/gaussian.hpp>
+//#include <boost/statistics/kernel/functional/nw_visitor_tuple.hpp>
+//#include <boost/statistics/kernel/functional/estimator.hpp>
+
+#include <boost/statistics/goodness_of_fit/include.hpp>
+
+#include <libs/statistics/kernel/example/benchmark_scalar.h>
+
+void example_benchmark_scalar(std::ostream& out){
+ out << "-> example_benchmark_scalar : ";
+ out.flush();
+
+ using namespace boost;
+ namespace stat = statistics;
+ namespace gof = stat::goodness_of_fit;
+
+ // Types
+ typedef double val_;
+ typedef function<val_(const val_&)> y_fun_;
+ typedef std::vector<val_> vals_;
+ typedef std::vector<unsigned> uints_;
+ typedef mt19937 urng_;
+ typedef normal_distribution<val_> rnorm_;
+ typedef variate_generator<urng_&,rnorm_> gen_;
+ typedef statistics::kernel::scalar::gaussian<val_> gauss_k_;
+ typedef statistics::kernel::rp_visitor<gauss_k_,val_> rp_visitor_;
+ typedef math::normal_distribution<val_> mnorm_;
+
+ // Constants
+ const unsigned max_j = 5; // # partitions, k = 2^j
+ const unsigned dim = 1;
+ typedef statistics::kernel::normal_distribution<dim,val_> opt_;
+
+ const unsigned total_n = 1e4;
+ const val_ mu = 0.0;
+ const val_ sigma = 2.0;
+
+ // Initialization
+ y_fun_ y_fun = fabs;
+
+ vals_ vec_x( total_n );
+ vals_ vec_y( total_n );
+ vals_ vec_pdf( total_n );
+ urng_ urng;
+ rnorm_ rnorm(mu,sigma);
+ gen_ gen(urng,rnorm);
+
+ { // Data generation
+ mnorm_ mnorm(mu,sigma);
+ generate_function_n<math::fun_wrap::pdf_>(
+ boost::begin(vec_x),
+ boost::begin(vec_pdf),
+ total_n,
+ mnorm,
+ urng
+ );
+ transform(
+ boost::begin(vec_x),
+ boost::end(vec_x),
+ boost::begin(vec_y),
+ y_fun
+ );
+ }
+ { // K-fold NW-estimation
+ typedef binary_op::tuple_range<const vals_&,const vals_&>
+ meta_range_xy_;
+ typedef meta_range_xy_::type range_xy_;
+ typedef range_value<range_xy_>::type xy_ref_cv_;
+ typedef binary_op::tuple_remove_ref_cv<
+ xy_ref_cv_
+ >::type xy_;
+
+ range_xy_ range_xy = meta_range_xy_::make(vec_x,vec_y);
+ typedef gof::k_fold_data< xy_ > k_fold_;
+ typedef k_fold_::range_train_data_type range_train_;
+
+ typedef statistics::kernel::estimator<
+ range_train_,
+ statistics::kernel::nw_visitor_tuple,
+ gauss_k_
+ > estimator_;
+
+ vals_ mse;
+ vals_ mae;
+ uints_ ns;
+ vals_ bandwidths;
+ vals_ vec_est_y( total_n );
+ unsigned k = 1;
+ for(unsigned j = 1; j<max_j; j++){
+ k *= 2;
+ k_fold_ k_fold(
+ k,
+ boost::begin(range_xy),
+ boost::end(range_xy)
+ );
+ ns.push_back( k_fold.n() );
+ // Warning: This bandwidth is not optimal for nw (but better than
+ // a wild guess).
+ bandwidths.push_back(
+ opt_::rp_bandwidth(sigma, k_fold.n())
+ );
+ estimator_ e(
+ bandwidths.back()
+ );
+ gof::regression_estimate(
+ k_fold,
+ e,
+ boost::begin( vec_est_y )
+ );
+ mse.push_back(
+ gof::sqrt_mse(
+ boost::begin( vec_est_y ),
+ boost::end( vec_est_y ),
+ boost::begin( vec_y )
+ )
+ );
+ mae.push_back(
+ gof::mean_abs_error(
+ boost::begin( vec_est_y ),
+ boost::end( vec_est_y ),
+ boost::begin( vec_y )
+ )
+ );
+ }
+
+ out << std::endl << "NW estimation";
+ out << std::endl <<"bandwidth : ";
+ std::copy(
+ boost::rbegin( bandwidths ),
+ boost::rend( bandwidths ),
+ std::ostream_iterator<val_>(out," ")
+ );
+ out << std::endl <<"n : ";
+ std::copy(
+ boost::rbegin( ns ),
+ boost::rend( ns ),
+ std::ostream_iterator<val_>(out," ")
+ );
+ out << std::endl <<"sqrt_mse : ";
+ std::copy(
+ boost::rbegin( mse ),
+ boost::rend( mse ),
+ std::ostream_iterator<val_>(out," ")
+ );
+ out << std::endl <<"mae : ";
+ std::copy(
+ boost::rbegin( mae ),
+ boost::rend( mae ),
+ std::ostream_iterator<val_>(out," ")
+ );
+
+ }
+ { // K-fold RP-estimation
+
+ typedef gof::k_fold_data< val_ > k_fold_;
+ typedef k_fold_::range_train_data_type range_train_;
+
+ typedef statistics::kernel::estimator<
+ range_train_,
+ statistics::kernel::rp_visitor,
+ gauss_k_
+ > estimator_;
+
+ vals_ mse;
+ vals_ mae;
+ uints_ ns;
+ vals_ bandwidths;
+ vals_ vec_est_pdf( total_n );
+ unsigned k = 1;
+ for(unsigned j = 1; j<max_j; j++){
+ k *= 2;
+ k_fold_ k_fold(
+ k,
+ boost::begin( vec_x ),
+ boost::end( vec_x )
+ );
+ ns.push_back( k_fold.n() );
+ // Note: This bandwidth is not optimal for density estimation
+ // using a gaussian kernel, under normality of data
+ bandwidths.push_back(
+ opt_::rp_bandwidth(sigma, k_fold.n())
+ );
+ estimator_ e(
+ bandwidths.back()
+ );
+ gof::marginal_estimate(
+ k_fold,
+ e,
+ boost::begin( vec_est_pdf )
+ );
+ mse.push_back(
+ gof::sqrt_mse(
+ boost::begin( vec_est_pdf ),
+ boost::end( vec_est_pdf ),
+ boost::begin( vec_pdf )
+ )
+ );
+ mae.push_back(
+ gof::mean_abs_error(
+ boost::begin( vec_est_pdf ),
+ boost::end( vec_est_pdf ),
+ boost::begin( vec_pdf )
+ )
+ );
+ }
+ out << std::endl << "RP estimation ";
+ out << std::endl <<"bandwidth : ";
+ std::copy(
+ boost::rbegin( bandwidths ),
+ boost::rend( bandwidths ),
+ std::ostream_iterator<val_>(out," ")
+ );
+ out << std::endl <<"n : ";
+ std::copy(
+ boost::rbegin( ns ),
+ boost::rend( ns ),
+ std::ostream_iterator<val_>(out," ")
+ );
+ out << std::endl <<"sqrt_mse : ";
+ std::copy(
+ boost::rbegin( mse ),
+ boost::rend( mse ),
+ std::ostream_iterator<val_>(out," ")
+ );
+ out << std::endl <<"mae : ";
+ std::copy(
+ boost::rbegin( mae ),
+ boost::rend( mae ),
+ std::ostream_iterator<val_>(out," ")
+ );
+
+ }
+ // DO the same for rp
+
+ out << "<-" << std::endl;
+}
+
+
+
Added: sandbox/statistics/kernel/libs/statistics/kernel/example/benchmark_scalar.h
==============================================================================
--- (empty file)
+++ sandbox/statistics/kernel/libs/statistics/kernel/example/benchmark_scalar.h 2009-08-28 18:46:30 EDT (Fri, 28 Aug 2009)
@@ -0,0 +1,14 @@
+///////////////////////////////////////////////////////////////////////////////
+// statistics::kernel::example::benchmark_scalar.h //
+// //
+// 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 LIBS_KERNEL_EXAMPLE_BENCHMARK_SCALAR_H_ER_2009
+#define LIBS_KERNEL_EXAMPLE_BENCHMARK_SCALAR_H_ER_2009
+#include <ostream>
+
+void example_benchmark_scalar(std::ostream&);
+
+#endif
\ No newline at end of file
Added: sandbox/statistics/kernel/libs/statistics/kernel/example/kernel_mono_rp.cpp
==============================================================================
--- (empty file)
+++ sandbox/statistics/kernel/libs/statistics/kernel/example/kernel_mono_rp.cpp 2009-08-28 18:46:30 EDT (Fri, 28 Aug 2009)
@@ -0,0 +1,94 @@
+///////////////////////////////////////////////////////////////////////////////
+// statistics::kernel::example::kernel_mono_rp.cpp //
+// //
+// 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) //
+///////////////////////////////////////////////////////////////////////////////
+#include <vector>
+#include <algorithm>
+#include <iterator>
+#include <boost/range.hpp>
+#include <boost/foreach.hpp>
+#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/normal_distribution.hpp>
+#include <boost/random/variate_generator.hpp>
+#include <boost/math/special_functions/fpclassify.hpp> //needed?
+#include <boost/math/tools/precision.hpp>
+
+#include <boost/statistics/kernel/scalar/gaussian.hpp>
+#include <boost/statistics/kernel/joint/kernel_mono.hpp>
+#include <boost/statistics/kernel/functional/rp_visitor.hpp>
+#include <boost/statistics/kernel/functional/estimator.hpp>
+#include <libs/statistics/kernel/example/scalar_rp.h>
+
+void example_kernel_mono_rp(std::ostream& out){
+ out << "-> example_kernel_mono_rp : ";
+ using namespace boost;
+
+ // Types
+ typedef double val_;
+ typedef std::vector<val_> vec_;
+ typedef std::vector<vec_> mat_;
+ typedef mt19937 urng_;
+ typedef normal_distribution<val_> norm_;
+ typedef variate_generator<urng_&,norm_> gen_;
+ typedef statistics::kernel::scalar::gaussian<val_> gauss_k_;
+
+ const unsigned dim = 2;
+ typedef statistics::kernel::joint::kernel_mono<gauss_k_,dim> kernel_mono_k_;
+ // NB const vec_&, not vec_
+ typedef statistics::kernel::rp_visitor<kernel_mono_k_,const vec_&> rp_visitor_;
+
+ // Constants
+ const val_ bandwidth = 0.5;
+ const val_ eps = math::tools::epsilon<val_>();
+ const unsigned n = 10;
+
+ // Generate sample
+ mat_ vec_x; vec_x.reserve(n);
+ vec_ vec_rp; vec_rp.reserve(n);
+ urng_ urng;
+ norm_ norm;
+ gen_ gen(urng,norm);
+ for(unsigned i = 0; i<n; i++){
+ vec_ tmp(dim);
+ std::generate_n(
+ begin(tmp),
+ dim,
+ gen
+ );
+ vec_x.push_back( tmp );
+ }
+
+ kernel_mono_k_ kernel_mono_k(bandwidth);
+
+ kernel_mono_k(vec_x[0],vec_x[1]);
+ // Density estimate for each x in vec_x using vec_x as the sample
+ BOOST_FOREACH(const vec_& x,vec_x){
+ val_ rp = std::for_each(
+ begin(vec_x),
+ end(vec_x),
+ rp_visitor_(bandwidth,x)
+ ).estimate();
+ vec_rp.push_back(rp);
+ }
+ typedef sub_range<mat_> sub_;
+ typedef statistics::kernel::estimator<
+ sub_,
+ statistics::kernel::rp_visitor,
+ kernel_mono_k_
+ > estimator_;
+ estimator_ estimator(bandwidth);
+ statistics::train(estimator,sub_(vec_x));
+ vec_ vec_rp2; vec_rp2.reserve(n);
+
+ // Same as previous but calls estimator instead of for_each
+ for(unsigned i = 0; i<n; i++){
+ vec_ x = vec_x[i];
+ val_ rp = vec_rp[i];
+ val_ rp2 = estimator(x).estimate();
+ BOOST_ASSERT(fabs(rp-rp2)<eps);
+ }
+ out << "<-" << std::endl;
+}
\ No newline at end of file
Added: sandbox/statistics/kernel/libs/statistics/kernel/example/kernel_mono_rp.h
==============================================================================
--- (empty file)
+++ sandbox/statistics/kernel/libs/statistics/kernel/example/kernel_mono_rp.h 2009-08-28 18:46:30 EDT (Fri, 28 Aug 2009)
@@ -0,0 +1,14 @@
+///////////////////////////////////////////////////////////////////////////////
+// statistics::kernel::example::kernel_mono_rp.h //
+// //
+// 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 LIBS_KERNEL_EXAMPLE_MULTI_RP_H_ER_2009
+#define LIBS_KERNEL_EXAMPLE_MULTI_RP_H_ER_2009
+#include <ostream>
+
+void example_kernel_mono_rp(std::ostream&);
+
+#endif
\ No newline at end of file
Added: sandbox/statistics/kernel/libs/statistics/kernel/example/scalar_nw.cpp
==============================================================================
--- (empty file)
+++ sandbox/statistics/kernel/libs/statistics/kernel/example/scalar_nw.cpp 2009-08-28 18:46:30 EDT (Fri, 28 Aug 2009)
@@ -0,0 +1,118 @@
+///////////////////////////////////////////////////////////////////////////////
+// statistics::kernel::example::scalar_nw.cpp //
+// //
+// 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) //
+///////////////////////////////////////////////////////////////////////////////
+#include <vector>
+#include <algorithm>
+#include <iterator>
+#include <boost/range.hpp>
+#include <boost/foreach.hpp>
+#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/normal_distribution.hpp>
+#include <boost/random/variate_generator.hpp>
+#include <boost/math/special_functions/fpclassify.hpp> //needed?
+#include <boost/math/tools/precision.hpp>
+#include <boost/typeof/typeof.hpp>
+
+#include <boost/binary_op/data/tuple_range.hpp>
+//#include <boost/statistics/estimator_concept/trainable_estimator/concept.hpp>
+#include <boost/statistics/kernel/scalar/gaussian.hpp>
+#include <boost/statistics/kernel/functional/nw_visitor_tuple.hpp>
+#include <boost/statistics/kernel/functional/estimator.hpp>
+#include <libs/statistics/kernel/example/scalar_nw.h>
+
+void example_scalar_nw(std::ostream& out){
+ out << "-> example_scalar_nw : ";
+ using namespace boost;
+
+ // Types
+ typedef double val_;
+ typedef std::vector<val_> vec_;
+ typedef mt19937 urng_;
+ typedef normal_distribution<val_> norm_;
+ typedef variate_generator<urng_&,norm_> gen_;
+ typedef statistics::kernel::scalar::gaussian<val_> gauss_k_;
+ typedef statistics::kernel::nw_visitor_tuple<gauss_k_,val_>
+ nw_visitor_tuple_;
+ typedef nw_visitor_tuple_::nw_visitor_type nw_visitor_;
+ typedef nw_visitor_tuple_::rp_visitor_type rp_visitor_;
+
+ // Constants
+ const val_ bandwidth = 0.5;
+ const val_ eps = math::tools::epsilon<val_>();
+ const unsigned n = 10;
+
+ // Initialization
+ vec_ vec_x(n);
+ vec_ vec_y(n,static_cast<val_>(1));
+ vec_ vec_rp; vec_rp.reserve(n);
+ vec_ vec_nw; vec_nw.reserve(n);
+ urng_ urng;
+ norm_ norm;
+ gen_ gen(urng,norm);
+ std::generate_n(
+ begin(vec_x),
+ n,
+ gen
+ );
+
+ // Computes a conditional mean estimate (nw) for each x in vec_x using
+ // a sequence of (x,y) tuples constructed from (vec_x,vec_y) as training
+ // sample. The density (rp) is computed as a by-product.
+ // Here, y = 1, so we should have rp = nw (un-normalized).
+ BOOST_FOREACH(val_& x,vec_x){
+ typedef binary_op::tuple_range<const vec_&,const vec_&> factory_;
+ typedef factory_::type range_tuple_;
+ range_tuple_ range_tuple = factory_::make(vec_x,vec_y);
+ nw_visitor_ nw_visitor = std::for_each(
+ begin(range_tuple),
+ end(range_tuple),
+ nw_visitor_tuple_(bandwidth,x)
+ ).nw_visitor();
+ val_ u_nw = nw_visitor.unnormalized_estimate();
+ vec_nw.push_back(u_nw);
+ rp_visitor_ rp_visitor = nw_visitor.rp_visitor();
+ val_ rp = rp_visitor.estimate();
+ BOOST_ASSERT(fabs(rp-u_nw)<eps);
+ }
+
+ typedef binary_op::tuple_range<const vec_&,const vec_&> factory_;
+ typedef factory_::type range_xy_;
+ range_xy_ range_xy = factory_::make(vec_x,vec_y);
+ // A pair of iterators is cheap to copy so no need to pass it by reference
+ typedef statistics::kernel::estimator<
+ range_xy_,
+ statistics::kernel::nw_visitor_tuple,
+ gauss_k_
+ > estimator_;
+ estimator_ estimator(bandwidth);
+ statistics::train(estimator,range_xy);
+
+ // Same as previous but calls estimator instead of for_each
+ BOOST_FOREACH(val_& x,vec_x){
+ // A local definition of nw_visitor_ is needed because x is passed
+ // by ref, not by value as in that outside the scope
+ typedef estimator_::result<val_>::type result_type;
+ typedef result_type::nw_visitor_type nw_visitor_;
+ typedef result_type::rp_visitor_type rp_visitor_;
+ nw_visitor_ nw_visitor = estimator(x).nw_visitor();
+ val_ u_nw = nw_visitor.unnormalized_estimate();
+ rp_visitor_ rp_visitor = nw_visitor.rp_visitor();
+ val_ rp = rp_visitor.estimate();
+ BOOST_ASSERT(fabs(rp-u_nw)<eps);
+ }
+
+ // Shorter version of the above
+ BOOST_FOREACH(val_& x,vec_x){
+ BOOST_AUTO( nw_visitor , estimator(x).nw_visitor() );
+ val_ u_nw = nw_visitor.unnormalized_estimate();
+ BOOST_AUTO( rp_visitor , nw_visitor.rp_visitor() );
+ val_ rp = rp_visitor.estimate();
+ BOOST_ASSERT(fabs(rp-u_nw)<eps);
+ }
+
+ out << "<-" << std::endl;
+}
\ No newline at end of file
Added: sandbox/statistics/kernel/libs/statistics/kernel/example/scalar_nw.h
==============================================================================
--- (empty file)
+++ sandbox/statistics/kernel/libs/statistics/kernel/example/scalar_nw.h 2009-08-28 18:46:30 EDT (Fri, 28 Aug 2009)
@@ -0,0 +1,15 @@
+///////////////////////////////////////////////////////////////////////////////
+// statistics::kernel::example::scalar_nw.h //
+// //
+// 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 LIBS_KERNEL_EXAMPLE_SCALAR_NW_H_ER_2009
+#define LIBS_KERNEL_EXAMPLE_SCALAR_NW_H_ER_2009
+#include <ostream>
+
+void example_scalar_nw(std::ostream&);
+
+#endif
+
Added: sandbox/statistics/kernel/libs/statistics/kernel/example/scalar_rp.cpp
==============================================================================
--- (empty file)
+++ sandbox/statistics/kernel/libs/statistics/kernel/example/scalar_rp.cpp 2009-08-28 18:46:30 EDT (Fri, 28 Aug 2009)
@@ -0,0 +1,90 @@
+///////////////////////////////////////////////////////////////////////////////
+// statistics::kernel::example::scalar_rp.cpp //
+// //
+// 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) //
+///////////////////////////////////////////////////////////////////////////////
+#include <vector>
+#include <algorithm>
+#include <iterator>
+#include <boost/range.hpp>
+#include <boost/foreach.hpp>
+#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/normal_distribution.hpp>
+#include <boost/random/variate_generator.hpp>
+#include <boost/math/special_functions/fpclassify.hpp> //needed?
+#include <boost/math/tools/precision.hpp>
+
+#include <boost/statistics/kernel/scalar/gaussian.hpp>
+#include <boost/statistics/kernel/functional/rp_visitor.hpp>
+#include <boost/statistics/kernel/functional/estimator.hpp>
+#include <libs/statistics/kernel/example/scalar_rp.h>
+
+void example_scalar_rp(std::ostream& out){
+ out << "-> example_scalar_rp : ";
+ using namespace boost;
+
+ //Types
+ typedef double val_;
+ typedef std::vector<val_> vec_;
+ typedef mt19937 urng_;
+ typedef normal_distribution<val_> norm_;
+ typedef variate_generator<urng_&,norm_> gen_;
+ typedef statistics::kernel::scalar::gaussian<val_> gauss_k_;
+ typedef statistics::kernel::rp_visitor<gauss_k_,val_> rp_visitor_;
+
+ // Contants
+ const val_ bandwidth = 0.5;
+ const val_ eps = math::tools::epsilon<val_>();
+ const unsigned n = 10;
+
+ // Initialization
+ vec_ vec_x(n);
+ vec_ vec_rp; vec_rp.reserve(n);
+ urng_ urng;
+ norm_ norm;
+ gen_ gen(urng,norm);
+ std::generate_n(
+ begin(vec_x),
+ n,
+ gen
+ );
+
+ // Computes a density estimate for each x in vec_x using vec_x as sample
+ BOOST_FOREACH(val_& x,vec_x){
+ val_ rp = for_each(
+ begin(vec_x),
+ end(vec_x),
+ rp_visitor_(bandwidth,x)
+ ).estimate();
+ vec_rp.push_back(rp);
+ }
+
+ std::copy(
+ begin(vec_rp),
+ end(vec_rp),
+ std::ostream_iterator<val_>(out," ")
+ );
+
+ typedef sub_range<vec_> sub_x_;
+ typedef
+ statistics::kernel::estimator<
+ sub_x_,
+ statistics::kernel::rp_visitor,
+ gauss_k_
+ > estimator_;
+ estimator_ estimator(bandwidth);
+ statistics::train(estimator,sub_x_(vec_x));
+ vec_ vec_rp2; vec_rp2.reserve(n);
+
+ // Same as previous but calls estimator instead of for_each
+ for(unsigned i = 0; i<n; i++){
+ val_ x = vec_x[i];
+ val_ rp = vec_rp[i];
+ val_ rp2 = estimator(x).estimate();
+ BOOST_ASSERT(fabs(rp-rp2)<eps);
+ }
+
+ out << "<-" << std::endl;
+}
\ No newline at end of file
Added: sandbox/statistics/kernel/libs/statistics/kernel/example/scalar_rp.h
==============================================================================
--- (empty file)
+++ sandbox/statistics/kernel/libs/statistics/kernel/example/scalar_rp.h 2009-08-28 18:46:30 EDT (Fri, 28 Aug 2009)
@@ -0,0 +1,15 @@
+///////////////////////////////////////////////////////////////////////////////
+// statistics::kernel::example::scalar_rp.h //
+// //
+// 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 LIBS_KERNEL_EXAMPLE_SCALAR_RP_H_ER_2009
+#define LIBS_KERNEL_EXAMPLE_SCALAR_RP_H_ER_2009
+#include <ostream>
+
+void example_scalar_rp(std::ostream&);
+
+#endif
+
Added: sandbox/statistics/kernel/libs/statistics/kernel/src/main.cpp
==============================================================================
--- (empty file)
+++ sandbox/statistics/kernel/libs/statistics/kernel/src/main.cpp 2009-08-28 18:46:30 EDT (Fri, 28 Aug 2009)
@@ -0,0 +1,22 @@
+///////////////////////////////////////////////////////////////////////////////
+// statistics::kernel::src::main.cpp //
+// //
+// 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) //
+///////////////////////////////////////////////////////////////////////////////
+#include <iostream>
+#include <libs/statistics/kernel/example/scalar_rp.h>
+#include <libs/statistics/kernel/example/scalar_nw.h>
+#include <libs/statistics/kernel/example/kernel_mono_rp.h>
+#include <libs/statistics/kernel/example/benchmark_scalar.h>
+
+int main(){
+
+ // example_scalar_rp(std::cout);
+ // example_scalar_nw(std::cout);
+ // example_kernel_mono_rp(std::cout);
+ example_benchmark_scalar(std::cout);
+
+ return 0;
+}
\ 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