- Main Page
- Related Pages
- Classes
- Files

- File List
- File Members

# particles.h

Go to the documentation of this file.

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