Changeset 283 for bdm/estim/libPF.h

Timestamp:

02/24/09 14:14:01 (15 years ago)

Author:

smidl

Message:

get rid of BMcond + adaptation in doprava/

Files:

• bdm/estim/libPF.h (modified) (9 diffs)

bdm/estim/libPF.h

@@ -40 +40 @@
         mpdf *obs;
 
+        //! which resampling method will be used
+        RESAMPLING_METHOD resmethod;
+
         //! \name Options
         //!@{
 
-        //! Log resampling times
-        bool opt_L_res;
         //! Log all samples
         bool opt_L_smp;
+        //! Log all samples
+        bool opt_L_wei;
         //!@}
 
@@ -52 +55 @@
         //! \name Constructors
         //!@{
-        PF ( ) :est(), _w ( est._w() ),_samples ( est._samples() ) {};
-        PF ( mpdf *par0, mpdf *obs0, epdf *epdf0, int n0 ) :
-                        est ( ),_w ( est._w() ),_samples ( est._samples() )
-        { set_parameters ( par0,obs0,n0 ); set_statistics ( ones ( n0 ),epdf0 ); };
-        void set_parameters ( mpdf *par0, mpdf *obs0, int n0 )
-        { par = par0; obs=obs0; n=n0; est.set_n ( n );};
-        void set_statistics ( const vec w0, epdf *epdf0 ) {est.set_parameters ( w0,epdf0 );};
+        PF ( ) :est(), _w ( est._w() ),_samples ( est._samples() ), opt_L_smp ( false ), opt_L_wei ( false ) {LIDs.set_size ( 5 );};
+        /*      PF ( mpdf *par0, mpdf *obs0, epdf *epdf0, int n0 ) :
+                                est ( ),_w ( est._w() ),_samples ( est._samples() ),opt_L_smp(false), opt_L_wei(false)
+                { set_parameters ( par0,obs0,n0 ); set_statistics ( ones ( n0 ),epdf0 ); };*/
+        void set_parameters ( mpdf *par0, mpdf *obs0, int n0, RESAMPLING_METHOD rm=SYSTEMATIC )
+        { par = par0; obs=obs0; n=n0; resmethod= rm;};
+        void set_statistics ( const vec w0, epdf *epdf0 ) {est.set_statistics ( w0,epdf0 );};
         //!@}
         //! Set posterior density by sampling from epdf0
