- Main Page
- Related Pages
- Modules
- Classes
- Files

# work/git/mixpp/bdm/estim/libPF.h

Go to the documentation of this file.

    00001
    00013 #ifndef PF_H
    00014 #define PF_H
    00015
    00016 #include <itpp/itbase.h>
    00017 #include "../stat/libEF.h"
    00018 #include "../math/libDC.h"
    00019
    00020 namespace bdm{
    00021
    00028 class PF : public BM {
    00029 protected:
    00031         int n;
    00033         eEmp est;
    00035         vec &_w;
    00037         Array<vec> &_samples;
    00039         mpdf &par;
    00041         mpdf &obs;
    00042 public:
    00044         PF ( const RV &rv0, mpdf &par0,  mpdf &obs0, int n0 ) :BM ( rv0 ),
    00045                         n ( n0 ),est ( rv0,n ),_w ( est._w() ),_samples ( est._samples() ),
    00046                         par ( par0 ), obs ( obs0 ) {};
    00047
    00049         void set_est ( const epdf &epdf0 );
    00050         void bayes ( const vec &dt );
    00052         vec* __w(){return &_w;}
    00053 };
    00054
    00061 template<class BM_T>
    00062
    00063 class MPF : public PF {
    00064         BM_T* Bms[10000];
    00065
    00067
    00068 class mpfepdf : public epdf  {
    00069         protected:
    00070                 eEmp &E;
    00071                 vec &_w;
    00072                 Array<const epdf*> Coms;
    00073         public:
    00074                 mpfepdf ( eEmp &E0, const RV &rvc ) :
    00075                                 epdf ( RV( ) ), E ( E0 ),  _w ( E._w() ),
    00076                                 Coms ( _w.length() ) {
    00077                         rv.add ( E._rv() );
    00078                         rv.add ( rvc );
    00079                 };
    00080
    00081                 void set_elements ( int &i, double wi, const epdf* ep )
    00082                 {_w ( i ) =wi; Coms ( i ) =ep;};
    00083
    00084                 vec mean() const {
    00085                         // ugly
    00086                         vec pom=zeros ( ( Coms ( 0 )->_rv() ).count() );
    00087                         for ( int i=0; i<_w.length(); i++ ) {pom += Coms ( i )->mean() * _w ( i );}
    00088                         return concat ( E.mean(),pom );
    00089                 }
    00090                 vec variance() const {
    00091                         // ugly
    00092                         vec pom=zeros ( ( Coms ( 0 )->_rv() ).count() );
    00093                         vec pom2=zeros ( ( Coms ( 0 )->_rv() ).count() );
    00094                         for ( int i=0; i<_w.length(); i++ ) {
    00095                                 pom += Coms ( i )->mean() * _w ( i );
    00096                                 pom2 += (Coms ( i )->variance() + pow(Coms(i)->mean(),2)) * _w ( i );}
    00097                         return concat ( E.variance(),pom2-pow(pom,2) );
    00098                 }
    00099
    00100                 vec sample() const {it_error ( "Not implemented" );return 0;}
    00101
    00102                 double evallog ( const vec &val ) const {it_error ( "not implemented" ); return 0.0;}
    00103         };
    00104
    00106         mpfepdf jest;
    00107
    00108 public:
    00110         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 ) {
    00111                 //
    00112                 //TODO test if rv and BMcond.rv are compatible.
    00113                 rv.add ( rvlin );
    00114                 //
    00115
    00116                 if ( n>10000 ) {it_error ( "increase 10000 here!" );}
    00117
    00118                 for ( int i=0;i<n;i++ ) {
    00119                         Bms[i] = new BM_T ( BMcond0 ); //copy constructor
    00120                         const epdf& pom=Bms[i]->_epdf();
    00121                         jest.set_elements ( i,1.0/n,&pom );
    00122                 }
    00123         };
    00124
    00125         ~MPF() {
    00126         }
    00127
    00128         void bayes ( const vec &dt );
    00129         const epdf& _epdf() const {return jest;}
    00130         const epdf* _e() const {return &jest;} //Fixme: is it useful?
    00132         void set_est ( const epdf& epdf0 ) {
    00133                 PF::set_est ( epdf0 );  // sample params in condition
    00134                 // copy conditions to BMs
    00135
    00136                 for ( int i=0;i<n;i++ ) {Bms[i]->condition ( _samples ( i ) );}
    00137         }
    00138
    00140         BM* _BM(int i){return Bms[i];}
    00141 };
    00142
    00143 template<class BM_T>
    00144 void MPF<BM_T>::bayes ( const vec &dt ) {
    00145         int i;
    00146         vec lls ( n );
    00147         vec llsP ( n );
    00148         ivec ind;
    00149         double mlls=-std::numeric_limits<double>::infinity();
    00150
    00151         #pragma omp parallel for
    00152         for ( i=0;i<n;i++ ) {
    00153                 //generate new samples from paramater evolution model;
    00154                 _samples ( i ) = par.samplecond ( _samples ( i ), llsP ( i ) );
    00155                 Bms[i]->condition ( _samples ( i ) );
    00156                 Bms[i]->bayes ( dt );
    00157                 lls ( i ) = Bms[i]->_ll(); // lls above is also in proposal her must be lls(i) =, not +=!!
    00158                 if ( lls ( i ) >mlls ) mlls=lls ( i ); //find maximum likelihood (for numerical stability)
    00159         }
    00160
    00161         double sum_w=0.0;
    00162         // compute weights
    00163         #pragma omp parallel for
    00164         for ( i=0;i<n;i++ ) {
    00165                 _w ( i ) *= exp ( lls ( i ) - mlls ); // multiply w by likelihood
    00166                 sum_w+=_w(i);
    00167         }
    00168
    00169         if ( sum_w  >0.0 ) {
    00170                 _w /=sum_w; //?
    00171         } else {
    00172                 cout<<"sum(w)==0"<<endl;
    00173         }
    00174
    00175
    00176         double eff = 1.0/ ( _w*_w );
    00177         if ( eff < ( 0.3*n ) ) {
    00178                 ind = est.resample();
    00179                 // Resample Bms!
    00180
    00181                 #pragma omp parallel for
    00182                 for ( i=0;i<n;i++ ) {
    00183                         if ( ind ( i ) !=i ) {//replace the current Bm by a new one
    00184                                 //fixme this would require new assignment operator
    00185                                 // *Bms[i] = *Bms[ind ( i ) ];
    00186
    00187                                 // poor-man's solution: replicate constructor here
    00188                                 // copied from MPF::MPF
    00189                                 delete Bms[i];
    00190                                 Bms[i] = new BM_T ( *Bms[ind ( i ) ] ); //copy constructor
    00191                                 const epdf& pom=Bms[i]->_epdf();
    00192                                 jest.set_elements ( i,1.0/n,&pom );
    00193                         }
    00194                 };
    00195                 cout << '.';
    00196         }
    00197 }
    00198
    00199 }
    00200 #endif // KF_H
    00201
    00202