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