[Boost-commit] svn:boost r55985 - in sandbox/statistics/survival_model/libs/statistics/survival/model: doc example src

Subject: [Boost-commit] svn:boost r55985 - in sandbox/statistics/survival_model/libs/statistics/survival/model: doc example src
From: erwann.rogard_at_[hidden]
Date: 2009-09-02 21:53:31

Author: e_r
Date: 2009-09-02 21:53:30 EDT (Wed, 02 Sep 2009)
New Revision: 55985
URL: http://svn.boost.org/trac/boost/changeset/55985

Log:
m
Text files modified: 
   sandbox/statistics/survival_model/libs/statistics/survival/model/doc/readme.txt                 |    25 ++++----                                
   sandbox/statistics/survival_model/libs/statistics/survival/model/example/exponential.cpp        |     4                                         
   sandbox/statistics/survival_model/libs/statistics/survival/model/example/exponential.h          |     2                                         
   sandbox/statistics/survival_model/libs/statistics/survival/model/example/posterior_analysis.cpp |   117 +++++++++++++++------------------------ 
   sandbox/statistics/survival_model/libs/statistics/survival/model/example/posterior_analysis.h   |     2                                         
   sandbox/statistics/survival_model/libs/statistics/survival/model/src/main.cpp                   |     9 ++                                      
   6 files changed, 70 insertions(+), 89 deletions(-)

Modified: sandbox/statistics/survival_model/libs/statistics/survival/model/doc/readme.txt
==============================================================================
--- sandbox/statistics/survival_model/libs/statistics/survival/model/doc/readme.txt	(original)
+++ sandbox/statistics/survival_model/libs/statistics/survival/model/doc/readme.txt	2009-09-02 21:53:30 EDT (Wed, 02 Sep 2009)
@@ -1,5 +1,5 @@
 //////////////////////////////////////////////////////////////////////////////
-// statistics::survival::model::doc::readme                                              //
+// statistics::survival::model::doc::readme                                 //
 //                                                                          //
 //  (C) Copyright 2009 Erwann Rogard                                        //
 //  Use, modification and distribution are subject to the                   //
@@ -54,7 +54,7 @@
 
 [ History ]
 
-July 2009 : Current version
+July 2009 : Creation
 
 [ models ]
 
@@ -100,17 +100,16 @@
 (-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)<-
+i = 0, t = 0.018712, pws = (222.593,1,0.0431399,0.861)
+i = 1000, t = 0.018008, pws = (195.971,1,0.0198646,0.982)
+i = 2000, t = 0.018061, pws = (190.532,1,0.0202048,0.996)
+i = 3000, t = 0.017771, pws = (197.043,1,0.0549215,0.866)
+i = 4000, t = 0.017818, pws = (195.51,1,0.0671974,0.73)
+i = 5000, t = 0.017887, pws = (199.301,1,0.0386544,0.976)
+i = 6000, t = 0.017974, pws = (205.461,1,0.0598528,0.866)
+i = 7000, t = 0.018313, pws = (182.456,1,0.0708064,0.866)
+i = 8000, t = 0.017948, pws = (181.937,1,0.0683872,0.861)
+i = 9000, t = 0.017947, pws = (188.914,1,0.00839656,0.248)<-
 
 

Modified: sandbox/statistics/survival_model/libs/statistics/survival/model/example/exponential.cpp
==============================================================================
--- sandbox/statistics/survival_model/libs/statistics/survival/model/example/exponential.cpp	(original)
+++ sandbox/statistics/survival_model/libs/statistics/survival/model/example/exponential.cpp	2009-09-02 21:53:30 EDT (Wed, 02 Sep 2009)
@@ -1,5 +1,5 @@
 ///////////////////////////////////////////////////////////////////////////////
-// statistics::survival::model::example::model::exponential.cpp                          //
+// statistics::survival::model::example::model::exponential.cpp              //
 //                                                                           //
 //  Copyright 2009 Erwann Rogard. Distributed under the Boost                //
 //  Software License, Version 1.0. (See accompanying file                    //
@@ -188,7 +188,7 @@
     // [ Posteriors ]
 
     vals_ lpos;
