- 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         shared_ptr<mpdf> par;
    00040         shared_ptr<mpdf> obs;
    00042         vec lls;
    00043
    00045         RESAMPLING_METHOD resmethod;
    00048         double res_threshold;
    00049
    00052
    00054         bool opt_L_smp;
    00056         bool opt_L_wei;
    00058
    00059 public:
    00062         PF ( ) : est(), _w ( est._w() ), _samples ( est._samples() ), opt_L_smp ( false ), opt_L_wei ( false ) {
    00063                 LIDs.set_size ( 5 );
    00064         };
    00065
    00066         void set_parameters (int n0, double res_th0=0.5, RESAMPLING_METHOD rm = SYSTEMATIC ) {
    00067                 n = n0;
    00068                 res_threshold = res_th0;
    00069                 resmethod = rm;
    00070         };
    00071         void set_model ( shared_ptr<mpdf> par0, shared_ptr<mpdf> obs0) {
    00072                 par = par0;
    00073                 obs = obs0;
    00074                 // set values for posterior
    00075                 est.set_rv(par->_rv());
    00076         };
    00077         void set_statistics ( const vec w0, const epdf &epdf0 ) {
    00078                 est.set_statistics ( w0, epdf0 );
    00079         };
    00080         void set_statistics ( const eEmp &epdf0 ) {
    00081                 bdm_assert_debug(epdf0._rv().equal(par->_rv()),"Incompatibel input");
    00082                 est=epdf0;
    00083         };
    00089         void set_options ( const string &opt ) {
    00090                 BM::set_options ( opt );
    00091                 opt_L_wei = ( opt.find ( "logweights" ) != string::npos );
    00092                 opt_L_smp = ( opt.find ( "logsamples" ) != string::npos );
    00093         }
    00095         virtual void bayes_gensmp();
    00097         virtual void bayes_weights();
    00099         virtual bool do_resampling(){
    00100                 double eff = 1.0 / ( _w * _w );
    00101                 return eff < ( res_threshold*n );
    00102         }
    00103         void bayes ( const vec &dt );
    00105         vec& __w() { return _w; }
    00107         vec& _lls() { return lls; }
    00109         RESAMPLING_METHOD _resmethod() const { return resmethod; }
    00111         const eEmp& posterior() const {return est;}
    00112
    00125         void from_setting(const Setting &set){
    00126                 par = UI::build<mpdf>(set,"parameter_pdf",UI::compulsory);
    00127                 obs = UI::build<mpdf>(set,"observation_pdf",UI::compulsory);
    00128
    00129                 prior_from_set(set);
    00130                 resmethod_from_set(set);
    00131                 // set resampling method
    00132                 //set drv
    00133                 //find potential input - what remains in rvc when we subtract rv
    00134                 RV u = par->_rvc().remove_time().subt( par->_rv() );
    00135                 //find potential input - what remains in rvc when we subtract x_t
    00136                 RV obs_u = obs->_rvc().remove_time().subt( par->_rv() );
    00137
    00138                 u.add(obs_u); // join both u, and check if they do not overlap
    00139
    00140                 set_drv(concat(obs->_rv(),u) );
    00141         }
    00143         void resmethod_from_set(const Setting &set){
    00144                 string resmeth;
    00145                 if (UI::get(resmeth,set,"resmethod",UI::optional)){
    00146                         if (resmeth=="systematic") {
    00147                                 resmethod= SYSTEMATIC;
    00148                         } else  {
    00149                                 if (resmeth=="multinomial"){
    00150                                         resmethod=MULTINOMIAL;
    00151                                 } else {
    00152                                         if (resmeth=="stratified"){
    00153                                                 resmethod= STRATIFIED;
    00154                                         } else {
    00155                                                 bdm_error("Unknown resampling method");
    00156                                         }
    00157                                 }
    00158                         }
    00159                 } else {
    00160                         resmethod=SYSTEMATIC;
    00161                 };
    00162                 if(!UI::get(res_threshold, set, "res_threshold", UI::optional)){
    00163                         res_threshold=0.5;
    00164                 }
    00165         }
    00167         void prior_from_set(const Setting & set){
    00168                 shared_ptr<epdf> pri = UI::build<epdf>(set,"prior",UI::compulsory);
    00169
    00170                 eEmp *test_emp=dynamic_cast<eEmp*>(&(*pri));
    00171                 if (test_emp) { // given pdf is sampled
    00172                         est=*test_emp;
    00173                 } else {
    00174                         int n;
    00175                         if (!UI::get(n,set,"n",UI::optional)){n=10;}
    00176                         // sample from prior
    00177                         set_statistics(ones(n)/n, *pri);
    00178                 }
    00179                 //validate();
    00180         }
    00181
    00182         void validate(){
    00183                 n=_w.length();
    00184                 lls=zeros(n);
    00185                 if (par->_rv()._dsize()>0) {
    00186                         bdm_assert(par->_rv()._dsize()==est.dimension(),"Mismatch of RV and dimension of posterior" );
    00187                 }
    00188         }
    00190         void resample(ivec &ind){
    00191                 est.resample(ind,resmethod);
    00192         }
    00194         Array<vec>& __samples(){return _samples;}
    00195 };
    00196 UIREGISTER(PF);
    00197
    00205 class MPF : public BM  {
    00206         protected:
    00208         shared_ptr<PF> pf;
    00210         Array<BM*> BMs;
    00211
    00213
    00214         class mpfepdf : public epdf  {
    00216                 shared_ptr<PF> &pf;
    00218                 Array<BM*> &BMs;
    00219         public:
    00221                 mpfepdf (shared_ptr<PF> &pf0, Array<BM*> &BMs0): epdf(), pf(pf0), BMs(BMs0) { };
    00223                 void read_parameters(){
    00224                         rv = concat(pf->posterior()._rv(), BMs(0)->posterior()._rv());
    00225                         dim = pf->posterior().dimension() + BMs(0)->posterior().dimension();
    00226                         bdm_assert_debug(dim == rv._dsize(), "Wrong name ");
    00227                 }
    00228                 vec mean() const {
    00229                         const vec &w = pf->posterior()._w();
    00230                         vec pom = zeros ( BMs(0)->posterior ().dimension() );
    00231                         //compute mean of BMs
    00232                         for ( int i = 0; i < w.length(); i++ ) {
    00233                                 pom += BMs ( i )->posterior().mean() * w ( i );
    00234                         }
    00235                         return concat ( pf->posterior().mean(), pom );
    00236                 }
    00237                 vec variance() const {
    00238                         const vec &w = pf->posterior()._w();
    00239
    00240                         vec pom = zeros ( BMs(0)->posterior ().dimension() );
    00241                         vec pom2 = zeros ( BMs(0)->posterior ().dimension() );
    00242                         vec mea;
    00243
    00244                         for ( int i = 0; i < w.length(); i++ ) {
    00245                                 // save current mean
    00246                                 mea = BMs ( i )->posterior().mean();
    00247                                 pom += mea * w ( i );
    00248                                 //compute variance
    00249                                 pom2 += ( BMs ( i )->posterior().variance() + pow ( mea, 2 ) ) * w ( i );
    00250                         }
    00251                         return concat ( pf->posterior().variance(), pom2 - pow ( pom, 2 ) );
    00252                 }
    00253
    00254                 void qbounds ( vec &lb, vec &ub, double perc = 0.95 ) const {
    00255                         //bounds on particles
    00256                         vec lbp;
    00257                         vec ubp;
    00258                         pf->posterior().qbounds ( lbp, ubp );
    00259
    00260                         //bounds on Components
    00261                         int dimC = BMs ( 0 )->posterior().dimension();
    00262                         int j;
    00263                         // temporary
    00264                         vec lbc ( dimC );
    00265                         vec ubc ( dimC );
    00266                         // minima and maxima
    00267                         vec Lbc ( dimC );
    00268                         vec Ubc ( dimC );
    00269                         Lbc = std::numeric_limits<double>::infinity();
    00270                         Ubc = -std::numeric_limits<double>::infinity();
    00271
    00272                         for ( int i = 0; i < BMs.length(); i++ ) {
    00273                                 // check Coms
    00274                                 BMs ( i )->posterior().qbounds ( lbc, ubc );
    00275                                 //save either minima or maxima
    00276                                 for ( j = 0; j < dimC; j++ ) {
    00277                                         if ( lbc ( j ) < Lbc ( j ) ) {
    00278                                                 Lbc ( j ) = lbc ( j );
    00279                                         }
    00280                                         if ( ubc ( j ) > Ubc ( j ) ) {
    00281                                                 Ubc ( j ) = ubc ( j );
    00282                                         }
    00283                                 }
    00284                         }
    00285                         lb = concat ( lbp, Lbc );
    00286                         ub = concat ( ubp, Ubc );
    00287                 }
    00288
    00289                 vec sample() const {
    00290                         bdm_error ( "Not implemented" );
    00291                         return vec();
    00292                 }
    00293
    00294                 double evallog ( const vec &val ) const {
    00295                         bdm_error ( "not implemented" );
    00296                         return 0.0;
    00297                 }
    00298         };
    00299
    00301         mpfepdf jest;
    00302
    00304         bool opt_L_mea;
    00305
    00306 public:
    00308         MPF () :  jest (pf,BMs) {};
    00310         void set_parameters ( shared_ptr<mpdf> par0, shared_ptr<mpdf> obs0, int n0, RESAMPLING_METHOD rm = SYSTEMATIC ) {
    00311                 pf->set_model ( par0, obs0);
    00312                 pf->set_parameters(n0, rm );
    00313                 BMs.set_length ( n0 );
    00314         }
    00316         void set_BM ( const BM &BMcond0 ) {
    00317
    00318                 int n=pf->__w().length();
    00319                 BMs.set_length(n);
    00320                 // copy
    00321                 //BMcond0 .condition ( pf->posterior()._sample ( 0 ) );
    00322                 for ( int i = 0; i < n; i++ ) {
    00323                         BMs ( i ) = BMcond0._copy_();
    00324                 }
    00325         };
    00326
    00327         void bayes ( const vec &dt );
    00328         const epdf& posterior() const {
    00329                 return jest;
    00330         }
    00333         void set_options ( const string &opt ) {
    00334                 BM::set_options(opt);
    00335                 opt_L_mea = ( opt.find ( "logmeans" ) != string::npos );
    00336         }
    00337
    00339         const BM* _BM ( int i ) {
    00340                 return BMs ( i );
    00341         }
    00342
    00354         void from_setting(const Setting &set){
    00355                 shared_ptr<mpdf> par = UI::build<mpdf>(set,"parameter_pdf",UI::compulsory);
    00356                 shared_ptr<mpdf> obs= new mpdf(); // not used!!
    00357
    00358                 pf = new PF;
    00359                 // rpior must be set before BM
    00360                 pf->prior_from_set(set);
    00361                 pf->resmethod_from_set(set);
    00362                 pf->set_model(par,obs);
    00363
    00364                 shared_ptr<BM> BM0 =UI::build<BM>(set,"BM",UI::compulsory);
    00365                 set_BM(*BM0);
    00366
    00367                 string opt;
    00368                 if (UI::get(opt,set,"options",UI::optional)){
    00369                         set_options(opt);
    00370                 }
    00371                 //set drv
    00372                 //find potential input - what remains in rvc when we subtract rv
    00373                 RV u = par->_rvc().remove_time().subt( par->_rv() );
    00374                 set_drv(concat(BM0->_drv(),u) );
    00375                 validate();
    00376         }
    00377         void validate(){
    00378                 try{
    00379                 pf->validate();
    00380                 } catch (std::exception &e){
    00381                         throw UIException("Error in PF part of MPF:");
    00382                 }
    00383                 jest.read_parameters();
    00384                 for ( int i = 0; i < pf->__w().length(); i++ ) {
    00385                         BMs ( i )->condition ( pf->posterior()._sample ( i ) );
    00386                 }
    00387         }
    00388
    00389 };
    00390 UIREGISTER(MPF);
    00391
    00392 }
    00393 #endif // KF_H
    00394
    00395