| | 627 | /*! Marginalized particle for state-space models with unknown parameters of distribuution of residues on \f$v_t\f$ and \f$ w_t \f$. |
| | 628 | |
| | 629 | \f{eqnarray*}{ |
| | 630 | x_t &=& g(x_{t-1}) + v_t,\\ |
| | 631 | y_t &= &h(x_t)+w_t, |
| | 632 | \f} |
| | 633 | |
| | 634 | This particle is a only a shell creating the residues calling internal estimator of their parameters. The internal estimator can be of any compatible type, e.g. ARX for Gaussian residues with unknown mean and variance. |
| | 635 | |
| | 636 | */ |
| | 637 | class NoiseParticleXY : public BM { |
| | 638 | //! discrte particle |
| | 639 | dirac est_emp; |
| | 640 | //! internal Bayes Model |
| | 641 | shared_ptr<BM> bmx; |
| | 642 | shared_ptr<BM> bmy; |
| | 643 | |
| | 644 | //! \brief Internal class for custom posterior - product of empirical and exact part |
| | 645 | class eprod_3:public eprod_base { |
| | 646 | protected: |
| | 647 | NoiseParticleXY ∓ |
| | 648 | public: |
| | 649 | eprod_3(NoiseParticleXY &m):mp(m) {} |
| | 650 | const epdf* factor(int i) const { |
| | 651 | if (i==0) return &mp.bmx->posterior() ; |
| | 652 | if(i==1) return &mp.bmy->posterior(); |
| | 653 | return &mp.est_emp; |
| | 654 | } |
| | 655 | const int no_factors() const { |
| | 656 | return 3; |
| | 657 | } |
| | 658 | } est; |
| | 659 | |
| | 660 | protected: |
| | 661 | //! function transforming xt, ut -> x_t+1 |
| | 662 | shared_ptr<fnc> g; // pdf for non-linear part |
| | 663 | //! function transforming xt,ut -> yt |
| | 664 | shared_ptr<fnc> h; // pdf for non-linear part |
| | 665 | |
| | 666 | RV rvx; |
| | 667 | RV rvxc; |
| | 668 | RV rvyc; |
| | 669 | |
| | 670 | //!link from condition to f |
| | 671 | datalink_part cond2g; |
| | 672 | //!link from condition to h |
| | 673 | datalink_part cond2h; |
| | 674 | //!link from xt to f |
| | 675 | datalink_part x2g; |
| | 676 | //!link from xt to h |
| | 677 | datalink_part x2h; |
| | 678 | |
| | 679 | public: |
| | 680 | NoiseParticleXY():est(*this) {}; |
| | 681 | NoiseParticleXY(const NoiseParticleXY &m2):BM(m2),est(*this),h(m2.h),g(m2.g), rvx(m2.rvx),rvxc(m2.rvxc),rvyc(m2.rvyc) { |
| | 682 | bmx = m2.bmx->_copy(); |
| | 683 | bmy = m2.bmy->_copy(); |
| | 684 | est_emp = m2.est_emp; |
| | 685 | //est.validate(); |
| | 686 | validate(); |
| | 687 | }; |
| | 688 | |
| | 689 | const eprod_3& posterior() const { |
| | 690 | return est; |
| | 691 | } |
| | 692 | |
| | 693 | void set_prior(const epdf *pdf0) { |
| | 694 | const eprod *ep=dynamic_cast<const eprod*>(pdf0); |
| | 695 | if (ep) { // full prior |
| | 696 | bdm_assert(ep->no_factors()==2,"Incompatible prod"); |
| | 697 | bmx->set_prior(ep->factor(0)); |
| | 698 | bmy->set_prior(ep->factor(1)); |
| | 699 | est_emp.set_point(ep->factor(2)->sample()); |
| | 700 | } else { |
| | 701 | // assume prior is only for emp; |
| | 702 | est_emp.set_point(pdf0->sample()); |
| | 703 | } |
| | 704 | } |
| | 705 | |
| | 706 | BM* _copy() const { |
| | 707 | return new NoiseParticleXY(*this); |
| | 708 | }; |
| | 709 | |
| | 710 | void bayes(const vec &dt, const vec &cond) { |
| | 711 | //shared_ptr<epdf> pred_v=bm->epredictor(); |
| | 712 | |
| | 713 | //vec vt=pred_v->sample(); |
| | 714 | vec vt = bmx->samplepred(); |
| | 715 | |
| | 716 | //new sample |
| | 717 | vec &xtm=est_emp.point; |
| | 718 | vec g_args(g->dimensionc()); |
| | 719 | x2g.filldown(xtm,g_args); |
| | 720 | cond2g.filldown(cond,g_args); |
| | 721 | vec xt = g->eval(g_args) + vt; |
| | 722 | est_emp.point=xt; |
| | 723 | |
| | 724 | // the vector [v_t] updates bm, |
| | 725 | bmx->bayes(vt); |
| | 726 | |
| | 727 | // residue of observation |
| | 728 | vec h_args(h->dimensionc()); |
| | 729 | x2h.filldown(xt,h_args); |
| | 730 | cond2h.filldown(cond,h_args); |
| | 731 | |
| | 732 | bmy->bayes(h->eval(h_args)-dt); |
| | 733 | ll= bmy->_ll(); |
| | 734 | } |
| | 735 | void from_setting(const Setting &set) { |
| | 736 | BM::from_setting(set); //reads bm, yrv,rvc, bm_rv, etc... |
| | 737 | bmx = UI::build<BM>(set,"bmx",UI::compulsory); |
| | 738 | |
| | 739 | bmy = UI::build<BM>(set,"bmy",UI::compulsory); |
| | 740 | g=UI::build<fnc>(set,"g",UI::compulsory); |
| | 741 | h=UI::build<fnc>(set,"h",UI::compulsory); |
| | 742 | UI::get(rvx,set,"rvx",UI::compulsory); |
| | 743 | est_emp.set_rv(rvx); |
| | 744 | |
| | 745 | UI::get(rvxc,set,"rvxc",UI::compulsory); |
| | 746 | UI::get(rvyc,set,"rvyc",UI::compulsory); |
| | 747 | |
| | 748 | } |
| | 749 | void to_setting (Setting &set) const { |
| | 750 | BM::to_setting(set); //reads bm, yrv,rvc, bm_rv, etc... |
| | 751 | UI::save(g,set,"g"); |
| | 752 | UI::save(h,set,"h"); |
| | 753 | UI::save(bmx,set,"bmx"); |
| | 754 | UI::save(bmy,set,"bmy"); |
| | 755 | } |
| | 756 | void validate() { |
| | 757 | BM::validate(); |
| | 758 | |
| | 759 | dimy = h->dimension(); |
| | 760 | // bmx->set_yrv(rvx); |
| | 761 | // bmy-> |
| | 762 | |
| | 763 | est_emp.set_rv(rvx); |
| | 764 | est_emp.set_dim(rvx._dsize()); |
| | 765 | est.validate(); |
| | 766 | // |
| | 767 | //check dimensions |
| | 768 | rvc = rvxc.subt(rvx.copy_t(-1)); |
| | 769 | rvc.add( rvyc); |
| | 770 | rvc=rvc.subt(rvx); |
| | 771 | |
| | 772 | bdm_assert(g->dimension()==rvx._dsize(),"rvx is not described"); |
| | 773 | bdm_assert(g->dimensionc()==rvxc._dsize(),"rvxc is not described"); |
| | 774 | bdm_assert(h->dimensionc()==rvyc._dsize(),"rvyc is not described"); |
| | 775 | |
| | 776 | bdm_assert(bmx->dimensiony()==g->dimension(), |
| | 777 | "Incompatible noise estimator of dimension " + |
| | 778 | num2str(bmx->dimensiony()) + " does not match dimension of g , " + |
| | 779 | num2str(g->dimension()) |
| | 780 | ); |
| | 781 | |
| | 782 | dimc = rvc._dsize(); |
| | 783 | |
| | 784 | //establish datalinks |
| | 785 | x2g.set_connection(rvxc, rvx.copy_t(-1)); |
| | 786 | cond2g.set_connection(rvxc, rvc); |
| | 787 | |
| | 788 | x2h.set_connection(rvyc, rvx); |
| | 789 | cond2h.set_connection(rvyc, rvc); |
| | 790 | } |
| | 791 | |
| | 792 | }; |
| | 793 | UIREGISTER(NoiseParticleXY); |
| | 794 | |
