particles.h
Go to the documentation of this file.00001 00013 #ifndef PARTICLES_H 00014 #define PARTICLES_H 00015 00016 00017 #include "../estim/arx_ext.h" 00018 #include "../stat/emix.h" 00019 00020 namespace bdm { 00021 00023 class MarginalizedParticleBase : public BM { 00024 protected: 00026 dirac est_emp; 00028 shared_ptr<BM> bm; 00029 00031 class eprod_2:public eprod_base { 00032 protected: 00033 MarginalizedParticleBase ∓ 00034 public: 00035 eprod_2(MarginalizedParticleBase &m):mp(m) {} 00036 const epdf* factor(int i) const { 00037 return (i==0) ? &mp.bm->posterior() : &mp.est_emp; 00038 } 00039 const int no_factors() const { 00040 return 2; 00041 } 00042 } est; 00043 00044 public: 00045 MarginalizedParticleBase():est(*this) {}; 00046 MarginalizedParticleBase(const MarginalizedParticleBase &m2):BM(m2),est(*this) { 00047 bm = m2.bm->_copy(); 00048 est_emp = m2.est_emp; 00049 est.validate(); 00050 validate(); 00051 }; 00052 void bayes(const vec &dt, const vec &cond) NOT_IMPLEMENTED_VOID; 00053 00054 const eprod_2& posterior() const { 00055 return est; 00056 } 00057 00058 void set_prior(const epdf *pdf0) { 00059 const eprod *ep=dynamic_cast<const eprod*>(pdf0); 00060 if (ep) { // full prior 00061 bdm_assert(ep->no_factors()==2,"Incompatible prod"); 00062 bm->set_prior(ep->factor(0)); 00063 est_emp.set_point(ep->factor(1)->sample()); 00064 } else { 00065 // assume prior is only for emp; 00066 est_emp.set_point(pdf0->sample()); 00067 } 00068 } 00069 00070 00080 void from_setting(const Setting &set) { 00081 BM::from_setting ( set ); 00082 bm = UI::build<BM> ( set, "bm", UI::compulsory ); 00083 } 00084 void validate() { 00085 BM::validate(); 00086 //est.validate(); --pdfs not known 00087 bdm_assert(bm,"Internal BM is not given"); 00088 } 00089 }; 00090 00092 class MarginalizedParticle : public MarginalizedParticleBase { 00093 protected: 00095 shared_ptr<pdf> par; // pdf for non-linear part 00097 shared_ptr<datalink_part> cond2bm; 00099 shared_ptr<datalink_part> cond2par; 00101 shared_ptr<datalink_part> emp2bm; 00103 shared_ptr<datalink_part> emp2par; 00104 00105 public: 00106 BM* _copy() const { 00107 return new MarginalizedParticle(*this); 00108 }; 00109 void bayes(const vec &dt, const vec &cond) { 00110 vec par_cond(par->dimensionc()); 00111 cond2par->filldown(cond,par_cond); // copy ut 00112 emp2par->filldown(est_emp._point(),par_cond); // copy xt-1 00113 00114 //sample new particle 00115 est_emp.set_point(par->samplecond(par_cond)); 00116 //if (evalll) 00117 vec bm_cond(bm->dimensionc()); 00118 cond2bm->filldown(cond, bm_cond);// set e.g. ut 00119 emp2bm->filldown(est_emp._point(), bm_cond);// set e.g. ut 00120 bm->bayes(dt, bm_cond); 00121 ll=bm->_ll(); 00122 } 00123 00133 void from_setting(const Setting &set) { 00134 MarginalizedParticleBase::from_setting ( set ); 00135 par = UI::build<pdf> ( set, "parameter_pdf", UI::compulsory ); 00136 } 00137 00138 void to_setting(Setting &set)const { 00139 MarginalizedParticleBase::to_setting(set); 00140 UI::save(par,set,"parameter_pdf"); 00141 UI::save(bm,set,"bm"); 00142 } 00143 void validate() { 00144 MarginalizedParticleBase::validate(); 00145 est_emp.set_rv(par->_rv()); 00146 if (est_emp.point.length()!=par->dimension()) 00147 est_emp.set_point(zeros(par->dimension())); 00148 est.validate(); 00149 00150 yrv = bm->_yrv(); 00151 dimy = bm->dimensiony(); 00152 set_rv( concat(bm->_rv(), par->_rv())); 00153 set_dim( par->dimension()+bm->dimension()); 00154 00155 rvc = par->_rvc(); 00156 rvc.add(bm->_rvc()); 00157 rvc=rvc.subt(par->_rv()); 00158 rvc=rvc.subt(par->_rv().copy_t(-1)); 00159 rvc=rvc.subt(bm->_rv().copy_t(-1)); // 00160 00161 cond2bm=new datalink_part; 00162 cond2par=new datalink_part; 00163 emp2bm =new datalink_part; 00164 emp2par =new datalink_part; 00165 cond2bm->set_connection(bm->_rvc(), rvc); 00166 cond2par->set_connection(par->_rvc(), rvc); 00167 emp2bm->set_connection(bm->_rvc(), par->_rv()); 00168 emp2par->set_connection(par->_rvc(), par->_rv().copy_t(-1)); 00169 00170 dimc = rvc._dsize(); 00171 }; 00172 }; 00173 UIREGISTER(MarginalizedParticle); 00174 00176 class BootstrapParticle : public BM { 00177 dirac est; 00178 shared_ptr<pdf> par; 00179 shared_ptr<pdf> obs; 00180 shared_ptr<datalink_part> cond2par; 00181 shared_ptr<datalink_part> cond2obs; 00182 shared_ptr<datalink_part> xt2obs; 00183 shared_ptr<datalink_part> xtm2par; 00184 public: 00185 BM* _copy() const { 00186 return new BootstrapParticle(*this); 00187 }; 00188 void bayes(const vec &dt, const vec &cond) { 00189 vec par_cond(par->dimensionc()); 00190 cond2par->filldown(cond,par_cond); // copy ut 00191 xtm2par->filldown(est._point(),par_cond); // copy xt-1 00192 00193 //sample new particle 00194 est.set_point(par->samplecond(par_cond)); 00195 //if (evalll) 00196 vec obs_cond(obs->dimensionc()); 00197 cond2obs->filldown(cond, obs_cond);// set e.g. ut 00198 xt2obs->filldown(est._point(), obs_cond);// set e.g. ut 00199 ll=obs->evallogcond(dt,obs_cond); 00200 } 00201 const dirac& posterior() const { 00202 return est; 00203 } 00204 00205 void set_prior(const epdf *pdf0) { 00206 est.set_point(pdf0->sample()); 00207 } 00208 00218 void from_setting(const Setting &set) { 00219 BM::from_setting ( set ); 00220 par = UI::build<pdf> ( set, "parameter_pdf", UI::compulsory ); 00221 obs = UI::build<pdf> ( set, "observation_pdf", UI::compulsory ); 00222 } 00223 00224 void validate() { 00225 yrv = obs->_rv(); 00226 dimy = obs->dimension(); 00227 set_rv( par->_rv()); 00228 set_dim( par->dimension()); 00229 00230 rvc = par->_rvc().subt(par->_rv().copy_t(-1)); 00231 rvc.add(obs->_rvc().subt(par->_rv())); // 00232 00233 cond2obs=new datalink_part; 00234 cond2par=new datalink_part; 00235 xt2obs =new datalink_part; 00236 xtm2par =new datalink_part; 00237 cond2obs->set_connection(obs->_rvc(), rvc); 00238 cond2par->set_connection(par->_rvc(), rvc); 00239 xt2obs->set_connection(obs->_rvc(), _rv()); 00240 xtm2par->set_connection(par->_rvc(), _rv().copy_t(-1)); 00241 00242 dimc = rvc._dsize(); 00243 }; 00244 }; 00245 UIREGISTER(BootstrapParticle); 00246 00247 00254 class PF : public BM { 00257 00260 LOG_LEVEL(PF,logweights,logmeans,logvars,logNeff); 00261 00262 class pf_mix: public emix_base { 00263 Array<BM*> &bms; 00264 public: 00265 pf_mix(vec &w0, Array<BM*> &bms0):emix_base(w0),bms(bms0) {} 00266 const epdf* component(const int &i)const { 00267 return &(bms(i)->posterior()); 00268 } 00269 int no_coms() const { 00270 return bms.length(); 00271 } 00272 }; 00273 protected: 00275 int n; 00277 pf_mix est; 00279 vec w; 00281 Array<BM*> particles; 00283 vec lls; 00284 00286 RESAMPLING_METHOD resmethod; 00289 double res_threshold; 00290 00294 00295 public: 00298 PF ( ) : est(w,particles) { }; 00299 00300 void set_parameters ( int n0, double res_th0 = 0.5, RESAMPLING_METHOD rm = SYSTEMATIC ) { 00301 n = n0; 00302 res_threshold = res_th0; 00303 resmethod = rm; 00304 }; 00305 void set_model ( const BM *particle0, const epdf *prior) { 00306 if (n>0) { 00307 particles.set_length(n); 00308 for (int i=0; i<n; i++) { 00309 particles(i) = particle0->_copy(); 00310 particles(i)->set_prior(prior); 00311 } 00312 } 00313 // set values for posterior 00314 est.set_rv ( particle0->posterior()._rv() ); 00315 }; 00316 void set_statistics ( const vec w0, const epdf &epdf0 ) { 00317 //est.set_statistics ( w0, epdf0 ); 00318 }; 00319 /* void set_statistics ( const eEmp &epdf0 ) { 00320 bdm_assert_debug ( epdf0._rv().equal ( par->_rv() ), "Incompatible input" ); 00321 est = epdf0; 00322 };*/ 00324 00326 virtual void bayes_weights(); 00328 virtual bool do_resampling() { 00329 double eff = 1.0 / ( w * w ); 00330 return eff < ( res_threshold*n ); 00331 } 00332 void bayes ( const vec &yt, const vec &cond ); 00334 vec& _lls() { 00335 return lls; 00336 } 00338 RESAMPLING_METHOD _resmethod() const { 00339 return resmethod; 00340 } 00342 const pf_mix& posterior() const { 00343 return est; 00344 } 00345 00358 void from_setting ( const Setting &set ) { 00359 BM::from_setting ( set ); 00360 UI::get ( log_level, set, "log_level", UI::optional ); 00361 00362 shared_ptr<BM> bm0 = UI::build<BM>(set, "particle",UI::compulsory); 00363 00364 n =0; 00365 UI::get(n,set,"n",UI::optional);; 00366 if (n>0) { 00367 particles.set_length(n); 00368 for(int i=0; i<n; i++) { 00369 particles(i)=bm0->_copy(); 00370 } 00371 w = ones(n)/n; 00372 } 00373 shared_ptr<epdf> pri = UI::build<epdf>(set,"prior"); 00374 set_prior(pri.get()); 00375 // set resampling method 00376 resmethod_from_set ( set ); 00377 //set drv 00378 00379 rvc = bm0->_rvc(); 00380 dimc = bm0->dimensionc(); 00381 BM::set_rv(bm0->_rv()); 00382 yrv=bm0->_yrv(); 00383 dimy = bm0->dimensiony(); 00384 } 00385 00386 void to_setting(Setting &set) const{ 00387 BM::to_setting(set); 00388 UI::save(particles, set,"particles"); 00389 UI::save(w,set,"w"); 00390 } 00391 00392 void log_register ( bdm::logger& L, const string& prefix ) { 00393 BM::log_register(L,prefix); 00394 if (log_level[logweights]) { 00395 L.add_vector( log_level, logweights, RV ( particles.length()), prefix); 00396 } 00397 if (log_level[logNeff]) { 00398 L.add_vector( log_level, logNeff, RV ( 1), prefix); 00399 } 00400 if (log_level[logmeans]) { 00401 for (int i=0; i<particles.length(); i++) { 00402 L.add_vector( log_level, logmeans, RV ( particles(i)->dimension() ), prefix , i); 00403 } 00404 } 00405 if (log_level[logvars]) { 00406 for (int i=0; i<particles.length(); i++) { 00407 L.add_vector( log_level, logvars, RV ( particles(i)->dimension() ), prefix , i); 00408 } 00409 } 00410 }; 00411 void log_write ( ) const { 00412 BM::log_write(); 00413 // weigths are before resamplign -- bayes 00414 if (log_level[logmeans]) { 00415 for (int i=0; i<particles.length(); i++) { 00416 log_level.store( logmeans, particles(i)->posterior().mean(), i); 00417 } 00418 } 00419 if (log_level[logvars]) { 00420 for (int i=0; i<particles.length(); i++) { 00421 log_level.store( logvars, particles(i)->posterior().variance(), i); 00422 } 00423 } 00424 00425 } 00426 00427 void set_prior(const epdf *pri) { 00428 const emix_base *emi=dynamic_cast<const emix_base*>(pri); 00429 if (emi) { 00430 bdm_assert(particles.length()>0, "initial particle is not assigned"); 00431 n = emi->_w().length(); 00432 int old_n = particles.length(); 00433 if (n!=old_n) { 00434 particles.set_length(n,true); 00435 } 00436 for(int i=old_n; i<n; i++) { 00437 particles(i)=particles(0)->_copy(); 00438 } 00439 00440 for (int i =0; i<n; i++) { 00441 particles(i)->set_prior(emi->_com(i)); 00442 } 00443 } else { 00444 // try to find "n" 00445 bdm_assert(n>0, "Field 'n' must be filled when prior is not of type emix"); 00446 for (int i =0; i<n; i++) { 00447 particles(i)->set_prior(pri); 00448 } 00449 00450 } 00451 } 00453 void resmethod_from_set ( const Setting &set ) { 00454 string resmeth; 00455 if ( UI::get ( resmeth, set, "resmethod", UI::optional ) ) { 00456 if ( resmeth == "systematic" ) { 00457 resmethod = SYSTEMATIC; 00458 } else { 00459 if ( resmeth == "multinomial" ) { 00460 resmethod = MULTINOMIAL; 00461 } else { 00462 if ( resmeth == "stratified" ) { 00463 resmethod = STRATIFIED; 00464 } else { 00465 bdm_error ( "Unknown resampling method" ); 00466 } 00467 } 00468 } 00469 } else { 00470 resmethod = SYSTEMATIC; 00471 }; 00472 if ( !UI::get ( res_threshold, set, "res_threshold", UI::optional ) ) { 00473 res_threshold = 0.9; 00474 } 00475 //validate(); 00476 } 00477 00478 void validate() { 00479 BM::validate(); 00480 est.validate(); 00481 bdm_assert ( n>0, "empty particle pool" ); 00482 n = w.length(); 00483 lls = zeros ( n ); 00484 00485 if ( particles(0)->_rv()._dsize() > 0 ) { 00486 bdm_assert ( particles(0)->_rv()._dsize() == est.dimension(), "MPF:: Mismatch of RV " +particles(0)->_rv().to_string() + 00487 " of size (" +num2str(particles(0)->_rv()._dsize())+") and dimension of posterior ("+num2str(est.dimension()) + ")" ); 00488 } 00489 } 00491 void resample ( ) { 00492 ivec ind = zeros_i ( n ); 00493 bdm::resample(w,ind,resmethod); 00494 // copy the internals according to ind 00495 for (int i = 0; i < n; i++ ) { 00496 if ( ind ( i ) != i ) { 00497 delete particles(i); 00498 particles( i ) = particles( ind ( i ) )->_copy(); 00499 } 00500 w ( i ) = 1.0 / n; 00501 } 00502 } 00504 Array<BM*>& _particles() { 00505 return particles; 00506 } 00507 ~PF() { 00508 for (int i=0; i<particles.length(); i++) { 00509 delete particles(i); 00510 } 00511 } 00512 00513 }; 00514 UIREGISTER ( PF ); 00515 00526 class NoiseParticleX : public MarginalizedParticleBase { 00527 protected: 00529 shared_ptr<fnc> g; // pdf for non-linear part 00531 shared_ptr<pdf> fy; // pdf for non-linear part 00532 00533 RV rvx; 00534 RV rvxc; 00535 RV rvyc; 00536 00538 datalink_part cond2g; 00540 datalink_part cond2fy; 00542 datalink_part x2g; 00544 datalink_part x2fy; 00545 00546 public: 00547 BM* _copy() const { 00548 return new NoiseParticleX(*this); 00549 }; 00550 void bayes(const vec &dt, const vec &cond) { 00551 //shared_ptr<epdf> pred_v=bm->epredictor(); 00552 00553 vec xt; 00554 vec xthat; 00555 //vec vt=pred_v->sample(); 00556 // vec vt = bm->samplepred(); 00557 if (1){ 00558 //vec vt=pred_v->sample(); 00559 int dim_x = bm->dimensiony(); 00560 int dim_y = dimensiony(); 00561 00562 vec post_mean =bm->posterior().mean(); 00563 mat SigV; 00564 if (post_mean.length()==dim_x){ 00565 SigV=diag(vec(post_mean._data(), dim_x)); 00566 } else { 00567 SigV=mat(post_mean._data(), dim_x, dim_x); 00568 } 00569 00570 vec &xtm=est_emp.point; 00571 vec g_args(g->dimensionc()); 00572 x2g.filldown(xtm,g_args); 00573 cond2g.filldown(cond,g_args); 00574 xthat = g->eval(g_args); 00575 00576 mat IC = eye(dim_x+dim_y); 00577 IC.set_submatrix(0,0,eye(dim_x)); 00578 mat SigJo=eye(dim_x + dim_y); 00579 SigJo.set_submatrix(0,0,SigV); 00580 mlnorm<ldmat>* mlfy=dynamic_cast<mlnorm<ldmat>* >(fy.get()); 00581 00582 SigJo.set_submatrix(dim_x,dim_x, mlfy->_R()); 00583 IC.set_submatrix(dim_x,0,-mlfy->_A()); 00584 00585 mat iIC=inv(IC); 00586 enorm<chmat> Jo; 00587 Jo._mu()=iIC*concat(xthat, zeros(dim_y)); 00588 Jo._R()=iIC*SigJo*iIC.T(); 00589 Jo.set_rv(concat(bm->_yrv(),yrv)); 00590 Jo.validate(); 00591 shared_ptr<pdf> Cond=Jo.condition(bm->_yrv()); 00592 //new sample 00593 00594 xt = Cond->samplecond(dt);// + vt; 00595 } else{ 00596 00597 00598 //new sample 00599 vec &xtm=est_emp.point; 00600 vec g_args(g->dimensionc()); 00601 x2g.filldown(xtm,g_args); 00602 cond2g.filldown(cond,g_args); 00603 xthat = g->eval(g_args); 00604 00605 xt = xthat + bm->samplepred(); 00606 00607 } 00608 est_emp.point=xt; 00609 00610 // the vector [v_t] updates bm, 00611 bm->bayes(xt-xthat); 00612 00613 // residue of observation 00614 vec fy_args(fy->dimensionc()); 00615 x2fy.filldown(xt,fy_args); 00616 cond2fy.filldown(cond,fy_args); 00617 00618 ll= fy->evallogcond(dt,fy_args); 00619 } 00620 void from_setting(const Setting &set) { 00621 MarginalizedParticleBase::from_setting(set); //reads bm, yrv,rvc, bm_rv, etc... 00622 00623 g=UI::build<fnc>(set,"g",UI::compulsory); 00624 fy=UI::build<pdf>(set,"fy",UI::compulsory); 00625 UI::get(rvx,set,"rvx",UI::compulsory); 00626 est_emp.set_rv(rvx); 00627 00628 UI::get(rvxc,set,"rvxc",UI::compulsory); 00629 UI::get(rvyc,set,"rvyc",UI::compulsory); 00630 00631 } 00632 void to_setting (Setting &set) const { 00633 MarginalizedParticleBase::to_setting(set); //reads bm, yrv,rvc, bm_rv, etc... 00634 UI::save(g,set,"g"); 00635 UI::save(fy,set,"fy"); 00636 UI::save(bm,set,"bm"); 00637 } 00638 void validate() { 00639 MarginalizedParticleBase::validate(); 00640 00641 dimy = fy->dimension(); 00642 bm->set_yrv(rvx); 00643 00644 est_emp.set_rv(rvx); 00645 est_emp.set_dim(rvx._dsize()); 00646 est.validate(); 00647 // 00648 //check dimensions 00649 rvc = rvxc.subt(rvx.copy_t(-1)); 00650 rvc.add( rvyc); 00651 rvc=rvc.subt(rvx); 00652 00653 bdm_assert(g->dimension()==rvx._dsize(),"rvx is not described"); 00654 bdm_assert(g->dimensionc()==rvxc._dsize(),"rvxc is not described"); 00655 bdm_assert(fy->dimensionc()==rvyc._dsize(),"rvyc is not described"); 00656 00657 bdm_assert(bm->dimensiony()==g->dimension(), 00658 "Incompatible noise estimator of dimension " + 00659 num2str(bm->dimensiony()) + " does not match dimension of g , " + 00660 num2str(g->dimension())); 00661 00662 dimc = rvc._dsize(); 00663 00664 //establish datalinks 00665 x2g.set_connection(rvxc, rvx.copy_t(-1)); 00666 cond2g.set_connection(rvxc, rvc); 00667 00668 x2fy.set_connection(rvyc, rvx); 00669 cond2fy.set_connection(rvyc, rvc); 00670 } 00671 }; 00672 UIREGISTER(NoiseParticleX); 00673 00674 00685 class NoiseParticleXY : public BM { 00686 protected: 00688 dirac est_emp; 00690 shared_ptr<BM> bmx; 00691 shared_ptr<BM> bmy; 00692 00694 class eprod_3:public eprod_base { 00695 protected: 00696 NoiseParticleXY ∓ 00697 public: 00698 eprod_3(NoiseParticleXY &m):mp(m) {} 00699 const epdf* factor(int i) const { 00700 if (i==0) return &mp.bmx->posterior() ; 00701 if(i==1) return &mp.bmy->posterior(); 00702 return &mp.est_emp; 00703 } 00704 const int no_factors() const { 00705 return 3; 00706 } 00707 } est; 00708 00709 protected: 00711 shared_ptr<fnc> g; // pdf for non-linear part 00713 shared_ptr<fnc> h; // pdf for non-linear part 00714 00715 RV rvx; 00716 RV rvxc; 00717 RV rvyc; 00718 00720 datalink_part cond2g; 00722 datalink_part cond2h; 00724 datalink_part x2g; 00726 datalink_part x2h; 00727 00728 public: 00729 NoiseParticleXY():est(*this) {}; 00730 NoiseParticleXY(const NoiseParticleXY &m2):BM(m2),est(*this),g(m2.g),h(m2.h), rvx(m2.rvx),rvxc(m2.rvxc),rvyc(m2.rvyc) { 00731 bmx = m2.bmx->_copy(); 00732 bmy = m2.bmy->_copy(); 00733 est_emp = m2.est_emp; 00734 //est.validate(); 00735 validate(); 00736 }; 00737 00738 const eprod_3& posterior() const { 00739 return est; 00740 } 00741 00742 void set_prior(const epdf *pdf0) { 00743 const eprod *ep=dynamic_cast<const eprod*>(pdf0); 00744 if (ep) { // full prior 00745 bdm_assert(ep->no_factors()==2,"Incompatible prod"); 00746 bmx->set_prior(ep->factor(0)); 00747 bmy->set_prior(ep->factor(1)); 00748 est_emp.set_point(ep->factor(2)->sample()); 00749 } else { 00750 // assume prior is only for emp; 00751 est_emp.set_point(pdf0->sample()); 00752 } 00753 } 00754 00755 BM* _copy() const { 00756 return new NoiseParticleXY(*this); 00757 }; 00758 00759 void bayes(const vec &dt, const vec &cond) { 00760 //shared_ptr<epdf> pred_v=bm->epredictor(); 00761 00762 //vec vt=pred_v->sample(); 00763 vec vt = bmx->samplepred(); 00764 00765 //new sample 00766 vec &xtm=est_emp.point; 00767 vec g_args(g->dimensionc()); 00768 x2g.filldown(xtm,g_args); 00769 cond2g.filldown(cond,g_args); 00770 vec xt = g->eval(g_args) + vt; 00771 est_emp.point=xt; 00772 00773 // the vector [v_t] updates bm, 00774 bmx->bayes(vt); 00775 00776 // residue of observation 00777 vec h_args(h->dimensionc()); 00778 x2h.filldown(xt,h_args); 00779 cond2h.filldown(cond,h_args); 00780 00781 vec z_y =h->eval(h_args)-dt; 00782 // ARX *abm = dynamic_cast<ARX*>(bmy.get()); 00783 /* double ll2; 00784 if (abm){ //ARX 00785 shared_ptr<epdf> pr_y(abm->epredictor_student(empty_vec)); 00786 ll2=pr_y->evallog(z_y); 00787 } else{ 00788 shared_ptr<epdf> pr_y(bmy->epredictor(empty_vec)); 00789 ll2=pr_y->evallog(z_y); 00790 }*/ 00791 00792 bmy->bayes(z_y); 00793 // test _lls 00794 ll= bmy->_ll(); 00795 } 00796 void from_setting(const Setting &set) { 00797 BM::from_setting(set); //reads bm, yrv,rvc, bm_rv, etc... 00798 bmx = UI::build<BM>(set,"bmx",UI::compulsory); 00799 00800 bmy = UI::build<BM>(set,"bmy",UI::compulsory); 00801 g=UI::build<fnc>(set,"g",UI::compulsory); 00802 h=UI::build<fnc>(set,"h",UI::compulsory); 00803 UI::get(rvx,set,"rvx",UI::compulsory); 00804 est_emp.set_rv(rvx); 00805 00806 UI::get(rvxc,set,"rvxc",UI::compulsory); 00807 UI::get(rvyc,set,"rvyc",UI::compulsory); 00808 00809 } 00810 void to_setting (Setting &set) const { 00811 BM::to_setting(set); //reads bm, yrv,rvc, bm_rv, etc... 00812 UI::save(g,set,"g"); 00813 UI::save(h,set,"h"); 00814 UI::save(bmx,set,"bmx"); 00815 UI::save(bmy,set,"bmy"); 00816 } 00817 void validate() { 00818 BM::validate(); 00819 00820 dimy = h->dimension(); 00821 // bmx->set_yrv(rvx); 00822 // bmy-> 00823 00824 est_emp.set_rv(rvx); 00825 est_emp.set_dim(rvx._dsize()); 00826 est.validate(); 00827 // 00828 //check dimensions 00829 rvc = rvxc.subt(rvx.copy_t(-1)); 00830 rvc.add( rvyc); 00831 rvc=rvc.subt(rvx); 00832 00833 bdm_assert(g->dimension()==rvx._dsize(),"rvx is not described"); 00834 bdm_assert(g->dimensionc()==rvxc._dsize(),"rvxc is not described"); 00835 bdm_assert(h->dimensionc()==rvyc._dsize(),"rvyc is not described"); 00836 00837 bdm_assert(bmx->dimensiony()==g->dimension(), 00838 "Incompatible noise estimator of dimension " + 00839 num2str(bmx->dimensiony()) + " does not match dimension of g , " + 00840 num2str(g->dimension()) 00841 ); 00842 00843 dimc = rvc._dsize(); 00844 00845 //establish datalinks 00846 x2g.set_connection(rvxc, rvx.copy_t(-1)); 00847 cond2g.set_connection(rvxc, rvc); 00848 00849 x2h.set_connection(rvyc, rvx); 00850 cond2h.set_connection(rvyc, rvc); 00851 } 00852 00853 }; 00854 UIREGISTER(NoiseParticleXY); 00855 00856 class NoiseParticleXYprop: public NoiseParticleXY{ 00857 public: 00858 BM* _copy() const { 00859 return new NoiseParticleXYprop(*this); 00860 }; 00861 void bayes(const vec &dt, const vec &cond) { 00862 int dim_x = bmx->dimensiony(); 00863 int dim_y = bmy->dimensiony(); 00864 00865 //vec vt=pred_v->sample(); 00866 mat SigV=mat(bmx->posterior().mean()._data(), dim_x, dim_x); 00867 mat SigW=mat(bmy->posterior().mean()._data(), dim_y, dim_y); 00868 00869 vec &xtm=est_emp.point; 00870 vec g_args(g->dimensionc()); 00871 x2g.filldown(xtm,g_args); 00872 cond2g.filldown(cond,g_args); 00873 vec xthat = g->eval(g_args); 00874 00875 mat IC = eye(dim_x+dim_y); 00876 IC.set_submatrix(0,0,eye(dim_x)); 00877 mat SigJo=eye(dim_x + dim_y); 00878 SigJo.set_submatrix(0,0,SigV); 00879 SigJo.set_submatrix(dim_x,dim_x, 10*SigW); 00880 00881 diffbifn* hb=dynamic_cast<bilinfn*>(h.get()); 00882 mat C=zeros(dim_y,dim_x); 00883 if (hb){ 00884 hb->dfdx_cond(xtm, cond, C,true); 00885 IC.set_submatrix(dim_x,0,-C); 00886 } 00887 00888 mat iIC=inv(IC); 00889 enorm<chmat> Jo; 00890 Jo._mu()=iIC*concat(xthat, zeros(dim_y)); 00891 Jo._R()=iIC*SigJo*iIC.T(); 00892 Jo.set_rv(concat(bmx->_yrv(), bmy->_yrv())); 00893 Jo.validate(); 00894 shared_ptr<pdf> Cond=Jo.condition(bmx->_yrv()); 00895 //new sample 00896 00897 vec xt = Cond->samplecond(dt);// + vt; 00898 00899 est_emp.point=xt; 00900 00901 // the vector [v_t] updates bm, 00902 bmx->bayes(xt-xthat); 00903 00904 // residue of observation 00905 vec h_args(h->dimensionc()); 00906 x2h.filldown(xt,h_args); 00907 cond2h.filldown(cond,h_args); 00908 00909 vec z_y =h->eval(h_args)-dt; 00910 // ARX *abm = dynamic_cast<ARX*>(bmy.get()); 00911 /* double ll2; 00912 i f (abm){ //ARX* 00913 shared_ptr<epdf> pr_y(abm->epredictor_student(empty_vec)); 00914 ll2=pr_y->evallog(z_y); 00915 } else{ 00916 shared_ptr<epdf> pr_y(bmy->epredictor(empty_vec)); 00917 ll2=pr_y->evallog(z_y); 00918 }*/ 00919 00920 bmy->bayes(z_y); 00921 // test _lls 00922 ll= bmy->_ll(); 00923 } 00924 }; 00925 UIREGISTER(NoiseParticleXYprop); 00926 00937 class NoiseParticle : public MarginalizedParticleBase { 00938 protected: 00940 shared_ptr<fnc> g; // pdf for non-linear part 00942 shared_ptr<fnc> h; // pdf for non-linear part 00943 00944 RV rvx; 00945 RV rvxc; 00946 RV rvyc; 00947 00949 datalink_part cond2g; 00951 datalink_part cond2h; 00953 datalink_part x2g; 00955 datalink_part x2h; 00956 00957 public: 00958 BM* _copy() const { 00959 return new NoiseParticle(*this); 00960 }; 00961 void bayes(const vec &dt, const vec &cond) { 00962 shared_ptr<epdf> pred_vw=bm->epredictor(); 00963 shared_ptr<epdf> pred_v = pred_vw->marginal(rvx); 00964 00965 vec vt=pred_v->sample(); 00966 00967 //new sample 00968 vec &xtm=est_emp.point; 00969 vec g_args(g->dimensionc()); 00970 x2g.filldown(xtm,g_args); 00971 cond2g.filldown(cond,g_args); 00972 vec xt = g->eval(g_args) + vt; 00973 est_emp.point=xt; 00974 00975 // residue of observation 00976 vec h_args(h->dimensionc()); 00977 x2h.filldown(xt,h_args); 00978 cond2h.filldown(cond,h_args); 00979 vec wt = dt-h->eval(h_args); 00980 // the vector [v_t,w_t] is now complete 00981 bm->bayes(concat(vt,wt)); 00982 ll=bm->_ll(); 00983 } 00984 void from_setting(const Setting &set) { 00985 MarginalizedParticleBase::from_setting(set); //reads bm, yrv,rvc, bm_rv, etc... 00986 00987 UI::get(g,set,"g",UI::compulsory); 00988 UI::get(h,set,"h",UI::compulsory); 00989 UI::get(rvx,set,"rvx",UI::compulsory); 00990 est_emp.set_rv(rvx); 00991 00992 UI::get(rvxc,set,"rvxc",UI::compulsory); 00993 UI::get(rvyc,set,"rvyc",UI::compulsory); 00994 00995 } 00996 void validate() { 00997 MarginalizedParticleBase::validate(); 00998 00999 dimy = h->dimension(); 01000 bm->set_yrv(concat(rvx,yrv)); 01001 01002 est_emp.set_rv(rvx); 01003 est_emp.set_dim(rvx._dsize()); 01004 est.validate(); 01005 // 01006 //check dimensions 01007 rvc = rvxc.subt(rvx.copy_t(-1)); 01008 rvc.add( rvyc); 01009 rvc=rvc.subt(rvx); 01010 01011 bdm_assert(g->dimension()==rvx._dsize(),"rvx is not described"); 01012 bdm_assert(g->dimensionc()==rvxc._dsize(),"rvxc is not described"); 01013 bdm_assert(h->dimension()==rvyc._dsize(),"rvyc is not described"); 01014 01015 bdm_assert(bm->dimensiony()==g->dimension()+h->dimension(), 01016 "Incompatible noise estimator of dimension " + 01017 num2str(bm->dimensiony()) + " does not match dimension of g and h, " + 01018 num2str(g->dimension())+" and "+ num2str(h->dimension()) ); 01019 01020 dimc = rvc._dsize(); 01021 01022 //establish datalinks 01023 x2g.set_connection(rvxc, rvx.copy_t(-1)); 01024 cond2g.set_connection(rvxc, rvc); 01025 01026 x2h.set_connection(rvyc, rvx); 01027 cond2h.set_connection(rvyc, rvc); 01028 } 01029 }; 01030 UIREGISTER(NoiseParticle); 01031 01032 01033 } 01034 #endif // KF_H 01035 01036
Generated on 2 Dec 2013 for mixpp by
1.4.7
