|
|
Boost-Commit : |
Subject: [Boost-commit] svn:boost r55946 - in sandbox/statistics/survival_model: . boost boost/statistics boost/statistics/survival boost/statistics/survival/model boost/statistics/survival/model/algorithm boost/statistics/survival/model/meta boost/statistics/survival/model/models boost/statistics/survival/model/models/exponential boost/statistics/survival/model/models/exponential/detail libs libs/statistics libs/statistics/survival libs/statistics/survival/model libs/statistics/survival/model/doc libs/statistics/survival/model/example libs/statistics/survival/model/src
From: erwann.rogard_at_[hidden]
Date: 2009-08-31 21:58:23
Author: e_r
Date: 2009-08-31 21:58:21 EDT (Mon, 31 Aug 2009)
New Revision: 55946
URL: http://svn.boost.org/trac/boost/changeset/55946
Log:
re-adding survival_model, but preceded by dir statistics
Added:
sandbox/statistics/survival_model/
sandbox/statistics/survival_model/boost/
sandbox/statistics/survival_model/boost/statistics/
sandbox/statistics/survival_model/boost/statistics/survival/
sandbox/statistics/survival_model/boost/statistics/survival/model/
sandbox/statistics/survival_model/boost/statistics/survival/model/algorithm/
sandbox/statistics/survival_model/boost/statistics/survival/model/algorithm/include.hpp (contents, props changed)
sandbox/statistics/survival_model/boost/statistics/survival/model/algorithm/prepare_weights.hpp (contents, props changed)
sandbox/statistics/survival_model/boost/statistics/survival/model/include.hpp (contents, props changed)
sandbox/statistics/survival_model/boost/statistics/survival/model/meta/
sandbox/statistics/survival_model/boost/statistics/survival/model/meta/data.hpp (contents, props changed)
sandbox/statistics/survival_model/boost/statistics/survival/model/meta/include.hpp (contents, props changed)
sandbox/statistics/survival_model/boost/statistics/survival/model/meta/model_data.hpp (contents, props changed)
sandbox/statistics/survival_model/boost/statistics/survival/model/meta/response.hpp (contents, props changed)
sandbox/statistics/survival_model/boost/statistics/survival/model/models/
sandbox/statistics/survival_model/boost/statistics/survival/model/models/exponential/
sandbox/statistics/survival_model/boost/statistics/survival/model/models/exponential/detail/
sandbox/statistics/survival_model/boost/statistics/survival/model/models/exponential/detail/log_likelihood.hpp (contents, props changed)
sandbox/statistics/survival_model/boost/statistics/survival/model/models/exponential/include.hpp (contents, props changed)
sandbox/statistics/survival_model/boost/statistics/survival/model/models/exponential/log_likelihood.hpp (contents, props changed)
sandbox/statistics/survival_model/boost/statistics/survival/model/models/exponential/meta_failure_distribution.hpp (contents, props changed)
sandbox/statistics/survival_model/boost/statistics/survival/model/models/exponential/model.hpp (contents, props changed)
sandbox/statistics/survival_model/libs/
sandbox/statistics/survival_model/libs/statistics/
sandbox/statistics/survival_model/libs/statistics/survival/
sandbox/statistics/survival_model/libs/statistics/survival/model/
sandbox/statistics/survival_model/libs/statistics/survival/model/doc/
sandbox/statistics/survival_model/libs/statistics/survival/model/doc/readme.txt (contents, props changed)
sandbox/statistics/survival_model/libs/statistics/survival/model/example/
sandbox/statistics/survival_model/libs/statistics/survival/model/example/exponential.cpp (contents, props changed)
sandbox/statistics/survival_model/libs/statistics/survival/model/example/exponential.h (contents, props changed)
sandbox/statistics/survival_model/libs/statistics/survival/model/example/posterior_analysis.cpp (contents, props changed)
sandbox/statistics/survival_model/libs/statistics/survival/model/example/posterior_analysis.h (contents, props changed)
sandbox/statistics/survival_model/libs/statistics/survival/model/src/
sandbox/statistics/survival_model/libs/statistics/survival/model/src/main.cpp (contents, props changed)
Added: sandbox/statistics/survival_model/boost/statistics/survival/model/algorithm/include.hpp
==============================================================================
--- (empty file)
+++ sandbox/statistics/survival_model/boost/statistics/survival/model/algorithm/include.hpp 2009-08-31 21:58:21 EDT (Mon, 31 Aug 2009)
@@ -0,0 +1,14 @@
+///////////////////////////////////////////////////////////////////////////////
+// statistics::survival::model::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_SURVIVAL_MODEL_FUNCTIONAL_INCLUDE_HPP_ER_2009
+#define BOOST_STATISTICS_SURVIVAL_MODEL_FUNCTIONAL_INCLUDE_HPP_ER_2009
+
+#include <boost/statistics/survival/model/functional/experiment_simulator.hpp>
+#include <boost/statistics/survival/model/functional/prepare_weights.hpp>
+
+#endif
\ No newline at end of file
Added: sandbox/statistics/survival_model/boost/statistics/survival/model/algorithm/prepare_weights.hpp
==============================================================================
--- (empty file)
+++ sandbox/statistics/survival_model/boost/statistics/survival/model/algorithm/prepare_weights.hpp 2009-08-31 21:58:21 EDT (Mon, 31 Aug 2009)
@@ -0,0 +1,88 @@
+///////////////////////////////////////////////////////////////////////////////
+// statistics::survival::model::algorithm::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_SURVIVAL_MODEL_FUNCTIONAL_PREPARE_WEIGHTS_HPP_ER_2009
+#define BOOST_STATISTICS_SURVIVAL_MODEL_FUNCTIONAL_PREPARE_WEIGHTS_HPP_ER_2009
+#include <ostream>
+#include <boost/importance_sampling/detail/prepare_weights_impl.hpp>
+#include <boost/model/wrap/aggregate/prior_model_dataset.hpp>
+#include <boost/model/algorithm/log_posteriors.hpp>
+
+namespace boost{
+namespace survival{
+namespace model{
+
+ // TODO is this deprecated?
+ //
+ // Maps proposal values and their log_pdfs to importance sampling weights.
+ //
+ // Warning: Elements in [b_p,e_p) will be rearranged
+ template<typename T>
+ class prepare_weights : public is::prepare_weights_impl<T> {
+ typedef is::prepare_weights_impl<T> super_;
+ public:
+ typedef T value_type;
+
+ prepare_weights();
+ prepare_weights(const T& max_log);
+
+ template<
+ typename D,typename M,typename Rx,typename Re,
+ typename ItP, typename ItLw
+ >
+ void operator()(
+ boost::model::prior_model_dataset_<D,M,Rx,Re> pmd,
+ ItP b_p, // Proposal values
+ ItP e_p,
+ ItLw b_lw, // Proposal log pdfs
+ ItLw o_w // Outputs weights - not back_inserter
+ );
+ };
+
+ // Implementation
+
+ template<typename T>
+ prepare_weights<T>::prepare_weights():super_(){}
+
+ template<typename T>
+ prepare_weights<T>::prepare_weights(const T& max_log):super_(max_log){}
+
+ template<typename T>
+ template<
+ typename D,typename M,typename Rx,typename Re,
+ typename ItP, typename ItLw
+ >
+ void
+ prepare_weights<T>::operator()(
+ boost::model::prior_model_dataset_<D,M,Rx,Re> pmd,
+ ItP b_p, // Proposal values
+ ItP e_p,
+ ItLw b_lw, // Proposal log pdfs
+ ItLw o_w // Outputs weights
+ ){
+ typedef boost::model::prior_model_dataset_<D,M,Rx,Re> pmd_;
+ ItLw b_w = o_w;
+ ItLw e_w = boost::model::log_posteriors<value_type>(
+ pmd,
+ b_p,
+ e_p,
+ b_lw,
+ o_w
+ );
+ return super_::impl(
+ b_w,
+ e_w,
+ b_p
+ ); // This affects the ordering of [b_w,e_w) and [b_p,e_p)
+ }
+
+
+}// model
+}// survival
+}// boost
+
+#endif
\ No newline at end of file
Added: sandbox/statistics/survival_model/boost/statistics/survival/model/include.hpp
==============================================================================
--- (empty file)
+++ sandbox/statistics/survival_model/boost/statistics/survival/model/include.hpp 2009-08-31 21:58:21 EDT (Mon, 31 Aug 2009)
@@ -0,0 +1,14 @@
+///////////////////////////////////////////////////////////////////////////////
+// statistics::survival::model::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_SURVIVAL_MODEL_INCLUDE_HPP_ER_2009
+#define BOOST_STATISTICS_SURVIVAL_MODEL_INCLUDE_HPP_ER_2009
+
+#include <boost/statistics/survival/model/meta/wrap.hpp>
+#include <boost/statistics/survival/model/functional/prepare_weights.hpp>
+
+#endif
\ No newline at end of file
Added: sandbox/statistics/survival_model/boost/statistics/survival/model/meta/data.hpp
==============================================================================
--- (empty file)
+++ sandbox/statistics/survival_model/boost/statistics/survival/model/meta/data.hpp 2009-08-31 21:58:21 EDT (Mon, 31 Aug 2009)
@@ -0,0 +1,32 @@
+///////////////////////////////////////////////////////////////////////////////
+// statistics::survival::model::meta::data.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_SURVIVAL_MODEL_META_DATA_HPP_ER_2009
+#define BOOST_STATISTICS_SURVIVAL_MODEL_META_DATA_HPP_ER_2009
+#include <boost/statistics/model/wrap/aggregate/data.hpp>
+#include <boost/statistics/survival/data/data/event.hpp>
+
+namespace boost{
+namespace statistics{
+namespace survival{
+namespace model{
+namespace meta{
+
+ // See statistics/model/libs/doc/readme
+ template<typename T,typename X>
+ struct data{
+ typedef boost::statistics::survival::data::event<T> y_;
+ typedef boost::statistics::model::data_< X, y_ > type;
+ };
+
+}// meta
+}// model
+}// survival
+}// statistics
+}// boost
+
+#endif
\ No newline at end of file
Added: sandbox/statistics/survival_model/boost/statistics/survival/model/meta/include.hpp
==============================================================================
--- (empty file)
+++ sandbox/statistics/survival_model/boost/statistics/survival/model/meta/include.hpp 2009-08-31 21:58:21 EDT (Mon, 31 Aug 2009)
@@ -0,0 +1,15 @@
+///////////////////////////////////////////////////////////////////////////////
+// statistics::survival::model::meta::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_SURVIVAL_MODEL_META_INCLUDE_HPP_ER_2009
+#define BOOST_STATISTICS_SURVIVAL_MODEL_META_INCLUDE_HPP_ER_2009
+
+#include <boost/statistics/survival/model/meta/data.hpp>
+#include <boost/statistics/survival/model/meta/response.hpp>
+#include <boost/statistics/survival/model/meta/model_data.hpp>
+
+#endif
\ No newline at end of file
Added: sandbox/statistics/survival_model/boost/statistics/survival/model/meta/model_data.hpp
==============================================================================
--- (empty file)
+++ sandbox/statistics/survival_model/boost/statistics/survival/model/meta/model_data.hpp 2009-08-31 21:58:21 EDT (Mon, 31 Aug 2009)
@@ -0,0 +1,32 @@
+///////////////////////////////////////////////////////////////////////////////
+// statistics::survival::model::meta::model_data.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_SURVIVAL_MODEL_META_MODEL_DATA_HPP_ER_2009
+#define BOOST_STATISTICS_SURVIVAL_MODEL_META_MODEL_DATA_HPP_ER_2009
+#include <boost/statistics/model/wrap/aggregate/model_data.hpp>
+#include <boost/statistics/survival/data/data/event.hpp>
+
+namespace boost{
+namespace statistics{
+namespace survival{
+namespace model{
+namespace meta{
+
+ // See statistics/model/libs/doc/readme
+ template<typename T,typename M,typename X>
+ struct model_data{
+ typedef boost::statistics::survival::data::event<T> y_;
+ typedef boost::statistics::model::model_data_< M, X, y_ > type;
+ };
+
+}// meta
+}// model
+}// survival
+}// statistics
+}// boost
+
+#endif
\ No newline at end of file
Added: sandbox/statistics/survival_model/boost/statistics/survival/model/meta/response.hpp
==============================================================================
--- (empty file)
+++ sandbox/statistics/survival_model/boost/statistics/survival/model/meta/response.hpp 2009-08-31 21:58:21 EDT (Mon, 31 Aug 2009)
@@ -0,0 +1,32 @@
+///////////////////////////////////////////////////////////////////////////////
+// statistics::survival::model::meta::response.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_SURVIVAL_MODEL_META_RESPONSE_HPP_ER_2009
+#define BOOST_STATISTICS_SURVIVAL_MODEL_META_RESPONSE_HPP_ER_2009
+#include <boost/statistics/model/wrap/unary/response.hpp>
+#include <boost/statistics/survival/data/data/event.hpp>
+
+namespace boost{
+namespace statistics{
+namespace survival{
+namespace model{
+namespace meta{
+
+ // See statistics/model/libs/doc/readme
+ template<typename T>
+ struct response{
+ typedef survival::data::event<T> y_;
+ typedef boost::statistics::model::response_<y_> type;
+ };
+
+}// meta
+}// model
+}// survival
+}// statistics
+}// boost
+
+#endif
\ No newline at end of file
Added: sandbox/statistics/survival_model/boost/statistics/survival/model/models/exponential/detail/log_likelihood.hpp
==============================================================================
--- (empty file)
+++ sandbox/statistics/survival_model/boost/statistics/survival/model/models/exponential/detail/log_likelihood.hpp 2009-08-31 21:58:21 EDT (Mon, 31 Aug 2009)
@@ -0,0 +1,62 @@
+//////////////////////////////////////////////////////////////////////////////////
+// statistics::survival::model::models::exponential::detail::log_likelihood.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_SURVIVAL_MODEL_MODELS_EXPONENTIAL_DETAIL_LOG_LIKELIHOOD_HPP_ER_2009
+#define BOOST_STATISTICS_SURVIVAL_MODEL_MODELS_EXPONENTIAL_DETAIL_LOG_LIKELIHOOD_HPP_ER_2009
+#include <cmath>
+#include <stdexcept>
+#include <boost/format.hpp>
+#include <boost/statistics/survival/data/data/event.hpp>
+
+namespace boost{
+namespace statistics{
+
+namespace survival{
+namespace model{
+namespace exponential{
+namespace detail{
+
+ template<typename T>
+ T
+ log_likelihood(
+ const T& log_rate,
+ const data::event<T>& e
+ ){
+ static const char* msg =
+ "survival::model::exponential::log_unnromalized_pdf(%1%,%2%)";
+ typedef T value_type;
+ value_type result = exp(log_rate);
+ result *= (- e.time());
+ try{
+ if( boost::math::isinf(result) ){
+ throw std::runtime_error("isinf(result)");
+ }
+ if( boost::math::isnan(result) ){
+ throw std::runtime_error("isnan(result)");
+ }
+ }catch(std::exception ex){
+ std::string str = msg;
+ str += ex.what();
+ format f(str); f % log_rate % e;
+ throw std::runtime_error(
+ f.str()
+ );
+ }
+ if(e.failure()){
+ result += log_rate;
+ }
+ return result;
+ }
+
+}// detail
+}// exponential
+}// model
+}// survival
+}// statistics
+}// boost
+
+#endif
\ No newline at end of file
Added: sandbox/statistics/survival_model/boost/statistics/survival/model/models/exponential/include.hpp
==============================================================================
--- (empty file)
+++ sandbox/statistics/survival_model/boost/statistics/survival/model/models/exponential/include.hpp 2009-08-31 21:58:21 EDT (Mon, 31 Aug 2009)
@@ -0,0 +1,15 @@
+///////////////////////////////////////////////////////////////////////////////
+// statistics::survival::model::models::exponential::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_SURVIVAL_MODEL_MODELS_EXPONENTIAL_INCLUDE_HPP_ER_2009
+#define BOOST_STATISTICS_SURVIVAL_MODEL_MODELS_EXPONENTIAL_INCLUDE_HPP_ER_2009
+
+#include <boost/statistics/survival/model/models/exponential/model.hpp>
+#include <boost/statistics/survival/model/models/exponential/log_likelihood.hpp>
+#include <boost/statistics/survival/model/models/exponential/meta_failure_distribution.hpp>
+
+#endif
\ No newline at end of file
Added: sandbox/statistics/survival_model/boost/statistics/survival/model/models/exponential/log_likelihood.hpp
==============================================================================
--- (empty file)
+++ sandbox/statistics/survival_model/boost/statistics/survival/model/models/exponential/log_likelihood.hpp 2009-08-31 21:58:21 EDT (Mon, 31 Aug 2009)
@@ -0,0 +1,63 @@
+///////////////////////////////////////////////////////////////////////////////
+// statistics::survival::model::models::exponential::log_likelihood.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_SURVIVAL_MODEL_MODELS_EXPONENTIAL_LOG_LIKELIHOOD_HPP_ER_2009
+#define BOOST_STATISTICS_SURVIVAL_MODEL_MODELS_EXPONENTIAL_LOG_LIKELIHOOD_HPP_ER_2009
+#include <boost/statistics/survival/model/meta/model_data.hpp>
+#include <boost/statistics/survival/model/models/exponential/model.hpp>
+#include <boost/statistics/survival/model/models/exponential/detail/log_likelihood.hpp>
+
+namespace boost{
+namespace statistics{
+namespace model{
+
+ // Models Model (sandbox/statistics/model).
+ //
+ // Intentionally not in namespace survival
+ template<typename T,typename X,typename B>
+ T log_likelihood(
+ typename survival::model::meta::model_data<
+ T,
+ survival::model::exponential::model<T>,
+ X
+ >::type md,
+ const B& beta
+ ){
+ typedef survival::model::exponential::model<T> model_;
+
+ T lr = model_::log_rate(
+ md.covariate(),
+ beta
+ );
+ return survival::model::exponential::detail::log_likelihood(
+ lr,
+ md.response()
+ );
+ }
+
+ // If B == X we need this overload, or else the compiler cannot find
+ // the above definition
+ template<typename T,typename X>
+ T log_likelihood(
+ typename survival::model::meta::model_data<
+ T,
+ survival::model::exponential::model<T>,
+ X
+ >::type md,
+ const X& beta
+ ){
+ return log_likelihood<T,X,X>(
+ md,
+ beta
+ );
+ }
+
+}// model
+}// statistics
+}// boost
+
+#endif
\ No newline at end of file
Added: sandbox/statistics/survival_model/boost/statistics/survival/model/models/exponential/meta_failure_distribution.hpp
==============================================================================
--- (empty file)
+++ sandbox/statistics/survival_model/boost/statistics/survival/model/models/exponential/meta_failure_distribution.hpp 2009-08-31 21:58:21 EDT (Mon, 31 Aug 2009)
@@ -0,0 +1,49 @@
+///////////////////////////////////////////////////////////////////////////////
+// statistics::survival::model::models::meta_failure_distribution.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_SURVIVAL_MODEL_MODELS_META_FAILURE_DISTRIBUTION_HPP_ER_2009
+#define BOOST_STATISTICS_SURVIVAL_MODEL_MODELS_META_FAILURE_DISTRIBUTION_HPP_ER_2009
+#include <cmath>
+#include <boost/statistics/model/wrap/aggregate/model_covariate_parameter.hpp>
+#include <boost/standard_distribution/distributions/exponential.hpp>
+#include <boost/statistics/survival/data/meta/failure_distribution.hpp>
+#include <boost/statistics/survival/model/models/exponential/model.hpp>
+
+namespace boost{
+namespace statistics{
+
+namespace survival{
+namespace data{
+
+ template<typename T>
+ struct meta_failure_distribution< survival::model::exponential::model<T> >{
+ typedef survival::model::exponential::model<T> model_;
+ typedef math::exponential_distribution<T> type;
+
+ template<typename X,typename P>
+ static type make(
+ boost::statistics::model::model_covariate_parameter_<model_,X,P> mcp
+ ){
+ T lambda = model_::log_rate(
+ mcp.covariate(),
+ mcp.parameter()
+ );
+ lambda = exp( lambda );
+ return type(
+ lambda
+ );
+ }
+ };
+
+
+
+}// data
+}// survival
+}// statistics
+}// boost
+
+#endif
\ No newline at end of file
Added: sandbox/statistics/survival_model/boost/statistics/survival/model/models/exponential/model.hpp
==============================================================================
--- (empty file)
+++ sandbox/statistics/survival_model/boost/statistics/survival/model/models/exponential/model.hpp 2009-08-31 21:58:21 EDT (Mon, 31 Aug 2009)
@@ -0,0 +1,83 @@
+///////////////////////////////////////////////////////////////////////////////
+// statistics::survival::model::exponential::model.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_SURVIVAL_MODEL_MODELS_EXPONENTIAL_MODEL_HPP_ER_2009
+#define BOOST_STATISTICS_SURVIVAL_MODEL_MODELS_EXPONENTIAL_MODEL_HPP_ER_2009
+#include <numeric>
+#include <boost/mpl/assert.hpp>
+#include <boost/archive/text_oarchive.hpp>
+#include <boost/archive/text_iarchive.hpp>
+#include <boost/type_traits/is_scalar.hpp>
+#include <boost/utility/enable_if.hpp>
+#include <boost/range.hpp>
+
+namespace boost{
+namespace statistics{
+
+namespace survival{
+namespace model{
+namespace exponential{
+
+ template<typename T>
+ class model{
+ public:
+ typedef T value_type;
+
+ model();
+
+ template<typename X,typename B>
+ static typename boost::enable_if<boost::is_scalar<X>,T>::type
+ log_rate(const X& x,const B& b);
+
+ template<typename X,typename B>
+ static typename boost::disable_if<boost::is_scalar<X>,T>::type
+ log_rate(const X& x,const B& b);
+
+ protected:
+ friend class boost::serialization::access;
+ template<class Archive>
+ void serialize(Archive & ar, const unsigned int version){
+ // no member variables
+ }
+
+ };
+
+ // Implementation //
+
+ template<typename T>
+ model<T>::model(){}
+
+ template<typename T>
+ template<typename X,typename B>
+ typename boost::enable_if<boost::is_scalar<X>,T>::type
+ model<T>::log_rate(const X& x,const B& b){
+ BOOST_MPL_ASSERT((
+ is_scalar<B>
+ ));
+ return ( static_cast<T>( x ) * static_cast<T>( b ) );
+ }
+
+ template<typename T>
+ template<typename X,typename B>
+ typename boost::disable_if<is_scalar<X>,T>::type
+ model<T>::log_rate(const X& x,const B& b){
+ BOOST_ASSERT(size(x) == size(b));
+ return std::inner_product(
+ boost::begin(x),
+ boost::end(x),
+ boost::begin(b),
+ static_cast<T>(0)
+ );
+ }
+
+}// exponential
+}// model
+}// survival
+}// statistics
+}// boost
+
+#endif
Added: sandbox/statistics/survival_model/libs/statistics/survival/model/doc/readme.txt
==============================================================================
--- (empty file)
+++ sandbox/statistics/survival_model/libs/statistics/survival/model/doc/readme.txt 2009-08-31 21:58:21 EDT (Mon, 31 Aug 2009)
@@ -0,0 +1,116 @@
+//////////////////////////////////////////////////////////////////////////////
+// statistics::survival::model::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 ]
+
+ In this C++ library, we provide models of Model (statistics/model) in the
+ survival modeling framework.
+
+[ Notation ]
+
+ See statistics/model/libs/model/doc/readme and the concept Model. Let mcp
+ denote an instance of model_covariate_parameter_<M, X, P>
+
+[ Bug ]
+
+ See libs/statistics/survival/example/posterior_analysis.cpp
+
+[ Compiler ]
+
+gcc version i686-apple-darwin9-gcc-4.0.1 (GCC) 4.0.1
+
+[ Dependencies ]
+
+/boost_1_39_0/
+/sandbox/statistics/arithmetic/
+/sandbox/statistics/non_par/
+/sandbox/statistics/importance_sampling/
+/sandbox/statistics/importance_weights/
+/sandbox/statistics/joint_dist/
+/sandbox/statistics/empirical_cdf/
+/sandbox/statistics/mpl/
+/sandbox/statistics/random/
+/sandbox/statistics/binary_op/
+/sandbox/statistics/functional/
+/sandbox/statistics/arithmetic/
+/sandbox/statistics/iterator/
+/sandbox/statistics/dist_random/
+/sandbox/statistics/scalar_dist/
+/sandbox/statistics/standard_distribution/
+/sandbox/statistics/survival_data/
+/sandbox/statistics/model/
+
+Link to : libboost_serialization-xgcc42-mt-1_39.a
+
+[ History ]
+
+July 2009 : Current version
+
+[ models ]
+
+These are models of [ SurvivalModel ] defined as:
+
+Let U/u denote a type/object modeling UniformaRandomNumberGenerator (see
+Boost.Random)
+
+M refines SurvivalModel if it models Model with the requirement that
+ Expression Returns
+ random::failure_time<T>(mcp,urng) Object of type T
+
+Modeled by :
+ model::exponential::model
+
+[ meta ]
+
+ - meta::failure_distribution maps M to a distribution in math/distributions
+ - meta::failure_random maps M to a RandomDistribution in boost/random
+
+ If these mappings are specified for M, random::failure_time is automatically
+ defined for M.
+
+[ Sources ]
+
+[1] Bayesian survival analysis By Joseph George Ibrahim, Ming-Hui Chen,
+Debajyoti Sinha
+
+[ Output ]
+-> example_model_exponential : size(x_cycle) = 100
+(-2.48172,-1.87264)
+log(prior,likelihood,posterior)
+(-2,-11.3373,-13.3373)
+(-0.5,-7.28159,-7.78159)
+(-0,-5.08452,-5.08452)
+(-0.5,-4.18529,-4.68529)
+(-2,-4.35436,-6.35436)
+
+log(prior,prior2)
+(-2,-2)
+(-0.5,-0.5)
+(-0,0)
+(-0.5,-0.5)
+(-2,-2)
+<-
+
+-> example_posterior_analysis :
+i = 0, t = 0.169766, pws = (222.593,1,0.047237,0.5896)
+i = 1000, t = 0.167177, pws = (195.97,1,0.0261956,0.7733)
+i = 2000, t = 0.166789, pws = (190.525,1,0.0163538,0.836)
+i = 3000, t = 0.167136, pws = (197.041,1,0.0519114,0.5842)
+i = 4000, t = 0.168067, pws = (195.495,1,0.05305,0.5729)
+i = 5000, t = 0.171331, pws = (199.301,1,0.0381148,0.6841)
+i = 6000, t = 0.167311, pws = (205.46,1,0.0556848,0.5621)
+i = 7000, t = 0.167732, pws = (182.455,1,0.0639737,0.5192)
+i = 8000, t = 0.166502, pws = (181.937,1,0.0661214,0.5163)
+i = 9000, t = 0.16732, pws = (188.906,1,0.0135898,0.8683)<-
+
+
Added: sandbox/statistics/survival_model/libs/statistics/survival/model/example/exponential.cpp
==============================================================================
--- (empty file)
+++ sandbox/statistics/survival_model/libs/statistics/survival/model/example/exponential.cpp 2009-08-31 21:58:21 EDT (Mon, 31 Aug 2009)
@@ -0,0 +1,258 @@
+///////////////////////////////////////////////////////////////////////////////
+// statistics::survival::model::example::model::exponential.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 <stdexcept>
+#include <string> //needed?
+#include <vector>
+#include <limits>
+#include <ostream>
+#include <fstream>
+#include <algorithm>
+#include <iterator>
+// #include <boost/archive/binary_oarchive.hpp>
+// #include <boost/archive/binary_iarchive.hpp>
+#include <boost/archive/text_oarchive.hpp>
+#include <boost/archive/text_iarchive.hpp>
+#include <boost/serialization/vector.hpp>
+
+#include <boost/arithmetic/equal.hpp>
+#include <boost/assign/std/vector.hpp>
+#include <boost/iterator/range_cycle.hpp>
+#include <boost/range.hpp>
+#include <boost/assert.hpp>
+#include <boost/foreach.hpp>
+#include <boost/assign/std/vector.hpp>
+#include <boost/iterator/range_cycle.hpp>
+
+#include <boost/standard_distribution/distributions/normal.hpp>
+
+#include <boost/statistics/survival/data/include.hpp>
+#include <boost/statistics/survival/model/models/exponential/include.hpp>
+#include <libs/statistics/survival/model/example/exponential.h>
+
+// Must come after the model to be used
+#include <boost/statistics/model/include.hpp>
+
+#include <libs/statistics/survival/model/example/exponential.h>
+
+void example_exponential(std::ostream& out){
+
+ out << "-> example_model_exponential : ";
+
+ // Steps shown in this example:
+ //
+ // Loads the first batch of a set of records
+ // Creates events at given time
+ // Evaluates the likelihoods and posteriors
+
+ using namespace boost;
+ using namespace statistics;
+ namespace surv = survival;
+
+ // [ Types ]
+ // Value
+ typedef double val_;
+ typedef std::vector<val_> vals_;
+ typedef surv::constant<val_> const_;
+
+ // I/O
+ typedef boost::archive::text_oarchive oa_;
+ typedef boost::archive::text_iarchive ia_;
+
+ // Records
+ typedef surv::data::record<val_> record_;
+ typedef std::vector<record_> records_;
+ typedef range_iterator<records_>::type it_record_;
+
+ // Events
+ typedef surv::data::event<val_> event_;
+ typedef std::vector<event_> events_;
+ typedef range_iterator<events_>::type it_event_;
+
+ // Covariates
+ typedef val_ x_;
+ typedef vals_ r_x_;
+ typedef range_cycle<> range_cycle_;
+ typedef range_cycle_::apply<r_x_>::type x_cycle_;
+
+ // Model
+ typedef surv::model::exponential::model<val_> model_;
+ typedef val_ par_;
+ typedef vals_ pars_;
+
+ // [ Constants ]
+ const val_ entry_bound = const_::inf_;
+ const val_ par = 2.0;
+ const char* batches_path
+ = "/Users/erwann/projets/2009/Xcode/survival/build/Release/batches";
+
+ // [ Variables ]
+ long n_record;
+
+ // [ Upload first batch of records ]
+ records_ records;
+ {
+ std::ifstream ifs(batches_path);
+ if(ifs.good()){
+ ia_ ia(ifs);
+ ia >> records;
+ }else{
+ std::string str = "error opening : ";
+ str.append( batches_path );
+ throw std::runtime_error(str);
+ }
+ ifs.close();
+ }
+ n_record = boost::size( records );
+
+ // [ Events ]
+ events_ events;
+ events.reserve( size(records) );
+ surv::data::events(
+ begin(records),
+ end(records),
+ entry_bound,
+ std::back_inserter(events)
+ );
+
+ // [ Covariates ]
+ r_x_ r_x;
+ {
+ using namespace boost::assign;
+ r_x += -0.5, 0.5;
+ }
+ x_cycle_ x_cycle = range_cycle_::make(r_x,0,n_record);
+ out << "size(x_cycle) = " << size(x_cycle) << std::endl;
+ BOOST_ASSERT( size(x_cycle)>=size(events) );
+ // Resize x_cycle to a size that matches that of events
+ x_cycle.advance_end( - (size(x_cycle) - size(events)) );
+ BOOST_ASSERT( size(x_cycle) == size(events) );
+
+ // Model
+ model_ model;
+
+ // Pars
+ pars_ pars;
+ {
+ using namespace assign;
+ pars += -2.0, -1.0, 0.0, 1.0, 2.0;
+ }
+
+ // [ Likelihood ]
+
+ typedef math::normal_distribution<val_> mprior_;
+ mprior_ mprior;
+
+ out << '(';
+ out << model::log_likelihood<val_>(
+ model::make_model_data(
+ model,
+ r_x[0],
+ events[0]
+ ),
+ par
+ );
+ out << ',';
+ out << model::log_likelihood<val_>(
+ model::make_model_data(
+ model,
+ r_x[1],
+ events[1]
+ ),
+ par
+ ) << ')';
+
+
+ // [ Likelihoods ]
+ vals_ lls;
+ model::log_likelihoods<val_>(
+ model::make_model_dataset(model,r_x,events),
+ boost::begin(pars),
+ boost::end(pars),
+ std::back_inserter(lls)
+ );
+
+ // [ Prior ]
+ vals_ lprs;
+ math::transform<math::fun_wrap::log_unnormalized_pdf_>(
+ mprior,
+ boost::begin(pars),
+ boost::end(pars),
+ std::back_inserter(lprs)
+ );
+
+ // [ Posteriors ]
+
+ vals_ lpos;
+ model::log_posteriors<val_>(
+ model::make_prior_model_dataset(mprior,model,r_x,events),
+ boost::begin(pars),
+ boost::end(pars),
+ std::back_inserter(lpos)
+ );
+
+ // Consistency check
+ typedef range_iterator<vals_>::type it_val_;
+ {
+
+ it_val_ i_lpr = boost::begin(lprs);
+ it_val_ i_lpo = boost::begin(lpos);
+ out << std::endl;
+ out << "log(prior,likelihood,posterior)" << std::endl;
+ for(
+ it_val_ i_ll = boost::begin(lls);
+ i_ll< boost::end(lls);
+ i_ll++,i_lpr++,i_lpo++
+ ){
+ out << '(';
+ val_ lpr = *i_lpr; out << lpr << ',';
+ val_ ll = *i_ll; out << ll << ',';
+ val_ lpo = *i_lpo; out << lpo << ')' << std::endl;
+ val_ lpo2 = lpr + ll;
+ BOOST_ASSERT(
+ arithmetic_tools::equal(
+ lpo,
+ lpo2
+ )
+ );
+ }
+ }
+
+ // Consistency check2
+
+ {
+ vals_ lpr2s ( size(pars) );
+ model::log_posteriors<val_>(
+ model::make_prior_model_dataset(mprior,model,r_x,events),
+ boost::begin(pars),
+ boost::end(pars),
+ boost::begin(lls), //subtracted
+ boost::begin(lpr2s)
+ );
+
+ it_val_ i_lpr = boost::begin(lprs);
+ out << std::endl;
+ out << "log(prior,prior2)" << std::endl;
+ for(
+ it_val_ i_lpr2 = boost::begin(lpr2s);
+ i_lpr2< boost::end(lpr2s);
+ i_lpr2++,i_lpr++
+ ){
+ out << '(';
+ val_ lpr = *i_lpr; out << lpr << ',';
+ val_ lpr2 = *i_lpr2; out << lpr2 << ')' << std::endl;
+ BOOST_ASSERT(
+ arithmetic_tools::equal(
+ lpr,
+ lpr2
+ )
+ );
+ }
+ }
+
+ out << "<-" << std::endl;
+}
Added: sandbox/statistics/survival_model/libs/statistics/survival/model/example/exponential.h
==============================================================================
--- (empty file)
+++ sandbox/statistics/survival_model/libs/statistics/survival/model/example/exponential.h 2009-08-31 21:58:21 EDT (Mon, 31 Aug 2009)
@@ -0,0 +1,14 @@
+///////////////////////////////////////////////////////////////////////////////
+// statistics::survival::model::example::model::exponential.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_SURVIVAL_EXAMPLE_MODEL_EXPONENTIAL_HPP_ER_2009
+#define LIBS_SURVIVAL_EXAMPLE_MODEL_EXPONENTIAL_HPP_ER_2009
+#include <ostream>
+
+void example_exponential(std::ostream&);
+
+#endif
\ No newline at end of file
Added: sandbox/statistics/survival_model/libs/statistics/survival/model/example/posterior_analysis.cpp
==============================================================================
--- (empty file)
+++ sandbox/statistics/survival_model/libs/statistics/survival/model/example/posterior_analysis.cpp 2009-08-31 21:58:21 EDT (Mon, 31 Aug 2009)
@@ -0,0 +1,390 @@
+///////////////////////////////////////////////////////////////////////////////
+// statistics::survival::model::example::model::posterior_analysis.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 <ostream>
+#include <ostream>
+#include <fstream>
+#include <stdexcept>
+#include <string> //needed?
+#include <vector>
+#include <limits>
+#include <algorithm>
+#include <iterator>
+// #include <boost/archive/binary_oarchive.hpp>
+// #include <boost/archive/binary_iarchive.hpp>
+#include <boost/archive/text_oarchive.hpp>
+#include <boost/archive/text_iarchive.hpp>
+#include <boost/serialization/vector.hpp>
+
+#include <boost/timer.hpp>
+#include <boost/assign/std/vector.hpp>
+#include <boost/iterator/range_cycle.hpp>
+#include <boost/range.hpp>
+#include <boost/assert.hpp>
+#include <boost/foreach.hpp>
+
+#include <boost/standard_distribution/distributions/normal.hpp>
+#include <boost/math/distributions/uniform.hpp>
+#include <boost/statistics/empirical_cdf/algorithm/sequential_kolmogorov_smirnov_distance.hpp>
+#include <boost/statistics/survival/data/include.hpp>
+#include <boost/statistics/survival/model/models/exponential/include.hpp>
+
+// Must come after the model to be used
+#include <boost/statistics/model/include.hpp>
+
+#include <boost/importance_weights/algorithm/prepare_weights.hpp>
+#include <boost/importance_sampling/include.hpp>
+
+#include <libs/statistics/survival/model/example/posterior_analysis.h>
+
+void example_posterior_analysis(std::ostream& out){
+
+ out << "-> example_posterior_analysis : ";
+ out.flush();
+
+ // Steps shown in this example:
+ //
+ // Samples proposal draws
+ // Successively loads batches of records, and for each
+ // Creates events at a given time
+ // Importance samples from the posterior as the target density
+ // Saves the Cook-Gelman (cg) statistics for each.
+ // Checks the cg data against U[0,1] using Kolmogorov Smirnov
+
+ using namespace boost;
+ using namespace statistics;
+ namespace surv = survival;
+ namespace iw = importance_weights;
+
+ // Types
+ typedef std::string str_;
+ typedef double val_;
+ typedef std::vector<val_> vals_;
+ typedef surv::constant<val_> const_;
+
+ // [ Input ]
+ typedef boost::archive::text_iarchive ia_;
+ typedef std::ifstream ifs_;
+
+ const str_ input_path
+ = "/Users/erwann/projets/2009/Xcode/survival/build/Release/";
+ const str_ prior_path = input_path + "/prior";
+ const str_ x_vals_path = input_path + "/x_vals";
+ const str_ pars_path = input_path + "/pars";
+ const str_ batches_path = input_path + "/batches";
+
+ // Covariates
+ typedef val_ x_;
+ typedef vals_ x_vals_;
+ typedef range_cycle<> range_cycle_;
+ typedef range_cycle_::apply<x_vals_>::type x_cycle_;
+ x_vals_ x_vals;
+ {
+ ifs_ ifs(x_vals_path.c_str());
+ ia_ ia(ifs);
+ ia >> x_vals;
+ }
+ x_cycle_ x_cycle;
+
+ // Prior
+ typedef math::normal_distribution<val_> mprior_;
+ mprior_ mprior;
+ {
+ typedef math::meta_distribution_primitives<mprior_>::type prim_;
+ ifs_ ifs(prior_path.c_str());
+ ia_ ia(ifs);
+ prim_ prim;
+ ia >> prim;
+ mprior = prim;
+ }
+
+ // Records
+ // TODO n_bath should be deduced from input (getline)
+ const long n_batch = 1e4;
+
+ typedef surv::data::record<val_> record_;
+ typedef std::vector<record_> records_;
+ typedef range_size<records_>::type size_type;
+ typedef range_iterator<records_>::type it_record_;
+
+ records_ records;
+ ifs_ ifs_batches(batches_path.c_str());
+ if(!ifs_batches.good()){
+ str_ str = "error opening : ";
+ str.append( batches_path );
+ throw std::runtime_error(str);
+ }
+ ia_ ia_batches(ifs_batches);
+
+ // True pars
+ typedef val_ par_;
+ typedef vals_ pars_;
+ const long n_true_pars_kss = -1;//n_batch;
+ pars_ true_pars;
+ true_pars.reserve(n_batch);
+ {
+ ifs_ ifs_pars(pars_path.c_str());
+ if(!ifs_pars.good()){
+ str_ str = "error opening : ";
+ str.append( pars_path );
+ throw std::runtime_error(str);
+ }
+ ia_ ia_pars(ifs_pars);
+
+ for(long i = 0; i< n_batch; i++){
+ val_ par;
+ ia_pars >> par;
+ true_pars.push_back(par);
+ }
+ }
+ if(n_true_pars_kss>0){ // Check that F_n(true_pars) agrees with mprior
+ BOOST_ASSERT( n_batch % n_true_pars_kss == 0);
+ vals_ true_pars_kss;
+ true_pars_kss.reserve(n_true_pars_kss);
+ statistics::empirical_cdf::sequential_kolmogorov_smirnov_distance(
+ mprior,
+ boost::begin(true_pars),
+ boost::end(true_pars),
+ n_true_pars_kss,
+ std::back_inserter(true_pars_kss)
+ );
+ out << "true pars kss : ";
+ std::copy(
+ boost::begin( true_pars_kss ),
+ boost::end( true_pars_kss ),
+ std::ostream_iterator<val_>(out," ")
+ );
+ out << std::endl;
+ }
+
+ // [ Output ]
+ typedef boost::archive::text_oarchive oa_;
+ typedef std::ofstream ofs_;
+
+ const char* t_pars_path = "./t_pars";
+ const char* cg_path = "./cg";
+ const char* ks_path = "./kss";
+
+ // Cg
+ ofs_ ofs_cg( cg_path );
+ vals_ cgs; cgs.reserve(n_batch);
+
+ // Targets
+ const long n_t_pars = 1e3;
+ pars_ t_pars;
+ t_pars.reserve( n_t_pars );
+ ofs_ ofs_t_pars(t_pars_path);
+ oa_ oa_t_pars(ofs_t_pars);
+
+ // [ Local ]
+
+ // Monitoring
+ const long n_batch_mod = 1e3;
+
+ // Events
+ const val_ entry_bound = const_::inf_;
+ typedef surv::data::event<val_> event_;
+ typedef std::vector<event_> events_;
+ typedef range_iterator<events_>::type it_event_;
+
+ // Urng
+ typedef boost::mt19937 urng_;
+ urng_ urng;
+
+ // Unif
+ typedef math::uniform_distribution<val_> munif_;
+
+ // Proposal
+ typedef mprior_ mproposal_;
+ pars_ p_pars;
+ const long n_proposal = 1e3; //1e4 recommended but takes longer
+ p_pars.reserve(n_proposal);
+
+ // This one size fits all proposal is likely to result in a small effective
+ // sample size. Check out<<pws below to determine n_proposal.
+ mproposal_ mproposal = mprior;
+ vals_ p_lpdfs;
+ p_lpdfs.reserve(n_proposal);
+
+ // Weights
+ const val_ max_log = 100.0;
+ typedef iw::prepare_weights<val_> pws_;
+ pws_ pws( max_log );
+ vals_ iws;
+ iws.reserve( n_proposal );
+
+ // Model
+ typedef surv::model::exponential::model<val_> model_;
+ model_ model;
+ typedef model::prior_model_dataset_<mprior_,model_,x_cycle_,events_> pmd_;
+
+ { // [ Proposal sample ]
+ const long n_p_pars_kss = -1;//n_batch;
+ boost::timer t;
+ p_pars.clear();
+ generate_n(
+ std::back_inserter(p_pars),
+ n_proposal,
+ mproposal,
+ urng
+ );
+
+ if(n_p_pars_kss>0){ // Check that F_n(p_pars) agrees with mproposal
+ BOOST_ASSERT( n_batch % n_p_pars_kss == 0);
+ vals_ p_pars_kss;
+ p_pars_kss.reserve(n_p_pars_kss);
+ statistics::empirical_cdf::sequential_kolmogorov_smirnov_distance(
+ mproposal,
+ boost::begin(p_pars),
+ boost::end(p_pars),
+ n_p_pars_kss,
+ std::back_inserter(p_pars_kss)
+ );
+ out << "p_pars kss : ";
+ std::copy(
+ boost::begin( p_pars_kss ),
+ boost::end( p_pars_kss ),
+ std::ostream_iterator<val_>(out," ")
+ );
+ out << std::endl;
+ }
+
+ p_lpdfs.clear();
+ math::transform<math::fun_wrap::log_unnormalized_pdf_>(
+ mproposal,
+ boost::begin( p_pars ),
+ boost::end( p_pars ),
+ std::back_inserter( p_lpdfs )
+ );
+ BOOST_ASSERT( size(p_pars) == size(p_lpdfs) );
+ }
+
+ // Loops over batches
+ for(unsigned i = 0; i<n_batch; i++){
+ records.clear();
+ ia_batches >> records;
+ events_ events;
+ events.reserve( size(records) );
+ surv::data::events(
+ begin(records),
+ end(records),
+ entry_bound,
+ std::back_inserter(events)
+ );
+ size_type n_event = size( events );
+
+ x_cycle = range_cycle_::make( x_vals, 0, n_event );
+ BOOST_ASSERT( size( x_cycle ) >= n_event );
+ x_cycle.advance_end( - ( size( x_cycle ) - n_event ) );
+ BOOST_ASSERT( size( x_cycle ) == n_event );
+
+ { // [ Posterior sample ]
+ boost::timer t;
+ pmd_ pmd( mprior, model, x_cycle, events);
+ iws.resize( size( p_pars ) );
+
+ // [Warning: Bug]
+ // The Cook-Gelman statistics come out higher at the ends of the
+ // [0,1] range than in the middle
+ //
+ // model::log_posteriors<val_>(
+ // pmd,
+ // boost::begin( p_pars ),
+ // boost::end( p_pars ),
+ // boost::begin( p_lpdfs ),
+ // boost::begin( iws )
+ // );
+ // out << std::endl << "iws:";
+ // copy(
+ // boost::begin( iws ),
+ // boost::next( boost::begin( iws ), 10),
+ // std::ostream_iterator<val_>(std::cout, " ")
+ // );
+ // For testing purposes only:
+ // We expect iws2 == iws. But small differences (precision error?)
+ // vals_ iws2 = iws;
+ // model::log_likelihoods<val_>(
+ // pmd,
+ // boost::begin( p_pars ),
+ // boost::end( p_pars ),
+ // boost::begin( iws2 )
+ // );
+ // out << std::endl << "iws2:";
+ // copy(
+ // boost::begin( iws2 ),
+ // boost::next( boost::begin( iws2 ), 10),
+ // std::ostream_iterator<val_>(std::cout, " ")
+ // );
+
+ // Temporary fix, until bug above is resolved.
+ BOOST_ASSERT(
+ make_distribution_primitives(mprior) ==
+ make_distribution_primitives(mproposal)
+ );
+ model::log_likelihoods<val_>(
+ pmd,
+ boost::begin( p_pars ),
+ boost::end( p_pars ),
+ boost::begin( iws )
+ );
+
+ pws(
+ boost::begin( iws ),
+ boost::end( iws ),
+ boost::begin( p_pars )
+ );
+ t_pars.clear();
+ t_pars.resize( n_t_pars );
+ is::generate(
+ urng,
+ boost::begin( iws ),
+ boost::end( iws ),
+ boost::begin( p_pars ),
+ boost::begin( t_pars ),
+ n_t_pars
+ );
+ oa_t_pars << t_pars;
+ val_ plt = statistics::empirical_cdf::proportion_less_than(
+ boost::begin( t_pars ),
+ boost::end( t_pars ),
+ true_pars[i]
+ );
+ cgs.push_back( plt );
+ ofs_cg << plt << ' ';
+ ofs_cg.flush();
+ const char* str = "i = %1%, t = %2%, pws = %3%";
+ if(i%n_batch_mod == 0){
+ format f(str);
+ f % i % t.elapsed() % pws;
+ out << std::endl << f.str();
+ }
+ }// Proposal sample
+ }// loop over batches
+
+ { // Kolmorov Smirnov of the cgs
+ long n_ks = n_batch;
+ BOOST_ASSERT( n_batch % n_ks == 0 );
+ vals_ kss;
+ kss.reserve( n_batch );
+ statistics::empirical_cdf::sequential_kolmogorov_smirnov_distance(
+ munif_(),
+ boost::begin( cgs ),
+ boost::end( cgs ),
+ n_ks,
+ std::back_inserter(kss)
+ );
+ ofs_ ofs_ks(ks_path);
+ std::copy(
+ boost::begin( kss ),
+ boost::end( kss ),
+ std::ostream_iterator<val_>(ofs_ks," ")
+ );
+ }
+
+ out << "<-" << std::endl;
+}
+
Added: sandbox/statistics/survival_model/libs/statistics/survival/model/example/posterior_analysis.h
==============================================================================
--- (empty file)
+++ sandbox/statistics/survival_model/libs/statistics/survival/model/example/posterior_analysis.h 2009-08-31 21:58:21 EDT (Mon, 31 Aug 2009)
@@ -0,0 +1,14 @@
+///////////////////////////////////////////////////////////////////////////////
+// statistics::survival::model::example::model::posterior_analysis.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_SURVIVAL_EXAMPLE_MODEL_POSTERIOR_ANALYSIS_HPP_ER_2009
+#define LIBS_SURVIVAL_EXAMPLE_MODEL_POSTERIOR_ANALYSIS_HPP_ER_2009
+#include <ostream>
+
+void example_posterior_analysis(std::ostream&);
+
+#endif
\ No newline at end of file
Added: sandbox/statistics/survival_model/libs/statistics/survival/model/src/main.cpp
==============================================================================
--- (empty file)
+++ sandbox/statistics/survival_model/libs/statistics/survival/model/src/main.cpp 2009-08-31 21:58:21 EDT (Mon, 31 Aug 2009)
@@ -0,0 +1,11 @@
+#include <iostream>
+#include <libs/statistics/survival/model/example/exponential.h>
+#include <libs/statistics/survival/model/example/posterior_analysis.h>
+
+int main(){
+
+ example_exponential(std::cout);
+ //example_posterior_analysis(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