-    model::log_posteriors<val_>(
+    model::log_posteriors2<val_>(
         model::make_prior_model_dataset(mprior,model,r_x,events),
         boost::begin(pars),
         boost::end(pars),

Modified: sandbox/statistics/survival_model/libs/statistics/survival/model/example/exponential.h
==============================================================================
--- sandbox/statistics/survival_model/libs/statistics/survival/model/example/exponential.h	(original)
+++ sandbox/statistics/survival_model/libs/statistics/survival/model/example/exponential.h	2009-09-02 21:53:30 EDT (Wed, 02 Sep 2009)
@@ -1,5 +1,5 @@
 ///////////////////////////////////////////////////////////////////////////////
-// statistics::survival::model::example::model::exponential.h                            //
+// statistics::survival::model::example::model::exponential.h                //
 //                                                                           //
 //  Copyright 2009 Erwann Rogard. Distributed under the Boost                //
 //  Software License, Version 1.0. (See accompanying file                    //

Modified: sandbox/statistics/survival_model/libs/statistics/survival/model/example/posterior_analysis.cpp
==============================================================================
--- sandbox/statistics/survival_model/libs/statistics/survival/model/example/posterior_analysis.cpp	(original)
+++ sandbox/statistics/survival_model/libs/statistics/survival/model/example/posterior_analysis.cpp	2009-09-02 21:53:30 EDT (Wed, 02 Sep 2009)
@@ -1,12 +1,12 @@
 ///////////////////////////////////////////////////////////////////////////////
-// statistics::survival::model::example::model::posterior_analysis.cpp                   //
+// 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 <iomanip>
 #include <fstream>
 #include <stdexcept>
 #include <string> //needed?
@@ -28,7 +28,8 @@
 #include <boost/foreach.hpp> 
 
 #include <boost/standard_distribution/distributions/normal.hpp>
-#include <boost/math/distributions/uniform.hpp>
+#include <boost/standard_distribution/distributions/uniform.hpp>
+#include <boost/dist_random/include.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>
@@ -104,7 +105,7 @@
 
     // Records
     // TODO n_bath should be deduced from input (getline)
-    const long n_batch                              = 1e4; 
+    const long n_batch                              = 5e4; 
 
     typedef surv::data::record<val_>                record_;
     typedef std::vector<record_>                    records_;
@@ -174,7 +175,7 @@
     vals_ cgs; cgs.reserve(n_batch);
 
     // Targets
-    const long n_t_pars = 1e3;
+    const long n_t_pars = 1e4;
     pars_ t_pars; 
     t_pars.reserve( n_t_pars );
     ofs_ ofs_t_pars(t_pars_path);
@@ -200,15 +201,15 @@
     
     // 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.
+    const long n_proposal = 1e4; 
+    // This one size fits all proposal is likely to result in a small ESS. You
+    // can monitor the ESS by calling out << pws. Ultimately, uniformity of
+    // the Cook-Gelman will determine if either the computations are wrong
+    // or if n_proposal is too small. n_proposal = 1e4 works for N(0,5).
     mproposal_  mproposal = mprior; 
-    vals_ p_lpdfs; 
-    p_lpdfs.reserve(n_proposal); 
+    pars_ p_pars (n_proposal); 
+    vals_ p_lpdfs(n_proposal);
+
 
     // Weights
     const val_ max_log = 100.0;
@@ -225,16 +226,17 @@
     {   // [ Proposal sample ]
         const long n_p_pars_kss = -1;//n_batch; 
         boost::timer t;
-        p_pars.clear();
-        generate_n(
-            std::back_inserter(p_pars),
+
+        generate_function_n<math::fun_wrap::log_unnormalized_pdf_>(
+            boost::begin( p_pars ), 
+            boost::begin( p_lpdfs ), 
             n_proposal,
-            mproposal,
+            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);
+            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(
@@ -253,13 +255,6 @@
             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) );
     }
 
@@ -287,50 +282,27 @@
             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 )  
-            );            
+            // The reason I uncomment the if statement is that log_posteriors2
+            // is supposed to work for either case. 
+            //if(
+            //    make_distribution_primitives( mprior ) 
+            //    ==make_distribution_primitives( mproposal )
+            //){
+            //    model::log_likelihoods<val_>(
+            //        pmd,
+            //        boost::begin( p_pars ),
+            //        boost::end( p_pars ),
+            //        boost::begin( iws )  
+            //    );            
+            //}else{
+                statistics::model::log_posteriors<val_>(
+                    pmd,
+                    boost::begin( p_pars ),
+                    boost::end( p_pars ),
+                    boost::begin( p_lpdfs ),   
+                    boost::begin( iws )  
+                ); 
+            //}
 
             pws(
                 boost::begin( iws ),
@@ -359,8 +331,11 @@
             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();  
+                f 
+                % i  
+                % t.elapsed() 
+                % pws; 
+                out << std::endl << std::setprecision( 10 ) <<  f.str();  
             }
         }// Proposal sample
     }// loop over batches

Modified: sandbox/statistics/survival_model/libs/statistics/survival/model/example/posterior_analysis.h
==============================================================================
--- sandbox/statistics/survival_model/libs/statistics/survival/model/example/posterior_analysis.h	(original)
+++ sandbox/statistics/survival_model/libs/statistics/survival/model/example/posterior_analysis.h	2009-09-02 21:53:30 EDT (Wed, 02 Sep 2009)
@@ -1,5 +1,5 @@
 ///////////////////////////////////////////////////////////////////////////////
-// statistics::survival::model::example::model::posterior_analysis.h                     //
+// statistics::survival::model::example::model::posterior_analysis.h         //
 //                                                                           //
 //  Copyright 2009 Erwann Rogard. Distributed under the Boost                //
 //  Software License, Version 1.0. (See accompanying file                    //

Modified: sandbox/statistics/survival_model/libs/statistics/survival/model/src/main.cpp
==============================================================================
--- sandbox/statistics/survival_model/libs/statistics/survival/model/src/main.cpp	(original)
+++ sandbox/statistics/survival_model/libs/statistics/survival/model/src/main.cpp	2009-09-02 21:53:30 EDT (Wed, 02 Sep 2009)
@@ -1,3 +1,10 @@
+///////////////////////////////////////////////////////////////////////////////
+// statistics::survival::model::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/survival/model/example/exponential.h>
 #include <libs/statistics/survival/model/example/posterior_analysis.h>
@@ -5,7 +12,7 @@
 int main(){
 
     example_exponential(std::cout);
-    //example_posterior_analysis(std::cout);
+    example_posterior_analysis(std::cout);
     
     return 0;
 }
\ No newline at end of file