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

new base for particle filtering

Files:
1 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