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

- File List
- File Members

# libPF.h

Go to the documentation of this file.

    00001
    00013 #ifndef PF_H
    00014 #define PF_H
    00015
    00016
    00017 #include "../stat/libEF.h"
    00018
    00019 namespace bdm {
    00020
    00027 class PF : public BM {
    00028 protected:
    00030         int n;
    00032         eEmp est;
    00034         vec &_w;
    00036         Array<vec> &_samples;
    00038         mpdf *par;
    00040         mpdf *obs;
    00041
    00043         RESAMPLING_METHOD resmethod;
    00044
    00047
    00049         bool opt_L_smp;
    00051         bool opt_L_wei;
    00053
    00054 public:
    00057         PF ( ) :est(), _w ( est._w() ),_samples ( est._samples() ), opt_L_smp ( false ), opt_L_wei ( false ) {LIDs.set_size ( 5 );};
    00058         /*      PF ( mpdf *par0, mpdf *obs0, epdf *epdf0, int n0 ) :
    00059                                 est ( ),_w ( est._w() ),_samples ( est._samples() ),opt_L_smp(false), opt_L_wei(false)
    00060                 { set_parameters ( par0,obs0,n0 ); set_statistics ( ones ( n0 ),epdf0 ); };*/
    00061         void set_parameters ( mpdf *par0, mpdf *obs0, int n0, RESAMPLING_METHOD rm=SYSTEMATIC )
    00062         { par = par0; obs=obs0; n=n0; resmethod= rm;};
    00063         void set_statistics ( const vec w0, epdf *epdf0 ) {est.set_statistics ( w0,epdf0 );};
    00066 //      void set_est ( const epdf &epdf0 );
    00067         void set_options ( const string &opt ) {
    00068                 BM::set_options(opt);
    00069                 opt_L_wei= ( opt.find ( "logweights" ) !=string::npos );
    00070                 opt_L_smp= ( opt.find ( "logsamples" ) !=string::npos );
    00071         }
    00072         void bayes ( const vec &dt );
    00074         vec* __w() {return &_w;}
    00075 };
    00076
    00083 template<class BM_T>
    00084 class MPF : public PF {
    00085         Array<BM_T*> BMs;
    00086
    00088
    00089 class mpfepdf : public epdf  {
    00090         protected:
    00091                 eEmp &E;
    00092                 vec &_w;
    00093                 Array<const epdf*> Coms;
    00094         public:
    00095                 mpfepdf ( eEmp &E0 ) :
    00096                                 epdf ( ), E ( E0 ),  _w ( E._w() ),
    00097                                 Coms ( _w.length() ) {
    00098                 };
    00100                 void read_statistics ( Array<BM_T*> &A ) {
    00101                         dim = E.dimension() +A ( 0 )->posterior().dimension();
    00102                         for ( int i=0; i<_w.length() ;i++ ) {Coms ( i ) = A ( i )->_e();}
    00103                 }
    00105                 void set_elements ( int &i, double wi, const epdf* ep )
    00106                 {_w ( i ) =wi; Coms ( i ) =ep;};
    00107
    00108                 void set_parameters ( int n ) {
    00109                         E.set_parameters ( n, false );
    00110                         Coms.set_length ( n );
    00111                 }
    00112                 vec mean() const {
    00113                         // ugly
    00114                         vec pom=zeros ( Coms ( 0 )->dimension() );
    00115                         for ( int i=0; i<_w.length(); i++ ) {pom += Coms ( i )->mean() * _w ( i );}
    00116                         return concat ( E.mean(),pom );
    00117                 }
    00118                 vec variance() const {
    00119                         // ugly
    00120                         vec pom=zeros ( Coms ( 0 )->dimension() );
    00121                         vec pom2=zeros ( Coms ( 0 )->dimension() );
    00122                         for ( int i=0; i<_w.length(); i++ ) {
    00123                                 pom += Coms ( i )->mean() * _w ( i );
    00124                                 pom2 += ( Coms ( i )->variance() + pow ( Coms ( i )->mean(),2 ) ) * _w ( i );
    00125                         }
    00126                         return concat ( E.variance(),pom2-pow ( pom,2 ) );
    00127                 }
    00128                 void qbounds ( vec &lb, vec &ub, double perc=0.95 ) const {
    00129                         //bounds on particles
    00130                         vec lbp;
    00131                         vec ubp;
    00132                         E.qbounds ( lbp,ubp );
    00133
    00134                         //bounds on Components
    00135                         int dimC=Coms ( 0 )->dimension();
    00136                         int j;
    00137                         // temporary
    00138                         vec lbc(dimC);
    00139                         vec ubc(dimC);
    00140                         // minima and maxima
    00141                         vec Lbc(dimC);
    00142                         vec Ubc(dimC);
    00143                         Lbc = std::numeric_limits<double>::infinity();
    00144                         Ubc = -std::numeric_limits<double>::infinity();
    00145
    00146                         for ( int i=0;i<_w.length();i++ ) {
    00147                                 // check Coms
    00148                                 Coms ( i )->qbounds ( lbc,ubc );
    00149                                 for ( j=0;j<dimC; j++ ) {
    00150                                         if ( lbc ( j ) <Lbc ( j ) ) {Lbc ( j ) =lbc ( j );}
    00151                                         if ( ubc ( j ) >Ubc ( j ) ) {Ubc ( j ) =ubc ( j );}
    00152                                 }
    00153                         }
    00154                         lb=concat(lbp,Lbc);
    00155                         ub=concat(ubp,Ubc);
    00156                 }
    00157
    00158                 vec sample() const {it_error ( "Not implemented" );return 0;}
    00159
    00160                 double evallog ( const vec &val ) const {it_error ( "not implemented" ); return 0.0;}
    00161         };
    00162
    00164         mpfepdf jest;
    00165
    00167         bool opt_L_mea;
    00168
    00169 public:
    00171         MPF () : PF (), jest ( est ) {};
    00172         void set_parameters ( mpdf *par0, mpdf *obs0, int n0, RESAMPLING_METHOD rm=SYSTEMATIC ) {
    00173                 PF::set_parameters ( par0, obs0, n0, rm );
    00174                 jest.set_parameters ( n0 );//duplication of rm
    00175                 BMs.set_length ( n0 );
    00176         }
    00177         void set_statistics ( epdf *epdf0, const BM_T* BMcond0 ) {
    00178
    00179                 PF::set_statistics ( ones ( n ) /n, epdf0 );
    00180                 // copy
    00181                 for ( int i=0;i<n;i++ ) { BMs ( i ) = new BM_T ( *BMcond0 ); BMs ( i )->condition ( _samples ( i ) );}
    00182
    00183                 jest.read_statistics ( BMs );
    00184                 //options
    00185         };
    00186
    00187         void bayes ( const vec &dt );
    00188         const epdf& posterior() const {return jest;}
    00189         const epdf* _e() const {return &jest;} //Fixme: is it useful?
    00191         /*      void set_est ( const epdf& epdf0 ) {
    00192                         PF::set_est ( epdf0 );  // sample params in condition
    00193                         // copy conditions to BMs
    00194
    00195                         for ( int i=0;i<n;i++ ) {BMs(i)->condition ( _samples ( i ) );}
    00196                 }*/
    00197         void set_options ( const string &opt ) {
    00198                 PF::set_options ( opt );
    00199                 opt_L_mea = ( opt.find ( "logmeans" ) !=string::npos );
    00200         }
    00201
    00203         BM* _BM ( int i ) {return BMs ( i );}
    00204 };
    00205
    00206 template<class BM_T>
    00207 void MPF<BM_T>::bayes ( const vec &dt ) {
    00208         int i;
    00209         vec lls ( n );
    00210         vec llsP ( n );
    00211         ivec ind;
    00212         double mlls=-std::numeric_limits<double>::infinity();
    00213
    00214 #pragma omp parallel for
    00215         for ( i=0;i<n;i++ ) {
    00216                 //generate new samples from paramater evolution model;
    00217                 _samples ( i ) = par->samplecond ( _samples ( i ) );
    00218                 llsP ( i ) = par->_e()->evallog ( _samples ( i ) );
    00219                 BMs ( i )->condition ( _samples ( i ) );
    00220                 BMs ( i )->bayes ( dt );
    00221                 lls ( i ) = BMs ( i )->_ll(); // lls above is also in proposal her must be lls(i) =, not +=!!
    00222                 if ( lls ( i ) >mlls ) mlls=lls ( i ); //find maximum likelihood (for numerical stability)
    00223         }
    00224
    00225         double sum_w=0.0;
    00226         // compute weights
    00227 #pragma omp parallel for
    00228         for ( i=0;i<n;i++ ) {
    00229                 _w ( i ) *= exp ( lls ( i ) - mlls ); // multiply w by likelihood
    00230                 sum_w+=_w ( i );
    00231         }
    00232
    00233         if ( sum_w  >0.0 ) {
    00234                 _w /=sum_w; //?
    00235         }
    00236         else {
    00237                 cout<<"sum(w)==0"<<endl;
    00238         }
    00239
    00240
    00241         double eff = 1.0/ ( _w*_w );
    00242         if ( eff < ( 0.3*n ) ) {
    00243                 ind = est.resample ( resmethod );
    00244                 // Resample Bms!
    00245
    00246 #pragma omp parallel for
    00247                 for ( i=0;i<n;i++ ) {
    00248                         if ( ind ( i ) !=i ) {//replace the current Bm by a new one
    00249                                 //fixme this would require new assignment operator
    00250                                 // *Bms[i] = *Bms[ind ( i ) ];
    00251
    00252                                 // poor-man's solution: replicate constructor here
    00253                                 // copied from MPF::MPF
    00254                                 delete BMs ( i );
    00255                                 BMs ( i ) = new BM_T ( *BMs ( ind ( i ) ) ); //copy constructor
    00256                                 const epdf& pom=BMs ( i )->posterior();
    00257                                 jest.set_elements ( i,1.0/n,&pom );
    00258                         }
    00259                 };
    00260                 cout << '.';
    00261         }
    00262 }
    00263
    00264 }
    00265 #endif // KF_H
    00266
    00267