-        void set_est ( const epdf &epdf0 );
-        void set_options(const string &opt){
-                opt_L_res= ( s.find ( "logres" ) !=string::npos );
-                opt_L_smp= ( s.find ( "logsmp" ) !=string::npos );              
+//      void set_est ( const epdf &epdf0 );
+        void set_options ( const string &opt ) {
+                BM::set_options(opt);
+                opt_L_wei= ( opt.find ( "logweights" ) !=string::npos );
+                opt_L_smp= ( opt.find ( "logsamples" ) !=string::npos );
         }
         void bayes ( const vec &dt );
@@ -78 +82 @@
 
 template<class BM_T>
-
 class MPF : public PF {
-        BM_T* Bms[10000];
+        Array<BM_T*> BMs;
 
         //! internal class for MPDF providing composition of eEmp with external components
@@ -94 +97 @@
                                 Coms ( _w.length() ) {
                 };
+                //! read statistics from MPF
+                void read_statistics ( Array<BM_T*> &A ) {
+                        dim = E.dimension() +A ( 0 )->posterior().dimension();
+                        for ( int i=0; i<_w.length() ;i++ ) {Coms ( i ) = A ( i )->_e();}
+                }
+                //! needed in resampling
                 void set_elements ( int &i, double wi, const epdf* ep )
                 {_w ( i ) =wi; Coms ( i ) =ep;};
 
-                void set_n ( int n ) {E.set_n ( n ); Coms.set_length ( n );}
+                void set_parameters ( int n ) {
+                        E.set_parameters ( n, false );
+                        Coms.set_length ( n );
+                }
                 vec mean() const {
                         // ugly
@@ -114 +126 @@
                         return concat ( E.variance(),pom2-pow ( pom,2 ) );
                 }
+                void qbounds ( vec &lb, vec &ub, double perc=0.95 ) const {
+                        //bounds on particles
+                        vec lbp;
+                        vec ubp;
+                        E.qbounds ( lbp,ubp );
+
+                        //bounds on Components
+                        int dimC=Coms ( 0 )->dimension();
+                        int j;
+                        // temporary
+                        vec lbc(dimC);
+                        vec ubc(dimC);
+                        // minima and maxima
+                        vec Lbc(dimC);
+                        vec Ubc(dimC);
+                        Lbc = std::numeric_limits<double>::infinity();
+                        Ubc = -std::numeric_limits<double>::infinity();
+
+                        for ( int i=0;i<_w.length();i++ ) {
+                                // check Coms
+                                Coms ( i )->qbounds ( lbc,ubc );
+                                for ( j=0;j<dimC; j++ ) {
+                                        if ( lbc ( j ) <Lbc ( j ) ) {Lbc ( j ) =lbc ( j );}
+                                        if ( ubc ( j ) >Ubc ( j ) ) {Ubc ( j ) =ubc ( j );}
+                                }
+                        }
+                        lb=concat(lbp,Lbc);
+                        ub=concat(ubp,Ubc);
+                }
 
                 vec sample() const {it_error ( "Not implemented" );return 0;}
@@ -125 +166 @@
         //! Log means of BMs
         bool opt_L_mea;
-        
+
 public:
         //! Default constructor.
-        MPF ( mpdf *par0, mpdf *obs0, int n, const BM_T &BMcond0 ) : PF (), jest ( est ) {
-                PF::set_parameters ( par0,obs0,n );
-                jest.set_n ( n );
-                //
-                //TODO test if rv and BMcond.rv are compatible.
-//              rv.add ( rvlin );
-                //
-
-                if ( n>10000 ) {it_error ( "increase 10000 here!" );}
-
-                for ( int i=0;i<n;i++ ) {
-                        Bms[i] = new BM_T ( BMcond0 ); //copy constructor
-                        const epdf& pom=Bms[i]->posterior();
-                        jest.set_elements ( i,1.0/n,&pom );
-                }
+        MPF () : PF (), jest ( est ) {};
+        void set_parameters ( mpdf *par0, mpdf *obs0, int n0, RESAMPLING_METHOD rm=SYSTEMATIC ) {
+                PF::set_parameters ( par0, obs0, n0, rm );
+                jest.set_parameters ( n0 );//duplication of rm
+                BMs.set_length ( n0 );
+        }
+        void set_statistics ( epdf *epdf0, const BM_T* BMcond0 ) {
+
+                PF::set_statistics ( ones ( n ) /n, epdf0 );
+                // copy
+                for ( int i=0;i<n;i++ ) { BMs ( i ) = new BM_T ( *BMcond0 ); BMs ( i )->condition ( _samples ( i ) );}
+
+                jest.read_statistics ( BMs );
+                //options
         };
-
-        ~MPF() {
-        }
 
         void bayes ( const vec &dt );
         const epdf& posterior() const {return jest;}
         const epdf* _e() const {return &jest;} //Fixme: is it useful?
-        //! Set postrior of \c rvc to samples from epdf0. Statistics of Bms are not re-computed! Use only for initialization!
-        void set_est ( const epdf& epdf0 ) {
-                PF::set_est ( epdf0 );  // sample params in condition
-                // copy conditions to BMs
-
-                for ( int i=0;i<n;i++ ) {Bms[i]->condition ( _samples ( i ) );}
-        }
-        void set_options(const string &opt){
-                PF:set_options(opt);
-                opt_L_mea = ( s.find ( "logmeans" ) !=string::npos );
-        }
-        
+        //! Set postrior of \c rvc to samples from epdf0. Statistics of BMs are not re-computed! Use only for initialization!
+        /*      void set_est ( const epdf& epdf0 ) {
+                        PF::set_est ( epdf0 );  // sample params in condition
+                        // copy conditions to BMs
+
+                        for ( int i=0;i<n;i++ ) {BMs(i)->condition ( _samples ( i ) );}
+                }*/
+        void set_options ( const string &opt ) {
+                PF::set_options ( opt );
+                opt_L_mea = ( opt.find ( "logmeans" ) !=string::npos );
+        }
+
         //!Access function
-        BM* _BM ( int i ) {return Bms[i];}
+        BM* _BM ( int i ) {return BMs ( i );}
 };
 
@@ -180 +217 @@
                 _samples ( i ) = par->samplecond ( _samples ( i ) );
                 llsP ( i ) = par->_e()->evallog ( _samples ( i ) );
-                Bms[i]->condition ( _samples ( i ) );
-                Bms[i]->bayes ( dt );
-                lls ( i ) = Bms[i]->_ll(); // lls above is also in proposal her must be lls(i) =, not +=!!
+                BMs ( i )->condition ( _samples ( i ) );
+                BMs ( i )->bayes ( dt );
+                lls ( i ) = BMs ( i )->_ll(); // lls above is also in proposal her must be lls(i) =, not +=!!
                 if ( lls ( i ) >mlls ) mlls=lls ( i ); //find maximum likelihood (for numerical stability)
         }
@@ -204 +241 @@
         double eff = 1.0/ ( _w*_w );
         if ( eff < ( 0.3*n ) ) {
-                ind = est.resample();
+                ind = est.resample ( resmethod );
                 // Resample Bms!
 
@@ -215 +252 @@
                                 // poor-man's solution: replicate constructor here
                                 // copied from MPF::MPF
-                                delete Bms[i];
-                                Bms[i] = new BM_T ( *Bms[ind ( i ) ] ); //copy constructor
-                                const epdf& pom=Bms[i]->posterior();
+                                delete BMs ( i );
+                                BMs ( i ) = new BM_T ( *BMs ( ind ( i ) ) ); //copy constructor
+                                const epdf& pom=BMs ( i )->posterior();
                                 jest.set_elements ( i,1.0/n,&pom );
                         }