work/mixpp/bdm/estim/libPF.h

Go to the documentation of this file.

#ifndef PF_H
#define PF_H

#include <itpp/itbase.h>
#include "../stat/libEF.h"
#include "../math/libDC.h"

using namespace itpp;

class PF : public BM {
protected:
        int n;
        eEmp ePdf;
        vec &w;
        Array<vec> &samples;
        mpdf &par;
        mpdf &obs;
public:
        PF ( const RV &rv0, mpdf &par0,  mpdf &obs0, int n ) :BM(rv0),
                        n ( 1 ),ePdf ( rv0 ),w ( ePdf._w() ),samples ( ePdf._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 );
};

template<class BM_T>
class MPF : public PF {
        BM_T* Bms[1000];
        RV rvc;
        emix est;
public:
        MPF ( const RV &rv0, mpdf &par0, mpdf &obs0, int n, const BM_T &BMcond0 ) :
                        PF ( rv0,par0,obs0,n ), est(rv0,w) {
                //
                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();
                        est.set_parameters(i,&pom);
                }
        };
        ~MPF() {
        }
        void bayes (const vec &dt );
        epdf& _epdf(){return est;}
};

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(), sum=0.0;
        
        for ( i=0;i<n;i++ ) {
        //generate new samples from paramater evolution model;
                samples(i) = par.samplecond( samples( i ), lls(i) );
                if (lls(i)>mlls) mlls=lls(i); //find maximum
        }
        // 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( i ) /=sum; //?
        
        ind = ePdf.resample();
        
}

#endif // KF_H

---