| 452 | | |
| 453 | | // class MPF : public BM { |
| 454 | | // //! Introduces new option |
| 455 | | // //! \li means - meaning TODO |
| 456 | | // LOG_LEVEL(MPF,means); |
| 457 | | // protected: |
| 458 | | // //! particle filter on non-linear variable |
| 459 | | // shared_ptr<PF> pf; |
| 460 | | // //! Array of Bayesian models |
| 461 | | // Array<BM*> BMs; |
| 462 | | // |
| 463 | | // //! internal class for pdf providing composition of eEmp with external components |
| 464 | | // |
| 465 | | // class mpfepdf : public epdf { |
| 466 | | // //! pointer to particle filter |
| 467 | | // shared_ptr<PF> &pf; |
| 468 | | // //! pointer to Array of BMs |
| 469 | | // Array<BM*> &BMs; |
| 470 | | // public: |
| 471 | | // //! constructor |
| 472 | | // mpfepdf ( shared_ptr<PF> &pf0, Array<BM*> &BMs0 ) : epdf(), pf ( pf0 ), BMs ( BMs0 ) { }; |
| 473 | | // //! a variant of set parameters - this time, parameters are read from BMs and pf |
| 474 | | // void read_parameters() { |
| 475 | | // rv = concat ( pf->posterior()._rv(), BMs ( 0 )->posterior()._rv() ); |
| 476 | | // dim = pf->posterior().dimension() + BMs ( 0 )->posterior().dimension(); |
| 477 | | // bdm_assert_debug ( dim == rv._dsize(), "Wrong name " ); |
| 478 | | // } |
| 479 | | // vec mean() const; |
| 480 | | // |
| 481 | | // vec variance() const; |
| 482 | | // |
| 483 | | // void qbounds ( vec &lb, vec &ub, double perc = 0.95 ) const; |
| 484 | | // |
| 485 | | // vec sample() const NOT_IMPLEMENTED(0); |
| 486 | | // |
| 487 | | // double evallog ( const vec &val ) const NOT_IMPLEMENTED(0); |
| 488 | | // }; |
| 489 | | // |
| 490 | | // //! Density joining PF.est with conditional parts |
| 491 | | // mpfepdf jest; |
| 492 | | // |
| 493 | | // //! datalink from global yt and cond (Up) to BMs yt and cond (Down) |
| 494 | | // datalink_m2m this2bm; |
| 495 | | // //! datalink from global yt and cond (Up) to PFs yt and cond (Down) |
| 496 | | // datalink_m2m this2pf; |
| 497 | | // //!datalink from PF part to BM |
| 498 | | // datalink_part pf2bm; |
| 499 | | // |
| 500 | | // public: |
| 501 | | // //! Default constructor. |
| 502 | | // MPF () : jest ( pf, BMs ) {}; |
| 503 | | // //! set all parameters at once |
| 504 | | // void set_pf ( shared_ptr<pdf> par0, int n0, RESAMPLING_METHOD rm = SYSTEMATIC ) { |
| 505 | | // if (!pf) pf=new PF; |
| 506 | | // pf->set_model ( par0, par0 ); // <=== nasty!!! |
| 507 | | // pf->set_parameters ( n0, rm ); |
| 508 | | // pf->set_rv(par0->_rv()); |
| 509 | | // BMs.set_length ( n0 ); |
| 510 | | // } |
| 511 | | // //! set a prototype of BM, copy it to as many times as there is particles in pf |
| 512 | | // void set_BM ( const BM &BMcond0 ) { |
| 513 | | // |
| 514 | | // int n = pf->__w().length(); |
| 515 | | // BMs.set_length ( n ); |
| 516 | | // // copy |
| 517 | | // //BMcond0 .condition ( pf->posterior()._sample ( 0 ) ); |
| 518 | | // for ( int i = 0; i < n; i++ ) { |
| 519 | | // BMs ( i ) = (BM*) BMcond0._copy(); |
| 520 | | // } |
| 521 | | // }; |
| 522 | | // |
| 523 | | // void bayes ( const vec &yt, const vec &cond ); |
| 524 | | // |
| 525 | | // const epdf& posterior() const { |
| 526 | | // return jest; |
| 527 | | // } |
| 528 | | // |
| 529 | | // //!Access function |
| 530 | | // const BM* _BM ( int i ) { |
| 531 | | // return BMs ( i ); |
| 532 | | // } |
| 533 | | // PF& _pf() {return *pf;} |
| 534 | | // |
| 535 | | // |
| 536 | | // virtual double logpred ( const vec &yt ) const NOT_IMPLEMENTED(0); |
| 537 | | // |
| 538 | | // virtual epdf* epredictor() const NOT_IMPLEMENTED(NULL); |
| 539 | | // |
| 540 | | // virtual pdf* predictor() const NOT_IMPLEMENTED(NULL); |
| 541 | | // |
| 542 | | // |
| 543 | | // /*! configuration structure for basic PF |
| 544 | | // \code |
| 545 | | // BM = BM_class; // Bayesian filtr for analytical part of the model |
| 546 | | // parameter_pdf = pdf_class; // transitional pdf for non-parametric part of the model |
| 547 | | // prior = epdf_class; // prior probability density |
| 548 | | // --- optional --- |
| 549 | | // n = 10; // number of particles |
| 550 | | // resmethod = 'systematic', or 'multinomial', or 'stratified' |
| 551 | | // // resampling method |
| 552 | | // \endcode |
| 553 | | // */ |
| 554 | | // void from_setting ( const Setting &set ) { |
| 555 | | // BM::from_setting( set ); |
| 556 | | // |
| 557 | | // shared_ptr<pdf> par = UI::build<pdf> ( set, "parameter_pdf", UI::compulsory ); |
| 558 | | // |
| 559 | | // pf = new PF; |
| 560 | | // // prior must be set before BM |
| 561 | | // pf->prior_from_set ( set ); |
| 562 | | // pf->resmethod_from_set ( set ); |
| 563 | | // pf->set_model ( par, par ); // too hackish! |
| 564 | | // |
| 565 | | // shared_ptr<BM> BM0 = UI::build<BM> ( set, "BM", UI::compulsory ); |
| 566 | | // set_BM ( *BM0 ); |
| 567 | | // |
| 568 | | // //set drv |
| 569 | | // //??set_yrv(concat(BM0->_yrv(),u) ); |
| 570 | | // set_yrv ( BM0->_yrv() ); |
| 571 | | // rvc = BM0->_rvc().subt ( par->_rv() ); |
| 572 | | // //find potential input - what remains in rvc when we subtract rv |
| 573 | | // RV u = par->_rvc().subt ( par->_rv().copy_t ( -1 ) ); |
| 574 | | // rvc.add ( u ); |
| 575 | | // dimc = rvc._dsize(); |
| 576 | | // validate(); |
| 577 | | // } |
| 578 | | // |
| 579 | | // void validate() { |
| 580 | | // BM::validate(); |
| 581 | | // try { |
| 582 | | // pf->validate(); |
| 583 | | // } catch ( std::exception ) { |
| 584 | | // throw UIException ( "Error in PF part of MPF:" ); |
| 585 | | // } |
| 586 | | // jest.read_parameters(); |
| 587 | | // this2bm.set_connection ( BMs ( 0 )->_yrv(), BMs ( 0 )->_rvc(), yrv, rvc ); |
| 588 | | // this2pf.set_connection ( pf->_yrv(), pf->_rvc(), yrv, rvc ); |
| 589 | | // pf2bm.set_connection ( BMs ( 0 )->_rvc(), pf->posterior()._rv() ); |
| 590 | | // } |
| 591 | | // }; |
| 592 | | // UIREGISTER ( MPF ); |
| 593 | | |
| 594 | | /*! ARXg for estimation of state-space variances |
| 595 | | */ |
| 596 | | // class MPF_ARXg :public BM{ |
| 597 | | // protected: |
| 598 | | // shared_ptr<PF> pf; |
| 599 | | // //! pointer to Array of BMs |
| 600 | | // Array<ARX*> BMso; |
| 601 | | // //! pointer to Array of BMs |
| 602 | | // Array<ARX*> BMsp; |
| 603 | | // //!parameter evolution |
| 604 | | // shared_ptr<fnc> g; |
| 605 | | // //!observation function |
| 606 | | // shared_ptr<fnc> h; |
| 607 | | // |
| 608 | | // public: |
| 609 | | // void bayes(const vec &yt, const vec &cond ); |
| 610 | | // void from_setting(const Setting &set) ; |
| 611 | | // void validate() { |
| 612 | | // bdm_assert(g->dimensionc()==g->dimension(),"not supported yet"); |
| 613 | | // bdm_assert(h->dimensionc()==g->dimension(),"not supported yet"); |
| 614 | | // } |
| 615 | | // |
| 616 | | // double logpred(const vec &cond) const NOT_IMPLEMENTED(0.0); |
| 617 | | // epdf* epredictor() const NOT_IMPLEMENTED(NULL); |
| 618 | | // pdf* predictor() const NOT_IMPLEMENTED(NULL); |
| 619 | | // const epdf& posterior() const {return pf->posterior();}; |
| 620 | | // |
| 621 | | // void log_register( logger &L, const string &prefix ){ |
| 622 | | // BM::log_register(L,prefix); |
| 623 | | // registered_logger->ids.set_size ( 3 ); |
| 624 | | // registered_logger->ids(1)= L.add_vector(RV("Q",dimension()*dimension()), prefix+L.prefix_sep()+"Q"); |
| 625 | | // registered_logger->ids(2)= L.add_vector(RV("R",dimensiony()*dimensiony()), prefix+L.prefix_sep()+"R"); |
| 626 | | // |
| 627 | | // }; |
| 628 | | // void log_write() const { |
| 629 | | // BM::log_write(); |
| 630 | | // mat mQ=zeros(dimension(),dimension()); |
| 631 | | // mat pom=zeros(dimension(),dimension()); |
| 632 | | // mat mR=zeros(dimensiony(),dimensiony()); |
| 633 | | // mat pom2=zeros(dimensiony(),dimensiony()); |
| 634 | | // mat dum; |
| 635 | | // const vec w=pf->posterior()._w(); |
| 636 | | // for (int i=0; i<w.length(); i++){ |
| 637 | | // BMsp(i)->posterior().mean_mat(dum,pom); |
| 638 | | // mQ += w(i) * pom; |
| 639 | | // BMso(i)->posterior().mean_mat(dum,pom2); |
| 640 | | // mR += w(i) * pom2; |
| 641 | | // |
| 642 | | // } |
| 643 | | // registered_logger->L.log_vector ( registered_logger->ids ( 1 ), cvectorize(mQ) ); |
| 644 | | // registered_logger->L.log_vector ( registered_logger->ids ( 2 ), cvectorize(mR) ); |
| 645 | | // |
| 646 | | // } |
| 647 | | // }; |
| 648 | | // UIREGISTER(MPF_ARXg); |
| | 472 | class NoiseParticle : public MarginalizedParticleBase{ |
| | 473 | protected: |
| | 474 | //! function transforming xt, ut -> x_t+1 |
| | 475 | shared_ptr<fnc> g; // pdf for non-linear part |
| | 476 | //! function transforming xt,ut -> yt |
| | 477 | shared_ptr<fnc> h; // pdf for non-linear part |
| | 478 | |
| | 479 | RV rvv; |
| | 480 | |
| | 481 | RV rvx; |
| | 482 | RV rvxc; |
| | 483 | RV rvyc; |
| | 484 | |
| | 485 | //!link from condition to f |
| | 486 | datalink_part cond2g; |
| | 487 | //!link from condition to h |
| | 488 | datalink_part cond2h; |
| | 489 | //!link from xt to f |
| | 490 | datalink_part x2g; |
| | 491 | //!link from xt to h |
| | 492 | datalink_part x2h; |
| | 493 | |
| | 494 | public: |
| | 495 | BM* _copy() const{return new NoiseParticle();}; |
| | 496 | void bayes(const vec &dt, const vec &cond){ |
| | 497 | epdf* pred_vw=bm->epredictor(); |
| | 498 | shared_ptr<epdf> pred_v = pred_vw->marginal(rvv); |
| | 499 | |
| | 500 | vec vt=pred_v->sample(); |
| | 501 | |
| | 502 | //new sample |
| | 503 | vec &xtm=est_emp.point; |
| | 504 | vec g_args(g->dimensionc()); |
| | 505 | x2g.filldown(xtm,g_args); |
| | 506 | cond2g.filldown(cond,g_args); |
| | 507 | vec xt = g->eval(g_args) + vt; |
| | 508 | |
| | 509 | // residue of observation |
| | 510 | vec h_args(h->dimensionc()); |
| | 511 | x2h.filldown(xt,h_args); |
| | 512 | cond2h.filldown(cond,h_args); |
| | 513 | vec wt = dt-h->eval(h_args); |
| | 514 | // the vector [v_t,w_t] is now complete |
| | 515 | bm->bayes(concat(vt,wt)); |
| | 516 | ll=bm->_ll(); |
| | 517 | } |
| | 518 | void from_setting(const Setting &set){ |
| | 519 | MarginalizedParticleBase::from_setting(set); //reads bm, yrv,rvc, bm_rv, etc... |
| | 520 | |
| | 521 | UI::get(g,set,"g",UI::compulsory); |
| | 522 | UI::get(h,set,"h",UI::compulsory); |
| | 523 | RV rvx; |
| | 524 | UI::get(rvx,set,"rvx",UI::compulsory); |
| | 525 | est_emp.set_rv(rvx); |
| | 526 | |
| | 527 | RV rvxc; |
| | 528 | UI::get(rvxc,set,"rvxc",UI::compulsory); |
| | 529 | RV rvyc; |
| | 530 | UI::get(rvyc,set,"rvyc",UI::compulsory); |
| | 531 | } |
| | 532 | void validate(){ |
| | 533 | MarginalizedParticleBase::validate(); |
| | 534 | //check dimensions |
| | 535 | rvc = rvxc; |
| | 536 | rvc.add( rvyc); |
| | 537 | rvc.subt(rvx); |
| | 538 | |
| | 539 | //establish datalinks |
| | 540 | x2g.set_connection(rvxc, rvx.copy_t(-1)); |
| | 541 | cond2g.set_connection(rvxc, rvc); |
| | 542 | |
| | 543 | x2h.set_connection(rvyc, rvx); |
| | 544 | cond2h.set_connection(rvyc, rvc); |
| | 545 | } |
| | 546 | }; |
| | 547 | UIREGISTER(NoiseParticle); |
| | 548 | |
