Changeset 887 for library/bdm

Show
Ignore:
Timestamp:
03/29/10 23:02:03 (14 years ago)
Author:
smidl
Message:

new base for particle filtering

Location:
library/bdm
Files:
4 modified

Legend:

Unmodified
Added
Removed
  • library/bdm/estim/particles.cpp

    r850 r887  
    55using std::endl; 
    66 
    7 void PF::bayes_gensmp ( const vec &ut ) { 
    8         if ( ut.length() > 0 ) { 
    9                 vec cond ( par->dimensionc() ); 
    10                 cond.set_subvector ( par->dimension(), ut ); 
    11 #pragma parallel for 
    12                 for ( int i = 0; i < n; i++ ) { 
    13                         cond.set_subvector ( 0, _samples ( i ) ); 
    14                         _samples ( i ) = par->samplecond ( cond ); 
    15                         lls ( i ) = 0; 
    16                 } 
    17         } else { 
    18 #pragma parallel for 
    19                 for ( int i = 0; i < n; i++ ) { 
    20                         _samples ( i ) = par->samplecond ( _samples ( i ) ); 
    21                         lls ( i ) = 0; 
    22                 } 
    23         } 
    24 } 
    257 
    268void PF::bayes_weights() { 
     
    2911        // compute weights 
    3012        for ( int  i = 0; i < n; i++ ) { 
    31                 _w ( i ) *= exp ( lls ( i ) - mlls ); // multiply w by likelihood 
     13                w ( i ) *= exp ( lls ( i ) - mlls ); // multiply w by likelihood 
    3214        } 
    3315 
    3416        //renormalize 
    35         double sw = sum ( _w ); 
     17        double sw = sum ( w ); 
    3618        if ( !std::isfinite ( sw ) ) { 
    3719                for ( int i = 0; i < n; i++ ) { 
    38                         if ( !std::isfinite ( _w ( i ) ) ) { 
    39                                 _w ( i ) = 0; 
     20                        if ( !std::isfinite ( w ( i ) ) ) { 
     21                                w ( i ) = 0; 
    4022                        } 
    4123                } 
    42                 sw = sum ( _w ); 
     24                sw = sum ( w ); 
    4325                if ( !std::isfinite ( sw ) || sw == 0.0 ) { 
    4426                        bdm_error ( "Particle filter is lost; no particle is good enough." ); 
    4527                } 
    4628        } 
    47         _w /= sw; 
     29        w /= sw; 
    4830} 
    4931 
    5032void PF::bayes ( const vec &yt, const vec &cond ) { 
    51         const vec &ut = cond; //todo 
    52  
    53         int i; 
    54         // generate samples - time step 
    55         bayes_gensmp ( ut ); 
     33 
     34/*      if (auxiliary){ 
     35                ... 
     36        }*/ 
    5637        // weight them - data step 
    57         for ( i = 0; i < n; i++ ) { 
    58                 lls ( i ) += obs->evallogcond ( yt, _samples ( i ) ); //+= because lls may have something from gensmp! 
     38        for (int i = 0; i < n; i++ ) { 
     39                particles(i)->bayes(yt,cond); 
     40                lls ( i ) = particles(i)->_ll(); //+= because lls may have something from gensmp! 
    5941        } 
    6042 
    6143        bayes_weights(); 
    6244 
    63         if ( do_resampling() ) { 
    64                 est.resample ( resmethod ); 
    65         } 
     45        if ( do_resampling() ) { resample ( ); } 
    6646 
    6747} 
     
    7555// } 
    7656 
    77 vec MPF::mpfepdf::mean() const { 
    78         const vec &w = pf->posterior()._w(); 
    79         vec pom = zeros ( BMs ( 0 )->posterior ().dimension() ); 
    80         //compute mean of BMs 
    81         for ( int i = 0; i < w.length(); i++ ) { 
    82                 pom += BMs ( i )->posterior().mean() * w ( i ); 
    83         } 
    84         return concat ( pf->posterior().mean(), pom ); 
    85 } 
    86  
    87 vec MPF::mpfepdf::variance() const { 
    88         const vec &w = pf->posterior()._w(); 
    89  
    90         vec pom = zeros ( BMs ( 0 )->posterior ().dimension() ); 
    91         vec pom2 = zeros ( BMs ( 0 )->posterior ().dimension() ); 
    92         vec mea; 
    93  
    94         for ( int i = 0; i < w.length(); i++ ) { 
    95                 // save current mean 
    96                 mea = BMs ( i )->posterior().mean(); 
    97                 pom += mea * w ( i ); 
    98                 //compute variance 
    99                 pom2 += ( BMs ( i )->posterior().variance() + pow ( mea, 2 ) ) * w ( i ); 
    100         } 
    101         return concat ( pf->posterior().variance(), pom2 - pow ( pom, 2 ) ); 
    102 } 
    103  
    104 void MPF::mpfepdf::qbounds ( vec &lb, vec &ub, double perc ) const { 
    105         //bounds on particles 
    106         vec lbp; 
    107         vec ubp; 
    108         pf->posterior().qbounds ( lbp, ubp ); 
    109  
    110         //bounds on Components 
    111         int dimC = BMs ( 0 )->posterior().dimension(); 
    112         int j; 
    113         // temporary 
    114         vec lbc ( dimC ); 
    115         vec ubc ( dimC ); 
    116         // minima and maxima 
    117         vec Lbc ( dimC ); 
    118         vec Ubc ( dimC ); 
    119         Lbc = std::numeric_limits<double>::infinity(); 
    120         Ubc = -std::numeric_limits<double>::infinity(); 
    121  
    122         for ( int i = 0; i < BMs.length(); i++ ) { 
    123                 // check Coms 
    124                 BMs ( i )->posterior().qbounds ( lbc, ubc ); 
    125                 //save either minima or maxima 
    126                 for ( j = 0; j < dimC; j++ ) { 
    127                         if ( lbc ( j ) < Lbc ( j ) ) { 
    128                                 Lbc ( j ) = lbc ( j ); 
    129                         } 
    130                         if ( ubc ( j ) > Ubc ( j ) ) { 
    131                                 Ubc ( j ) = ubc ( j ); 
    132                         } 
    133                 } 
    134         } 
    135         lb = concat ( lbp, Lbc ); 
    136         ub = concat ( ubp, Ubc ); 
    137 } 
    138  
    139  
    140  
    141 void MPF::bayes ( const vec &yt, const vec &cond ) { 
    142         // follows PF::bayes in most places!!! 
    143         int i; 
    144         int n = pf->__w().length(); 
    145         vec &lls = pf->_lls(); 
    146         Array<vec> &samples = pf->__samples(); 
    147  
    148         // generate samples - time step 
    149         pf->bayes_gensmp ( vec ( 0 ) ); 
    150         // weight them - data step 
    151 #pragma parallel for 
    152         for ( i = 0; i < n; i++ ) { 
    153                 vec bm_cond ( BMs ( i )->dimensionc() ); 
    154                 this2bm.fill_cond ( yt, cond, bm_cond ); 
    155                 pf2bm.filldown ( samples ( i ), bm_cond ); 
    156                 BMs ( i ) -> bayes ( this2bm.pushdown ( yt ), bm_cond ); 
    157                 lls ( i ) += BMs ( i )->_ll(); 
    158         } 
    159  
    160         pf->bayes_weights(); 
    161  
    162         ivec ind; 
    163         if ( pf->do_resampling() ) { 
    164                 pf->resample ( ind ); 
    165  
    166 #pragma omp parallel for 
    167                 for ( i = 0; i < n; i++ ) { 
    168                         if ( ind ( i ) != i ) {//replace the current Bm by a new one 
    169                                 delete BMs ( i ); 
    170                                 BMs ( i ) = (BM*) BMs ( ind ( i ) )->_copy(); //copy constructor 
    171                         } 
    172                 }; 
    173         } 
    174 }; 
    175  
    176  
    177 void MPF_ARXg::bayes ( const vec &yt, const vec &cond ) { 
    178         // follows PF::bayes in most places!!! 
    179         int i; 
    180         int n = pf->__w().length(); 
    181         vec &lls = pf->_lls(); 
    182         Array<vec> &samples = pf->__samples(); 
    183          
    184         // generate samples - time step 
    185         for (int i =0;i<n; i++){ 
    186                 mat M; chmat R; 
    187                 BMsp(i)->posterior().sample_mat(M,R); 
    188                 vec diff=R._Ch().T()*randn(samples(i).length()); 
    189                 samples(i) = g->eval(samples(i)) + diff; 
    190                 //////////// 
    191 //              samples(i) = cond; 
    192                 ///////// 
    193                 BMsp(i)->bayes(diff); 
    194         } 
    195         // weight them - data step 
    196         enorm<ldmat> ooo; 
    197         ooo.set_parameters(zeros(2),0.1*eye(2)); 
    198         ooo.validate(); 
    199          
    200         #pragma parallel for 
    201         for ( i = 0; i < n; i++ ) { 
    202                 vec tmp=yt - h->eval(samples(i)); 
    203                 BMso ( i ) -> bayes ( tmp ); 
    204                 lls ( i ) += BMso ( i )->_ll(); 
    205                 ///// 
    206                 ooo._mu()=h->eval(samples(i)); 
    207                 lls(i) = ooo.evallog_nn(yt); 
    208         } 
    209          
    210         pf->bayes_weights(); 
    211          
    212         ivec ind; 
    213         if ( pf->do_resampling() ) { 
    214                 pf->resample ( ind ); 
    215                  
    216                 #pragma omp parallel for 
    217                 for ( i = 0; i < n; i++ ) { 
    218                         if ( ind ( i ) != i ) {//replace the current Bm by a new one 
    219                                 delete BMso ( i ); 
    220                                 BMso ( i ) = (ARX*) BMso ( ind ( i ) )->_copy(); //copy constructor 
    221                                 delete BMsp ( i ); 
    222                                 BMsp ( i ) = (ARX*) BMsp ( ind ( i ) )->_copy(); //copy constructor 
    223                         } 
    224                 }; 
    225         } 
    226 }; 
    227  
    228 void MPF_ARXg::from_setting(const libconfig::Setting& set) { 
    229         BM::from_setting(set); 
    230          
    231         pf = new PF; 
    232         // prior must be set before BM 
    233         pf->prior_from_set ( set ); 
    234         pf->resmethod_from_set ( set ); 
    235                  
    236         // read functions g and h 
    237         g=UI::build<fnc>(set,"g"); 
    238         h=UI::build<fnc>(set,"h"); 
    239          
    240         set_dim( g->dimension()); 
    241         dimy = h->dimension(); 
    242          
    243         shared_ptr<ARX> arxo=UI::build<ARX>(set,"arxo"); 
    244         shared_ptr<ARX> arxp=UI::build<ARX>(set,"arxp"); 
    245         int n = pf->__samples().length(); 
    246         BMso.set_length(n); 
    247         BMsp.set_length(n); 
    248         for(int i=0; i<n; i++){ 
    249                 BMso(i)=arxo->_copy(); 
    250                 BMsp(i)=arxp->_copy(); 
    251         } 
    252          
    253         yrv = arxo->_yrv(); 
    254         //rvc = arxo->_rvc(); 
    255          
    256         validate(); 
    257          
    258 } 
     57// vec MPF::mpfepdf::mean() const { 
     58//      const vec &w = pf->posterior()._w(); 
     59//      vec pom = zeros ( BMs ( 0 )->posterior ().dimension() ); 
     60//      //compute mean of BMs 
     61//      for ( int i = 0; i < w.length(); i++ ) { 
     62//              pom += BMs ( i )->posterior().mean() * w ( i ); 
     63//      } 
     64//      return concat ( pf->posterior().mean(), pom ); 
     65// } 
     66//  
     67// vec MPF::mpfepdf::variance() const { 
     68//      const vec &w = pf->posterior()._w(); 
     69//  
     70//      vec pom = zeros ( BMs ( 0 )->posterior ().dimension() ); 
     71//      vec pom2 = zeros ( BMs ( 0 )->posterior ().dimension() ); 
     72//      vec mea; 
     73//  
     74//      for ( int i = 0; i < w.length(); i++ ) { 
     75//              // save current mean 
     76//              mea = BMs ( i )->posterior().mean(); 
     77//              pom += mea * w ( i ); 
     78//              //compute variance 
     79//              pom2 += ( BMs ( i )->posterior().variance() + pow ( mea, 2 ) ) * w ( i ); 
     80//      } 
     81//      return concat ( pf->posterior().variance(), pom2 - pow ( pom, 2 ) ); 
     82// } 
     83//  
     84// void MPF::mpfepdf::qbounds ( vec &lb, vec &ub, double perc ) const { 
     85//      //bounds on particles 
     86//      vec lbp; 
     87//      vec ubp; 
     88//      pf->posterior().qbounds ( lbp, ubp ); 
     89//  
     90//      //bounds on Components 
     91//      int dimC = BMs ( 0 )->posterior().dimension(); 
     92//      int j; 
     93//      // temporary 
     94//      vec lbc ( dimC ); 
     95//      vec ubc ( dimC ); 
     96//      // minima and maxima 
     97//      vec Lbc ( dimC ); 
     98//      vec Ubc ( dimC ); 
     99//      Lbc = std::numeric_limits<double>::infinity(); 
     100//      Ubc = -std::numeric_limits<double>::infinity(); 
     101//  
     102//      for ( int i = 0; i < BMs.length(); i++ ) { 
     103//              // check Coms 
     104//              BMs ( i )->posterior().qbounds ( lbc, ubc ); 
     105//              //save either minima or maxima 
     106//              for ( j = 0; j < dimC; j++ ) { 
     107//                      if ( lbc ( j ) < Lbc ( j ) ) { 
     108//                              Lbc ( j ) = lbc ( j ); 
     109//                      } 
     110//                      if ( ubc ( j ) > Ubc ( j ) ) { 
     111//                              Ubc ( j ) = ubc ( j ); 
     112//                      } 
     113//              } 
     114//      } 
     115//      lb = concat ( lbp, Lbc ); 
     116//      ub = concat ( ubp, Ubc ); 
     117// } 
     118//  
     119//  
     120//  
     121// void MPF::bayes ( const vec &yt, const vec &cond ) { 
     122//      // follows PF::bayes in most places!!! 
     123//      int i; 
     124//      int n = pf->__w().length(); 
     125//      vec &lls = pf->_lls(); 
     126//      Array<vec> &samples = pf->__samples(); 
     127//  
     128//      // generate samples - time step 
     129//      pf->bayes_gensmp ( vec ( 0 ) ); 
     130//      // weight them - data step 
     131// #pragma parallel for 
     132//      for ( i = 0; i < n; i++ ) { 
     133//              vec bm_cond ( BMs ( i )->dimensionc() ); 
     134//              this2bm.fill_cond ( yt, cond, bm_cond ); 
     135//              pf2bm.filldown ( samples ( i ), bm_cond ); 
     136//              BMs ( i ) -> bayes ( this2bm.pushdown ( yt ), bm_cond ); 
     137//              lls ( i ) += BMs ( i )->_ll(); 
     138//      } 
     139//  
     140//      pf->bayes_weights(); 
     141//  
     142//      ivec ind; 
     143//      if ( pf->do_resampling() ) { 
     144//              pf->resample ( ind ); 
     145//  
     146// #pragma omp parallel for 
     147//              for ( i = 0; i < n; i++ ) { 
     148//                      if ( ind ( i ) != i ) {//replace the current Bm by a new one 
     149//                              delete BMs ( i ); 
     150//                              BMs ( i ) = (BM*) BMs ( ind ( i ) )->_copy(); //copy constructor 
     151//                      } 
     152//              }; 
     153//      } 
     154// }; 
     155//  
     156//  
     157// void MPF_ARXg::bayes ( const vec &yt, const vec &cond ) { 
     158//      // follows PF::bayes in most places!!! 
     159//      int i; 
     160//      int n = pf->__w().length(); 
     161//      vec &lls = pf->_lls(); 
     162//      Array<vec> &samples = pf->__samples(); 
     163//       
     164//      // generate samples - time step 
     165//      for (int i =0;i<n; i++){ 
     166//              mat M; chmat R; 
     167//              BMsp(i)->posterior().sample_mat(M,R); 
     168//              vec diff=R._Ch().T()*randn(samples(i).length()); 
     169//              samples(i) = g->eval(samples(i)) + diff; 
     170//              //////////// 
     171// //           samples(i) = cond; 
     172//              ///////// 
     173//              BMsp(i)->bayes(diff); 
     174//      } 
     175//      // weight them - data step 
     176//      enorm<ldmat> ooo; 
     177//      ooo.set_parameters(zeros(2),0.1*eye(2)); 
     178//      ooo.validate(); 
     179//       
     180//      #pragma parallel for 
     181//      for ( i = 0; i < n; i++ ) { 
     182//              vec tmp=yt - h->eval(samples(i)); 
     183//              BMso ( i ) -> bayes ( tmp ); 
     184//              lls ( i ) += BMso ( i )->_ll(); 
     185//              ///// 
     186//              ooo._mu()=h->eval(samples(i)); 
     187//              lls(i) = ooo.evallog_nn(yt); 
     188//      } 
     189//       
     190//      pf->bayes_weights(); 
     191//       
     192//      ivec ind; 
     193//      if ( pf->do_resampling() ) { 
     194//              pf->resample ( ind ); 
     195//               
     196//              #pragma omp parallel for 
     197//              for ( i = 0; i < n; i++ ) { 
     198//                      if ( ind ( i ) != i ) {//replace the current Bm by a new one 
     199//                              delete BMso ( i ); 
     200//                              BMso ( i ) = (ARX*) BMso ( ind ( i ) )->_copy(); //copy constructor 
     201//                              delete BMsp ( i ); 
     202//                              BMsp ( i ) = (ARX*) BMsp ( ind ( i ) )->_copy(); //copy constructor 
     203//                      } 
     204//              }; 
     205//      } 
     206// }; 
     207//  
     208// void MPF_ARXg::from_setting(const libconfig::Setting& set) { 
     209//      BM::from_setting(set); 
     210//       
     211//      pf = new PF; 
     212//      // prior must be set before BM 
     213//      pf->prior_from_set ( set ); 
     214//      pf->resmethod_from_set ( set ); 
     215//               
     216//      // read functions g and h 
     217//      g=UI::build<fnc>(set,"g"); 
     218//      h=UI::build<fnc>(set,"h"); 
     219//       
     220//      set_dim( g->dimension()); 
     221//      dimy = h->dimension(); 
     222//       
     223//      shared_ptr<ARX> arxo=UI::build<ARX>(set,"arxo"); 
     224//      shared_ptr<ARX> arxp=UI::build<ARX>(set,"arxp"); 
     225//      int n = pf->__samples().length(); 
     226//      BMso.set_length(n); 
     227//      BMsp.set_length(n); 
     228//      for(int i=0; i<n; i++){ 
     229//              BMso(i)=arxo->_copy(); 
     230//              BMsp(i)=arxp->_copy(); 
     231//      } 
     232//       
     233//      yrv = arxo->_yrv(); 
     234//      //rvc = arxo->_rvc(); 
     235//       
     236//      validate(); 
     237//       
     238// } 
    259239 
    260240 
  • library/bdm/estim/particles.h

    r870 r887  
    1616 
    1717#include "../estim/arx_ext.h" 
     18#include "../stat/emix.h" 
    1819 
    1920namespace bdm { 
     21         
     22//! class used in PF 
     23class MarginalizedParticle : public BM{ 
     24        protected: 
     25        //! discrte particle 
     26        dirac est_emp; 
     27        //! internal Bayes Model 
     28        shared_ptr<BM> bm;  
     29        //! pdf with for transitional par 
     30        shared_ptr<pdf> par; // pdf for non-linear part 
     31        //! link from this to bm 
     32        shared_ptr<datalink_part> cond2bm; 
     33        //! link from cond to par 
     34        shared_ptr<datalink_part> cond2par; 
     35        //! link from emp 2 par 
     36        shared_ptr<datalink_part> emp2bm; 
     37        //! link from emp 2 par 
     38        shared_ptr<datalink_part> emp2par; 
     39         
     40        //! custom posterior - product 
     41        class eprod_2:public eprod_base { 
     42                protected: 
     43                MarginalizedParticle &mp; 
     44                public: 
     45                eprod_2(MarginalizedParticle &m):mp(m){} 
     46                const epdf* factor(int i) const {return (i==0) ? &mp.bm->posterior() : &mp.est_emp;} 
     47                const int no_factors() const {return 2;} 
     48        } est; 
     49         
     50        public: 
     51                MarginalizedParticle():est(*this){}; 
     52                BM* _copy() const{return new MarginalizedParticle(*this);}; 
     53                void bayes(const vec &dt, const vec &cond){ 
     54                        vec par_cond(par->dimensionc()); 
     55                        cond2par->filldown(cond,par_cond); // copy ut 
     56                        emp2par->filldown(est_emp._point(),par_cond); // copy xt-1 
     57                         
     58                        //sample new particle 
     59                        est_emp.set_point(par->samplecond(par_cond)); 
     60                        //if (evalll) 
     61                        vec bm_cond(bm->dimensionc()); 
     62                        cond2bm->filldown(cond, bm_cond);// set e.g. ut 
     63                        emp2bm->filldown(est_emp._point(), bm_cond);// set e.g. ut 
     64                        bm->bayes(dt, bm_cond); 
     65                        ll=bm->_ll(); 
     66                } 
     67                const eprod_2& posterior() const {return est;} 
     68         
     69                void set_prior(const epdf *pdf0){ 
     70                        const eprod *ep=dynamic_cast<const eprod*>(pdf0); 
     71                        if (ep){ // full prior 
     72                                bdm_assert(ep->no_factors()==2,"Incompatible prod"); 
     73                                bm->set_prior(ep->factor(0)); 
     74                                est_emp.set_point(ep->factor(1)->sample()); 
     75                        } else { 
     76                                // assume prior is only for emp; 
     77                                est_emp.set_point(pdf0->sample()); 
     78                        } 
     79                } 
     80 
     81                /*! parse structure 
     82                \code 
     83                class = "BootstrapParticle"; 
     84                parameter_pdf = {class = 'epdf_offspring', ...}; 
     85                observation_pdf = {class = 'epdf_offspring',...}; 
     86                \endcode 
     87                If rvs are set, then it checks for compatibility. 
     88                */ 
     89                void from_setting(const Setting &set){ 
     90                        BM::from_setting ( set ); 
     91                        par = UI::build<pdf> ( set, "parameter_pdf", UI::compulsory ); 
     92                        bm = UI::build<BM> ( set, "bm", UI::compulsory ); 
     93                } 
     94                void validate(){ 
     95                        yrv = bm->_rv(); 
     96                        dimy = bm->dimension(); 
     97                        set_rv( par->_rv()); 
     98                        set_dim( par->dimension()); 
     99 
     100                        rvc = par->_rvc().subt(par->_rv().copy_t(-1)); 
     101                        rvc.add(bm->_rvc()); // 
     102                         
     103                        cond2bm=new datalink_part; 
     104                        cond2par=new datalink_part; 
     105                        emp2bm  =new datalink_part; 
     106                        emp2par =new datalink_part; 
     107                        cond2bm->set_connection(bm->_rvc(), rvc); 
     108                        cond2par->set_connection(par->_rvc(), rvc); 
     109                        emp2bm->set_connection(bm->_rvc(), _rv()); 
     110                        emp2par->set_connection(par->_rvc(), _rv().copy_t(-1)); 
     111                         
     112                        dimc = rvc._dsize(); 
     113                }; 
     114}; 
     115UIREGISTER(MarginalizedParticle); 
     116 
     117//! class used in PF 
     118class BootstrapParticle : public BM{ 
     119        dirac est; 
     120        shared_ptr<pdf> par; 
     121        shared_ptr<pdf> obs; 
     122        shared_ptr<datalink_part> cond2par; 
     123        shared_ptr<datalink_part> cond2obs; 
     124        shared_ptr<datalink_part> xt2obs; 
     125        shared_ptr<datalink_part> xtm2par; 
     126        public: 
     127                BM* _copy() const{return new BootstrapParticle(*this);}; 
     128                void bayes(const vec &dt, const vec &cond){ 
     129                        vec par_cond(par->dimensionc()); 
     130                        cond2par->filldown(cond,par_cond); // copy ut 
     131                        xtm2par->filldown(est._point(),par_cond); // copy xt-1 
     132                         
     133                        //sample new particle 
     134                        est.set_point(par->samplecond(par_cond)); 
     135                        //if (evalll) 
     136                        vec obs_cond(obs->dimensionc()); 
     137                        cond2obs->filldown(cond, obs_cond);// set e.g. ut 
     138                        xt2obs->filldown(est._point(), obs_cond);// set e.g. ut 
     139                        ll=obs->evallogcond(dt,obs_cond); 
     140                } 
     141                const dirac& posterior() const {return est;} 
     142                 
     143                void set_prior(const epdf *pdf0){est.set_point(pdf0->sample());} 
     144                 
     145                /*! parse structure 
     146                \code 
     147                class = "BootstrapParticle"; 
     148                parameter_pdf = {class = 'epdf_offspring', ...}; 
     149                observation_pdf = {class = 'epdf_offspring',...}; 
     150                \endcode 
     151                If rvs are set, then it checks for compatibility. 
     152                */ 
     153                void from_setting(const Setting &set){ 
     154                        BM::from_setting ( set ); 
     155                        par = UI::build<pdf> ( set, "parameter_pdf", UI::compulsory ); 
     156                        obs = UI::build<pdf> ( set, "observation_pdf", UI::compulsory ); 
     157                } 
     158                void validate(){ 
     159                        yrv = obs->_rv(); 
     160                        dimy = obs->dimension(); 
     161                        set_rv( par->_rv()); 
     162                        set_dim( par->dimension()); 
     163                         
     164                        rvc = par->_rvc().subt(par->_rv().copy_t(-1)); 
     165                        rvc.add(obs->_rvc()); // 
     166                         
     167                        cond2obs=new datalink_part; 
     168                        cond2par=new datalink_part; 
     169                        xt2obs  =new datalink_part; 
     170                        xtm2par =new datalink_part; 
     171                        cond2obs->set_connection(obs->_rvc(), rvc); 
     172                        cond2par->set_connection(par->_rvc(), rvc); 
     173                        xt2obs->set_connection(obs->_rvc(), _rv()); 
     174                        xtm2par->set_connection(par->_rvc(), _rv().copy_t(-1)); 
     175                         
     176                        dimc = rvc._dsize(); 
     177                }; 
     178}; 
     179UIREGISTER(BootstrapParticle); 
     180 
    20181 
    21182/*! 
     
    33194        LOG_LEVEL(PF,weights,samples); 
    34195         
     196        class pf_mix: public emix_base{ 
     197                Array<BM*> &bms; 
     198                public: 
     199                        pf_mix(vec &w0, Array<BM*> &bms0):emix_base(w0),bms(bms0){} 
     200                        const epdf* component(const int &i)const{return &(bms(i)->posterior());} 
     201                        int no_coms() const {return bms.length();} 
     202        }; 
    35203protected: 
    36204        //!number of particles; 
    37205        int n; 
    38206        //!posterior density 
    39         eEmp est; 
    40         //! pointer into \c eEmp 
    41         vec &_w; 
    42         //! pointer into \c eEmp 
    43         Array<vec> &_samples; 
    44         //! Parameter evolution model 
    45         shared_ptr<pdf> par; 
    46         //! Observation model 
    47         shared_ptr<pdf> obs; 
     207        pf_mix est; 
     208        //! weights; 
     209        vec w; 
     210        //! particles 
     211        Array<BM*> particles; 
    48212        //! internal structure storing loglikelihood of predictions 
    49213        vec lls; 
     
    62226        //! \name Constructors 
    63227        //!@{ 
    64         PF ( ) : est(), _w ( est._w() ), _samples ( est._samples() ) { 
    65         }; 
     228        PF ( ) : est(w,particles) { }; 
    66229 
    67230        void set_parameters ( int n0, double res_th0 = 0.5, RESAMPLING_METHOD rm = SYSTEMATIC ) { 
     
    70233                resmethod = rm; 
    71234        }; 
    72         void set_model ( shared_ptr<pdf> par0, shared_ptr<pdf> obs0 ) { 
    73                 par = par0; 
    74                 obs = obs0; 
     235        void set_model ( const BM *particle0, const epdf *prior) { 
     236                if (n>0){ 
     237                        particles.set_length(n); 
     238                        for (int i=0; i<n;i++){ 
     239                                particles(i) = particle0->_copy(); 
     240                                particles(i)->set_prior(prior); 
     241                        } 
     242                } 
    75243                // set values for posterior 
    76                 est.set_rv ( par->_rv() ); 
     244                est.set_rv ( particle0->posterior()._rv() ); 
    77245        }; 
    78246        void set_statistics ( const vec w0, const epdf &epdf0 ) { 
    79                 est.set_statistics ( w0, epdf0 ); 
     247                //est.set_statistics ( w0, epdf0 ); 
    80248        }; 
    81         void set_statistics ( const eEmp &epdf0 ) { 
     249/*      void set_statistics ( const eEmp &epdf0 ) { 
    82250                bdm_assert_debug ( epdf0._rv().equal ( par->_rv() ), "Incompatible input" ); 
    83251                est = epdf0; 
    84         }; 
     252        };*/ 
    85253        //!@} 
    86254 
    87         //! Set posterior density by sampling from epdf0 
    88         //! bayes I - generate samples and add their weights to lls 
    89         virtual void bayes_gensmp ( const vec &ut ); 
    90         //! bayes II - compute weights of the 
     255        //! bayes compute weights of the 
    91256        virtual void bayes_weights(); 
    92257        //! important part of particle filtering - decide if it is time to perform resampling 
    93258        virtual bool do_resampling() { 
    94                 double eff = 1.0 / ( _w * _w ); 
     259                double eff = 1.0 / ( w * w ); 
    95260                return eff < ( res_threshold*n ); 
    96261        } 
    97262        void bayes ( const vec &yt, const vec &cond ); 
    98         //!access function 
    99         vec& __w() { 
    100                 return _w; 
    101         } 
    102263        //!access function 
    103264        vec& _lls() { 
     
    109270        } 
    110271        //! return correctly typed posterior (covariant return) 
    111         const eEmp& posterior() const { 
     272        const pf_mix& posterior() const { 
    112273                return est; 
    113274        } 
     
    127288        void from_setting ( const Setting &set ) { 
    128289                BM::from_setting ( set ); 
    129                 par = UI::build<pdf> ( set, "parameter_pdf", UI::compulsory ); 
    130                 obs = UI::build<pdf> ( set, "observation_pdf", UI::compulsory ); 
    131  
    132                 prior_from_set ( set ); 
     290                 
     291                shared_ptr<BM> bm0 = UI::build<BM>(set, "particle",UI::compulsory); 
     292                 
     293                shared_ptr<epdf> pri = UI::build<epdf> ( set, "prior", UI::compulsory ); 
     294                n =0; 
     295                UI::get(n,set,"n",UI::optional);; 
     296                if (n>0){ 
     297                        particles.set_length(n); 
     298                        for(int i=0;i<n;i++){particles(i)=bm0->_copy();} 
     299                        w = ones(n)/n; 
     300                } 
     301                set_prior(pri.get()); 
     302                // set resampling method 
    133303                resmethod_from_set ( set ); 
    134                 // set resampling method 
    135304                //set drv 
    136                 //find potential input - what remains in rvc when we subtract rv 
    137                 RV u = par->_rvc().remove_time().subt ( par->_rv() ); 
    138                 //find potential input - what remains in rvc when we subtract x_t 
    139                 RV obs_u = obs->_rvc().remove_time().subt ( par->_rv() ); 
    140  
    141                 u.add ( obs_u ); // join both u, and check if they do not overlap 
    142  
    143                 set_yrv ( obs->_rv() ); 
    144                 rvc = u; 
     305 
     306                set_yrv ( bm0->_rv() ); 
     307                rvc = bm0->_rvc(); 
     308                BM::set_rv(bm0->_rv()); 
     309                yrv=bm0->_yrv(); 
     310        } 
     311         
     312        void set_prior(const epdf *pri){ 
     313                const emix_base *emi=dynamic_cast<const emix_base*>(pri); 
     314                if (emi) { 
     315                        bdm_assert(particles.length()>0, "initial particle is not assigned"); 
     316                        n = emi->_w().length(); 
     317                        int old_n = particles.length(); 
     318                        if (n!=old_n){ 
     319                                particles.set_length(n,true); 
     320                        }  
     321                        for(int i=old_n;i<n;i++){particles(i)=particles(0)->_copy();} 
     322                         
     323                        for (int i =0; i<n; i++){ 
     324                                particles(i)->set_prior(emi->_com(i)); 
     325                        } 
     326                } else { 
     327                        // try to find "n" 
     328                        bdm_assert(n>0, "Field 'n' must be filled when prior is not of type emix"); 
     329                        for (int i =0; i<n; i++){ 
     330                                particles(i)->set_prior(pri); 
     331                        } 
     332                         
     333                } 
    145334        } 
    146335        //! auxiliary function reading parameter 'resmethod' from configuration file 
     
    165354                }; 
    166355                if ( !UI::get ( res_threshold, set, "res_threshold", UI::optional ) ) { 
    167                         res_threshold = 0.5; 
     356                        res_threshold = 0.9; 
    168357                } 
    169358                //validate(); 
    170         } 
    171         //! load prior information from set and set internal structures accordingly 
    172         void prior_from_set ( const Setting & set ) { 
    173                 shared_ptr<epdf> pri = UI::build<epdf> ( set, "prior", UI::compulsory ); 
    174  
    175                 eEmp *test_emp = dynamic_cast<eEmp*> ( & ( *pri ) ); 
    176                 if ( test_emp ) { // given pdf is sampled 
    177                         est = *test_emp; 
    178                 } else { 
    179                         //int n; 
    180                         if ( !UI::get ( n, set, "n", UI::optional ) ) { 
    181                                 n = 10; 
    182                         } 
    183                         // sample from prior 
    184                         set_statistics ( ones ( n ) / n, *pri ); 
    185                 } 
    186                 n = est._w().length(); 
    187                 lls=zeros(n); 
    188359        } 
    189360 
    190361        void validate() { 
    191362                BM::validate(); 
    192                 bdm_assert ( par, "PF::parameter_pdf not specified." ); 
    193                 n = _w.length(); 
     363                est.validate(); 
     364                bdm_assert ( n>0, "empty particle pool" ); 
     365                n = w.length(); 
    194366                lls = zeros ( n ); 
    195367 
    196                 if ( par->_rv()._dsize() > 0 ) { 
    197                         bdm_assert ( par->_rv()._dsize() == est.dimension(), "Mismatch of RV and dimension of posterior" ); 
     368                if ( particles(0)->_rv()._dsize() > 0 ) { 
     369                        bdm_assert (  particles(0)->_rv()._dsize() == est.dimension(), "Mismatch of RV and dimension of posterior" ); 
    198370                } 
    199371        } 
    200372        //! resample posterior density (from outside - see MPF) 
    201         void resample ( ivec &ind ) { 
    202                 est.resample ( ind, resmethod ); 
     373        void resample ( ) { 
     374                ivec ind = zeros_i ( n ); 
     375                bdm::resample(w,ind,resmethod); 
     376                // copy the internals according to ind 
     377                for (int i = 0; i < n; i++ ) { 
     378                        if ( ind ( i ) != i ) { 
     379                                particles( i ) = particles( ind ( i ) )->_copy(); 
     380                        } 
     381                        w ( i ) = 1.0 / n; 
     382                } 
    203383        } 
    204384        //! access function 
    205         Array<vec>& __samples() { 
    206                 return _samples; 
    207         } 
    208  
    209         virtual double logpred ( const vec &yt ) const NOT_IMPLEMENTED(0); 
    210  
    211         virtual epdf* epredictor() const NOT_IMPLEMENTED(NULL); 
    212          
    213         virtual pdf* predictor() const NOT_IMPLEMENTED(NULL); 
     385        Array<BM*>& _particles() { 
     386                return particles; 
     387        } 
     388 
    214389}; 
    215390UIREGISTER ( PF ); 
     
    222397*/ 
    223398 
    224 class MPF : public BM  { 
    225         //! \var log_level_enums means  
    226         //! TODO DOPLNIT 
    227         LOG_LEVEL(MPF,means); 
    228  
    229 protected: 
    230         //! particle filter on non-linear variable 
    231         shared_ptr<PF> pf; 
    232         //! Array of Bayesian models 
    233         Array<BM*> BMs; 
    234  
    235         //! internal class for pdf providing composition of eEmp with external components 
    236  
    237         class mpfepdf : public epdf  { 
    238                 //! pointer to particle filter 
    239                 shared_ptr<PF> &pf; 
    240                 //! pointer to Array of BMs 
    241                 Array<BM*> &BMs; 
    242         public: 
    243                 //! constructor 
    244                 mpfepdf ( shared_ptr<PF> &pf0, Array<BM*> &BMs0 ) : epdf(), pf ( pf0 ), BMs ( BMs0 ) { }; 
    245                 //! a variant of set parameters - this time, parameters are read from BMs and pf 
    246                 void read_parameters() { 
    247                         rv = concat ( pf->posterior()._rv(), BMs ( 0 )->posterior()._rv() ); 
    248                         dim = pf->posterior().dimension() + BMs ( 0 )->posterior().dimension(); 
    249                         bdm_assert_debug ( dim == rv._dsize(), "Wrong name " ); 
    250                 } 
    251                 vec mean() const; 
    252  
    253                 vec variance() const; 
    254  
    255                 void qbounds ( vec &lb, vec &ub, double perc = 0.95 ) const; 
    256  
    257                 vec sample() const NOT_IMPLEMENTED(0); 
    258  
    259                 double evallog ( const vec &val ) const NOT_IMPLEMENTED(0);              
    260         }; 
    261  
    262         //! Density joining PF.est with conditional parts 
    263         mpfepdf jest; 
    264  
    265         //! datalink from global yt and cond (Up) to BMs yt and cond (Down) 
    266         datalink_m2m this2bm; 
    267         //! datalink from global yt and cond (Up) to PFs yt and cond (Down) 
    268         datalink_m2m this2pf; 
    269         //!datalink from PF part to BM 
    270         datalink_part pf2bm; 
    271  
    272 public: 
    273         //! Default constructor. 
    274         MPF () :  jest ( pf, BMs ) {}; 
    275         //! set all parameters at once 
    276         void set_pf ( shared_ptr<pdf> par0, int n0, RESAMPLING_METHOD rm = SYSTEMATIC ) { 
    277                 if (!pf) pf=new PF; 
    278                 pf->set_model ( par0, par0 ); // <=== nasty!!! 
    279                 pf->set_parameters ( n0, rm ); 
    280                 pf->set_rv(par0->_rv()); 
    281                 BMs.set_length ( n0 ); 
    282         } 
    283         //! set a prototype of BM, copy it to as many times as there is particles in pf 
    284         void set_BM ( const BM &BMcond0 ) { 
    285  
    286                 int n = pf->__w().length(); 
    287                 BMs.set_length ( n ); 
    288                 // copy 
    289                 //BMcond0 .condition ( pf->posterior()._sample ( 0 ) ); 
    290                 for ( int i = 0; i < n; i++ ) { 
    291                         BMs ( i ) = (BM*) BMcond0._copy(); 
    292                 } 
    293         }; 
    294  
    295         void bayes ( const vec &yt, const vec &cond ); 
    296  
    297         const epdf& posterior() const { 
    298                 return jest; 
    299         } 
    300  
    301         //!Access function 
    302         const BM* _BM ( int i ) { 
    303                 return BMs ( i ); 
    304         } 
    305         PF& _pf() {return *pf;} 
    306  
    307  
    308         virtual double logpred ( const vec &yt ) const NOT_IMPLEMENTED(0);  
    309                  
    310         virtual epdf* epredictor() const NOT_IMPLEMENTED(NULL);  
    311          
    312         virtual pdf* predictor() const NOT_IMPLEMENTED(NULL);  
    313  
    314  
    315         /*! configuration structure for basic PF 
    316         \code 
    317         BM              = BM_class;           // Bayesian filtr for analytical part of the model 
    318         parameter_pdf   = pdf_class;         // transitional pdf for non-parametric part of the model 
    319         prior           = epdf_class;         // prior probability density 
    320         --- optional --- 
    321         n               = 10;                 // number of particles 
    322         resmethod       = 'systematic', or 'multinomial', or 'stratified' 
    323                                                                                   // resampling method 
    324         \endcode 
    325         */ 
    326         void from_setting ( const Setting &set ) { 
    327                 BM::from_setting( set ); 
    328  
    329                 shared_ptr<pdf> par = UI::build<pdf> ( set, "parameter_pdf", UI::compulsory ); 
    330  
    331                 pf = new PF; 
    332                 // prior must be set before BM 
    333                 pf->prior_from_set ( set ); 
    334                 pf->resmethod_from_set ( set ); 
    335                 pf->set_model ( par, par ); // too hackish! 
    336  
    337                 shared_ptr<BM> BM0 = UI::build<BM> ( set, "BM", UI::compulsory ); 
    338                 set_BM ( *BM0 ); 
    339  
    340                 //set drv 
    341                 //??set_yrv(concat(BM0->_yrv(),u) ); 
    342                 set_yrv ( BM0->_yrv() ); 
    343                 rvc = BM0->_rvc().subt ( par->_rv() ); 
    344                 //find potential input - what remains in rvc when we subtract rv 
    345                 RV u = par->_rvc().subt ( par->_rv().copy_t ( -1 ) ); 
    346                 rvc.add ( u ); 
    347                 dimc = rvc._dsize(); 
    348                 validate(); 
    349         } 
    350  
    351         void validate() { 
    352                 BM::validate(); 
    353                 try { 
    354                         pf->validate(); 
    355                 } catch ( std::exception ) { 
    356                         throw UIException ( "Error in PF part of MPF:" ); 
    357                 } 
    358                 jest.read_parameters(); 
    359                 this2bm.set_connection ( BMs ( 0 )->_yrv(), BMs ( 0 )->_rvc(), yrv, rvc ); 
    360                 this2pf.set_connection ( pf->_yrv(), pf->_rvc(), yrv, rvc ); 
    361                 pf2bm.set_connection ( BMs ( 0 )->_rvc(), pf->posterior()._rv() ); 
    362         } 
    363 }; 
    364 UIREGISTER ( MPF ); 
     399// class MPF : public BM  { 
     400//      //! Introduces new option  
     401//      //! \li means - meaning TODO 
     402//      LOG_LEVEL(MPF,means); 
     403// protected: 
     404//      //! particle filter on non-linear variable 
     405//      shared_ptr<PF> pf; 
     406//      //! Array of Bayesian models 
     407//      Array<BM*> BMs; 
     408//  
     409//      //! internal class for pdf providing composition of eEmp with external components 
     410//  
     411//      class mpfepdf : public epdf  { 
     412//              //! pointer to particle filter 
     413//              shared_ptr<PF> &pf; 
     414//              //! pointer to Array of BMs 
     415//              Array<BM*> &BMs; 
     416//      public: 
     417//              //! constructor 
     418//              mpfepdf ( shared_ptr<PF> &pf0, Array<BM*> &BMs0 ) : epdf(), pf ( pf0 ), BMs ( BMs0 ) { }; 
     419//              //! a variant of set parameters - this time, parameters are read from BMs and pf 
     420//              void read_parameters() { 
     421//                      rv = concat ( pf->posterior()._rv(), BMs ( 0 )->posterior()._rv() ); 
     422//                      dim = pf->posterior().dimension() + BMs ( 0 )->posterior().dimension(); 
     423//                      bdm_assert_debug ( dim == rv._dsize(), "Wrong name " ); 
     424//              } 
     425//              vec mean() const; 
     426//  
     427//              vec variance() const; 
     428//  
     429//              void qbounds ( vec &lb, vec &ub, double perc = 0.95 ) const; 
     430//  
     431//              vec sample() const NOT_IMPLEMENTED(0); 
     432//  
     433//              double evallog ( const vec &val ) const NOT_IMPLEMENTED(0);              
     434//      }; 
     435//  
     436//      //! Density joining PF.est with conditional parts 
     437//      mpfepdf jest; 
     438//  
     439//      //! datalink from global yt and cond (Up) to BMs yt and cond (Down) 
     440//      datalink_m2m this2bm; 
     441//      //! datalink from global yt and cond (Up) to PFs yt and cond (Down) 
     442//      datalink_m2m this2pf; 
     443//      //!datalink from PF part to BM 
     444//      datalink_part pf2bm; 
     445//  
     446// public: 
     447//      //! Default constructor. 
     448//      MPF () :  jest ( pf, BMs ) {}; 
     449//      //! set all parameters at once 
     450//      void set_pf ( shared_ptr<pdf> par0, int n0, RESAMPLING_METHOD rm = SYSTEMATIC ) { 
     451//              if (!pf) pf=new PF; 
     452//              pf->set_model ( par0, par0 ); // <=== nasty!!! 
     453//              pf->set_parameters ( n0, rm ); 
     454//              pf->set_rv(par0->_rv()); 
     455//              BMs.set_length ( n0 ); 
     456//      } 
     457//      //! set a prototype of BM, copy it to as many times as there is particles in pf 
     458//      void set_BM ( const BM &BMcond0 ) { 
     459//  
     460//              int n = pf->__w().length(); 
     461//              BMs.set_length ( n ); 
     462//              // copy 
     463//              //BMcond0 .condition ( pf->posterior()._sample ( 0 ) ); 
     464//              for ( int i = 0; i < n; i++ ) { 
     465//                      BMs ( i ) = (BM*) BMcond0._copy(); 
     466//              } 
     467//      }; 
     468//  
     469//      void bayes ( const vec &yt, const vec &cond ); 
     470//  
     471//      const epdf& posterior() const { 
     472//              return jest; 
     473//      } 
     474//  
     475//      //!Access function 
     476//      const BM* _BM ( int i ) { 
     477//              return BMs ( i ); 
     478//      } 
     479//      PF& _pf() {return *pf;} 
     480//  
     481//  
     482//      virtual double logpred ( const vec &yt ) const NOT_IMPLEMENTED(0);  
     483//               
     484//      virtual epdf* epredictor() const NOT_IMPLEMENTED(NULL);  
     485//       
     486//      virtual pdf* predictor() const NOT_IMPLEMENTED(NULL);  
     487//  
     488//  
     489//      /*! configuration structure for basic PF 
     490//      \code 
     491//      BM              = BM_class;           // Bayesian filtr for analytical part of the model 
     492//      parameter_pdf   = pdf_class;         // transitional pdf for non-parametric part of the model 
     493//      prior           = epdf_class;         // prior probability density 
     494//      --- optional --- 
     495//      n               = 10;                 // number of particles 
     496//      resmethod       = 'systematic', or 'multinomial', or 'stratified' 
     497//                                                                                // resampling method 
     498//      \endcode 
     499//      */ 
     500//      void from_setting ( const Setting &set ) { 
     501//              BM::from_setting( set ); 
     502//  
     503//              shared_ptr<pdf> par = UI::build<pdf> ( set, "parameter_pdf", UI::compulsory ); 
     504//  
     505//              pf = new PF; 
     506//              // prior must be set before BM 
     507//              pf->prior_from_set ( set ); 
     508//              pf->resmethod_from_set ( set ); 
     509//              pf->set_model ( par, par ); // too hackish! 
     510//  
     511//              shared_ptr<BM> BM0 = UI::build<BM> ( set, "BM", UI::compulsory ); 
     512//              set_BM ( *BM0 ); 
     513//  
     514//              //set drv 
     515//              //??set_yrv(concat(BM0->_yrv(),u) ); 
     516//              set_yrv ( BM0->_yrv() ); 
     517//              rvc = BM0->_rvc().subt ( par->_rv() ); 
     518//              //find potential input - what remains in rvc when we subtract rv 
     519//              RV u = par->_rvc().subt ( par->_rv().copy_t ( -1 ) ); 
     520//              rvc.add ( u ); 
     521//              dimc = rvc._dsize(); 
     522//              validate(); 
     523//      } 
     524//  
     525//      void validate() { 
     526//              BM::validate(); 
     527//              try { 
     528//                      pf->validate(); 
     529//              } catch ( std::exception ) { 
     530//                      throw UIException ( "Error in PF part of MPF:" ); 
     531//              } 
     532//              jest.read_parameters(); 
     533//              this2bm.set_connection ( BMs ( 0 )->_yrv(), BMs ( 0 )->_rvc(), yrv, rvc ); 
     534//              this2pf.set_connection ( pf->_yrv(), pf->_rvc(), yrv, rvc ); 
     535//              pf2bm.set_connection ( BMs ( 0 )->_rvc(), pf->posterior()._rv() ); 
     536//      } 
     537// }; 
     538// UIREGISTER ( MPF ); 
    365539 
    366540/*! ARXg for estimation of state-space variances 
    367541*/ 
    368 class MPF_ARXg :public BM{ 
    369         protected: 
    370         shared_ptr<PF> pf; 
    371         //! pointer to Array of BMs 
    372         Array<ARX*> BMso; 
    373         //! pointer to Array of BMs 
    374         Array<ARX*> BMsp; 
    375         //!parameter evolution 
    376         shared_ptr<fnc> g; 
    377         //!observation function 
    378         shared_ptr<fnc> h; 
    379          
    380         public: 
    381                 void bayes(const vec &yt, const vec &cond ); 
    382                 void from_setting(const Setting &set) ; 
    383                 void validate() { 
    384                         bdm_assert(g->dimensionc()==g->dimension(),"not supported yet"); 
    385                         bdm_assert(h->dimensionc()==g->dimension(),"not supported yet");                         
    386                 } 
    387  
    388                 double logpred(const vec &cond) const NOT_IMPLEMENTED(0.0); 
    389                 epdf* epredictor() const NOT_IMPLEMENTED(NULL); 
    390                 pdf* predictor() const NOT_IMPLEMENTED(NULL); 
    391                 const epdf& posterior() const {return pf->posterior();}; 
    392                  
    393                 void log_register( logger &L, const string &prefix ){ 
    394                         BM::log_register(L,prefix); 
    395                         logrec->ids.set_size ( 3 ); 
    396                         logrec->ids(1)= L.add_vector(RV("Q",dimension()*dimension()), prefix+L.prefix_sep()+"Q"); 
    397                         logrec->ids(2)= L.add_vector(RV("R",dimensiony()*dimensiony()), prefix+L.prefix_sep()+"R"); 
    398                          
    399                 }; 
    400                 void log_write() const { 
    401                         BM::log_write(); 
    402                         mat mQ=zeros(dimension(),dimension()); 
    403                         mat pom=zeros(dimension(),dimension()); 
    404                         mat mR=zeros(dimensiony(),dimensiony()); 
    405                         mat pom2=zeros(dimensiony(),dimensiony()); 
    406                         mat dum; 
    407                         const vec w=pf->posterior()._w(); 
    408                         for (int i=0; i<w.length(); i++){ 
    409                                 BMsp(i)->posterior().mean_mat(dum,pom); 
    410                                 mQ += w(i) * pom; 
    411                                 BMso(i)->posterior().mean_mat(dum,pom2); 
    412                                 mR += w(i) * pom2; 
    413                                  
    414                         } 
    415                         logrec->L.log_vector ( logrec->ids ( 1 ), cvectorize(mQ) ); 
    416                         logrec->L.log_vector ( logrec->ids ( 2 ), cvectorize(mR) ); 
    417                          
    418                 } 
    419 }; 
    420 UIREGISTER(MPF_ARXg); 
     542// class MPF_ARXg :public BM{ 
     543//      protected: 
     544//      shared_ptr<PF> pf; 
     545//      //! pointer to Array of BMs 
     546//      Array<ARX*> BMso; 
     547//      //! pointer to Array of BMs 
     548//      Array<ARX*> BMsp; 
     549//      //!parameter evolution 
     550//      shared_ptr<fnc> g; 
     551//      //!observation function 
     552//      shared_ptr<fnc> h; 
     553//       
     554//      public: 
     555//              void bayes(const vec &yt, const vec &cond ); 
     556//              void from_setting(const Setting &set) ; 
     557//              void validate() { 
     558//                      bdm_assert(g->dimensionc()==g->dimension(),"not supported yet"); 
     559//                      bdm_assert(h->dimensionc()==g->dimension(),"not supported yet");                         
     560//              } 
     561//  
     562//              double logpred(const vec &cond) const NOT_IMPLEMENTED(0.0); 
     563//              epdf* epredictor() const NOT_IMPLEMENTED(NULL); 
     564//              pdf* predictor() const NOT_IMPLEMENTED(NULL); 
     565//              const epdf& posterior() const {return pf->posterior();}; 
     566//               
     567//              void log_register( logger &L, const string &prefix ){ 
     568//                      BM::log_register(L,prefix); 
     569//                      logrec->ids.set_size ( 3 ); 
     570//                      logrec->ids(1)= L.add_vector(RV("Q",dimension()*dimension()), prefix+L.prefix_sep()+"Q"); 
     571//                      logrec->ids(2)= L.add_vector(RV("R",dimensiony()*dimensiony()), prefix+L.prefix_sep()+"R"); 
     572//                       
     573//              }; 
     574//              void log_write() const { 
     575//                      BM::log_write(); 
     576//                      mat mQ=zeros(dimension(),dimension()); 
     577//                      mat pom=zeros(dimension(),dimension()); 
     578//                      mat mR=zeros(dimensiony(),dimensiony()); 
     579//                      mat pom2=zeros(dimensiony(),dimensiony()); 
     580//                      mat dum; 
     581//                      const vec w=pf->posterior()._w(); 
     582//                      for (int i=0; i<w.length(); i++){ 
     583//                              BMsp(i)->posterior().mean_mat(dum,pom); 
     584//                              mQ += w(i) * pom; 
     585//                              BMso(i)->posterior().mean_mat(dum,pom2); 
     586//                              mR += w(i) * pom2; 
     587//                               
     588//                      } 
     589//                      logrec->L.log_vector ( logrec->ids ( 1 ), cvectorize(mQ) ); 
     590//                      logrec->L.log_vector ( logrec->ids ( 2 ), cvectorize(mR) ); 
     591//                       
     592//              } 
     593// }; 
     594// UIREGISTER(MPF_ARXg); 
    421595 
    422596 
  • library/bdm/stat/exp_family.cpp

    r878 r887  
    423423} 
    424424 
    425 void eEmp::resample ( ivec &ind, RESAMPLING_METHOD method ) { 
     425void eEmp::resample ( RESAMPLING_METHOD method ) { 
     426        ivec ind = zeros_i ( n ); 
     427        bdm::resample(w,ind,method); 
     428        // copy the internals according to ind 
     429        for (int i = 0; i < n; i++ ) { 
     430                if ( ind ( i ) != i ) { 
     431                        samples ( i ) = samples ( ind ( i ) ); 
     432                } 
     433                w ( i ) = 1.0 / n; 
     434        } 
     435} 
     436 
     437void resample ( const vec &w, ivec &ind, RESAMPLING_METHOD method ) { 
     438        int n = w.length(); 
    426439        ind = zeros_i ( n ); 
    427440        ivec N_babies = zeros_i ( n ); 
     
    430443        int i, j, parent; 
    431444        double u0; 
    432  
     445         
    433446        switch ( method ) { 
    434         case MULTINOMIAL: 
    435                 u ( n - 1 ) = pow ( UniRNG.sample(), 1.0 / n ); 
    436  
    437                 for ( i = n - 2; i >= 0; i-- ) { 
    438                         u ( i ) = u ( i + 1 ) * pow ( UniRNG.sample(), 1.0 / ( i + 1 ) ); 
    439                 } 
    440  
    441                 break; 
    442  
    443         case STRATIFIED: 
    444  
    445                 for ( i = 0; i < n; i++ ) { 
    446                         u ( i ) = ( i + UniRNG.sample() ) / n; 
    447                 } 
    448  
    449                 break; 
    450  
    451         case SYSTEMATIC: 
    452                 u0 = UniRNG.sample(); 
    453  
    454                 for ( i = 0; i < n; i++ ) { 
    455                         u ( i ) = ( i + u0 ) / n; 
    456                 } 
    457  
    458                 break; 
    459  
    460         default: 
    461                 bdm_error ( "PF::resample(): Unknown resampling method" ); 
    462         } 
    463  
     447                case MULTINOMIAL: 
     448                        u ( n - 1 ) = pow ( UniRNG.sample(), 1.0 / n ); 
     449                         
     450                        for ( i = n - 2; i >= 0; i-- ) { 
     451                                u ( i ) = u ( i + 1 ) * pow ( UniRNG.sample(), 1.0 / ( i + 1 ) ); 
     452                        } 
     453                         
     454                        break; 
     455                         
     456                case STRATIFIED: 
     457                         
     458                        for ( i = 0; i < n; i++ ) { 
     459                                u ( i ) = ( i + UniRNG.sample() ) / n; 
     460                        } 
     461                         
     462                        break; 
     463                         
     464                case SYSTEMATIC: 
     465                        u0 = UniRNG.sample(); 
     466                         
     467                        for ( i = 0; i < n; i++ ) { 
     468                                u ( i ) = ( i + u0 ) / n; 
     469                        } 
     470                         
     471                        break; 
     472                         
     473                default: 
     474                        bdm_error ( "PF::resample(): Unknown resampling method" ); 
     475        } 
     476         
    464477        // U is now full 
    465478        j = 0; 
    466  
     479         
    467480        for ( i = 0; i < n; i++ ) { 
    468481                while ( u ( i ) > cumDist ( j ) ) j++; 
    469  
     482                 
    470483                N_babies ( j ) ++; 
    471484        } 
     
    474487        // * particles with at least one baby should not move * 
    475488        // This assures that reassignment can be done inplace; 
    476  
     489         
    477490        // find the first parent; 
    478491        parent = 0; 
    479492        while ( N_babies ( parent ) == 0 ) parent++; 
    480  
     493         
    481494        // Build index 
    482495        for ( i = 0; i < n; i++ ) { 
     
    488501                        // if yes, find the next one 
    489502                        while ( ( N_babies ( parent ) == 0 ) || ( N_babies ( parent ) == 1 && parent > i ) ) parent++; 
    490  
     503                         
    491504                        // Replicate parent 
    492505                        ind ( i ) = parent; 
    493  
     506                         
    494507                        N_babies ( parent ) --; //this index was now replicated; 
    495508                } 
    496  
    497         } 
    498  
    499 // copy the internals according to ind 
    500         for ( i = 0; i < n; i++ ) { 
    501                 if ( ind ( i ) != i ) { 
    502                         samples ( i ) = samples ( ind ( i ) ); 
    503                 } 
    504                 w ( i ) = 1.0 / n; 
     509                 
    505510        } 
    506511} 
  • library/bdm/stat/exp_family.h

    r878 r887  
    9797        virtual void flatten ( const BMEF * B ) NOT_IMPLEMENTED_VOID; 
    9898 
    99         double logpred ( const vec &yt ) const NOT_IMPLEMENTED(0); 
    100          
    101         virtual epdf* epredictor() const NOT_IMPLEMENTED(NULL); 
    102  
    103         virtual pdf* predictor() const NOT_IMPLEMENTED(NULL); 
    10499 
    105100        void to_setting ( Setting &set ) const 
     
    267262UIREGISTER2 ( enorm, fsqmat ); 
    268263SHAREDPTR2 ( enorm, fsqmat ); 
     264 
     265typedef enorm<ldmat> egauss; 
     266UIREGISTER(egauss); 
     267 
    269268 
    270269//forward declaration 
     
    15961595//! Switch between various resampling methods. 
    15971596enum RESAMPLING_METHOD { MULTINOMIAL = 0, STRATIFIED = 1, SYSTEMATIC = 3 }; 
     1597 
     1598//! Shortcut for multinomial.sample(int n). Various simplifications are allowed see RESAMPLING_METHOD  
     1599void resample(const vec &w, ivec &ind, RESAMPLING_METHOD=SYSTEMATIC); 
     1600 
    15981601/*! 
    15991602\brief Weighted empirical density 
     
    16641667        }; 
    16651668        //! Function performs resampling, i.e. removal of low-weight samples and duplication of high-weight samples such that the new samples represent the same density. 
    1666         //! The vector with indeces of new samples is returned in variable \c index. 
    1667         void resample ( ivec &index, RESAMPLING_METHOD method = SYSTEMATIC ); 
    1668  
    1669         //! Resampling without returning index of new particles. 
    1670         void resample ( RESAMPLING_METHOD method = SYSTEMATIC ) { 
    1671                 ivec ind; 
    1672                 resample ( ind, method ); 
    1673         }; 
     1669        void resample ( RESAMPLING_METHOD method = SYSTEMATIC ); 
    16741670 
    16751671        //! inherited operation : NOT implemented 
     
    18741870} 
    18751871 
     1872/*! Dirac delta function distribution */ 
     1873class dirac: public epdf{ 
     1874        protected:  
     1875                vec point; 
     1876        public: 
     1877                double evallog (const vec &dt) const {return -inf;} 
     1878                vec mean () const {return point;} 
     1879                vec variance () const {return pow(point,2);} 
     1880                void qbounds ( vec &lb, vec &ub, double percentage = 0.95 ) const { lb = point; ub = point;} 
     1881                //! access 
     1882                const vec& _point() {return point;} 
     1883                void set_point(const vec& p){point =p; dim=p.length();} 
     1884                vec sample() const {return point;} 
     1885        }; 
     1886 
    18761887//// 
    18771888///////