Changeset 32 for bdm/estim/libPF.h

Timestamp:

03/03/08 13:00:32 (16 years ago)

Author:

smidl

Message:

test KF : estimation of R in KF is not possible! Likelihood of y_t is growing when R -> 0

Files:

• bdm/estim/libPF.h (modified) (1 diff)

bdm/estim/libPF.h

@@ -15 +15 @@
 
 #include <itpp/itbase.h>
-#include "../stat/libBM.h"
+#include "../stat/libEF.h"
 #include "../math/libDC.h"
 
 using namespace itpp;
 
-enum RESAMPLING_METHOD { MULTINOMIAL = 0, STRATIFIED = 1, SYSTEMATIC = 3 };
+/*!
+* \brief Trivial particle filter with proposal density equal to parameter evolution model.
 
-/*!
-* \brief A Particle Filter prototype
+Posterior density is represented by a weighted empirical density (\c eEmp ).
+*/
 
-Bayesian Filtering equations hold.
-*/
-class PF : public BM { 
+class PF : public BM {
 protected:
-        int n; //number of particles
-        vec w; //particle weights
-        Uniform_RNG URNG; //used for resampling
-        
+        //!number of particles;
+        int n;
+        //!posterior density
+        eEmp est;
+        //! pointer into \c eEmp
+        vec &_w;
+        //! pointer into \c eEmp
+        Array<vec> &_samples;
+        //! Parameter evolution model
+        mpdf &par;
+        //! Observation model
+        mpdf &obs;
 public:
-        //! Returns indexes of particles that should be resampled. The ordering MUST guarantee inplace replacement. (Important for MPF.)
-        ivec resample(RESAMPLING_METHOD method = SYSTEMATIC);
-        PF (vec w);
+        PF ( const RV &rv0, mpdf &par0,  mpdf &obs0, int n0 ) :BM ( rv0 ),
+                        n ( n0 ),est ( rv0,n ),_w ( est._w() ),_samples ( est._samples() ),
+                        par ( par0 ), obs ( obs0 ) {};
+
+//      void set_parametres(mpdf &par0, mpdf &obs0) {par=&par0;obs=&obs0;};
+        void set_est ( const epdf &epdf0 );
         //TODO remove or implement bayes()!
-        void bayes(const vec &dt){};
-        epdf* _epdf(){return NULL;}
+        void bayes ( const vec &dt );
 };
 
 /*!
-* \brief Trivial particle filter with proposal density that is not conditioned on the data.
+\brief Marginalized Particle filter
 
-
+Trivial version: proposal $=$ parameter evolution, observation model is not used. (it is assumed to be part of BM).
 */
 
-class TrivialPF : public PF {
-        Array<vec> ptcls;
-        
-        bool is_proposal;
-        BM *prop;
-        mpdf *par;
-        mpdf *obs;
-        
+template<class BM_T>
+
+class MPF : public PF {
+        BM_T* Bms[1000];
+
+        //! internal class for MPDF providing composition of eEmp with external components
+
+class mpfepdf : public epdf  {
+        protected:
+                eEmp &E;
+                vec &_w;
+                Array<epdf*> Coms;
         public:
-        TrivialPF(mpdf &par, mpdf &obs, BM &prop, int n0);
-        TrivialPF(mpdf &par, mpdf &obs, int n0);
-        void bayes(const vec &dt, bool evalll);
+                mpfepdf ( eEmp &E0, const RV &rvc ) :
+                                epdf ( RV( ) ), E ( E0 ),  _w ( E._w() ),
+                                Coms ( _w.length() ) {
+                        rv.add ( E._rv() );
+                        rv.add ( rvc );
+                };
+
+                void set_elements ( int &i, double wi, epdf* ep )
+                {_w ( i ) =wi; Coms ( i ) =ep;};
+
+                vec mean() const {
+                        // ugly
+                        vec pom=zeros ( ( Coms ( 0 )->_rv() ).count() );
+
+                        for ( int i=0; i<_w.length(); i++ ) {pom += Coms ( i )->mean() * _w ( i );}
+
+                        return concat ( E.mean(),pom );
+                }
+
+                vec sample() const {it_error ( "Not implemented" );return 0;}
+
+                double evalpdflog ( const vec &val ) const {it_error ( "not implemented" ); return 0.0;}
+        };
+
+        //! estimate joining PF.est with conditional
+        mpfepdf jest;
+
+public:
+        //! Default constructor.
+        MPF ( const RV &rvlin, const RV &rvpf, mpdf &par0, mpdf &obs0, int n, const BM_T &BMcond0 ) : PF ( rvpf ,par0,obs0,n ),jest ( est,rvlin ) {
+                //
+                //TODO test if rv and BMcond.rv are compatible.
+                rv.add ( rvlin );
+                //
+
+                if ( n>1000 ) it_error ( "increase 1000 here!" );
+
+                for ( int i=0;i<n;i++ ) {
+                        Bms[i] = new BM_T ( BMcond0 ); //copy constructor
+                        epdf& pom=Bms[i]->_epdf();
+                        jest.set_elements ( i,1.0/n,&pom );
+                }
+        };
+
+        ~MPF() {
+        }
+
+        void bayes ( const vec &dt );
+        epdf& _epdf() {return jest;}
+
+        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 ) );}
+        }
 };
 
-class MPF : public TrivialPF {
-        Array<BM> Bms;
-        public:
-        MPF(BM &B, BM &prop, mpdf &obs, mpdf &par);
-        void bayes(vec &dt);    
-};
+template<class BM_T>
+void MPF<BM_T>::bayes ( const vec &dt ) {
+        int i;
+        vec lls ( n );
+        ivec ind;
+        double mlls=-std::numeric_limits<double>::infinity();
+        double sum=0.0;
+
+        for ( i=0;i<n;i++ ) {
+                //generate new samples from paramater evolution model;
+                _samples ( i ) = par.samplecond ( _samples ( i ), lls ( i ) );
+                Bms[i]->condition ( _samples ( i ) );
+                Bms[i]->bayes ( dt );
+                lls ( i ) += Bms[i]->_ll();
+
+                if ( lls ( i ) >mlls ) mlls=lls ( i ); //find maximum likelihood (for numerical stability)
+        }
+
+        // compute weights
+        for ( i=0;i<n;i++ ) {
+                _w ( i ) *= exp ( lls ( i ) - mlls ); // multiply w by likelihood
+        }
+
+        //renormalize
+        for ( i=0;i<n;i++ ) {sum+=_w ( i );};
+
+        _w /=sum; //?
+
+        if ( ( _w*_w ) < ( 0.5*n ) ) {
+                ind = est.resample();
+                // Resample Bms!
+
+                for ( i=0;i<n;i++ ) {
+                        if ( ind ( i ) !=i ) {//replace the current Bm by a new one
+                                //fixme this would require new assignment operator
+                                // *Bms[i] = *Bms[ind ( i ) ];
+                                
+                                // poor-man's solution: replicate constructor here
+                                // copied from MPF::MPF
+                                delete Bms[i];
+                                Bms[i] = new BM_T ( *Bms[ind ( i ) ] ); //copy constructor
+                                epdf& pom=Bms[i]->_epdf();
+                                jest.set_elements ( i,1.0/n,&pom );
+                        }
+                };
+        }
+}
 
 #endif // KF_H