| 639 | //! EKF with covariance R estimated by the VB method |
| 640 | class EKFvbR: public EKFfull{ |
| 641 | LOG_LEVEL(EKFvbR,logR,logQ); |
| 642 | |
| 643 | //! Statistics of the R estimator |
| 644 | mat PsiR; |
| 645 | //! Statistics of the Q estimator |
| 646 | mat PsiQ; |
| 647 | //! forgetting factor |
| 648 | double phi; |
| 649 | //! number of VB iterations |
| 650 | int niter; |
| 651 | //! degrees of freedom |
| 652 | double nu; |
| 653 | //! stabilizing element |
| 654 | mat PsiR0; |
| 655 | //! stabilizing element |
| 656 | mat PsiQ0; |
| 657 | |
| 658 | void from_setting(const Setting &set){ |
| 659 | EKFfull::from_setting(set); |
| 660 | if (!UI::get(phi,set,"phi",UI::optional)){ |
| 661 | phi = 0.99; |
| 662 | } |
| 663 | if (!UI::get(niter,set,"niter",UI::optional)){ |
| 664 | niter = 3; |
| 665 | } |
| 666 | PsiQ = Q; |
| 667 | PsiQ0 = Q; |
| 668 | PsiR = R; |
| 669 | PsiR0 = R; |
| 670 | nu = 3; |
| 671 | } |
| 672 | void log_register ( logger &L, const string &prefix ) { |
| 673 | EKFfull::log_register(L,prefix); |
| 674 | L.add_vector(log_level, logR, RV("{R }",vec_1<int>(4)), prefix); |
| 675 | L.add_vector(log_level, logQ, RV("{Q }",vec_1<int>(4)), prefix); |
| 676 | } |
| 677 | void bayes ( const vec &val, const vec &cond ) { |
| 678 | vec diffx, diffy; |
| 679 | mat Psi_vbQ; |
| 680 | mat Psi_vbR; |
| 681 | nu = phi*nu + (1-phi)*2 + 1.0; |
| 682 | |
| 683 | //save initial values of posterior |
| 684 | vec mu0=est._mu(); |
| 685 | fsqmat P0=est._R(); |
| 686 | vec xpred = pfxu->eval(mu0,cond); |
| 687 | |
| 688 | for (int i=0; i<niter; i++){ |
| 689 | est._mu() = mu0; |
| 690 | est._R() = P0; |
| 691 | |
| 692 | EKFfull::bayes(val,cond); |
| 693 | |
| 694 | diffy = val - fy._mu(); |
| 695 | Psi_vbR = phi*PsiR + (1-phi)*PsiR0+ outer_product(diffy,diffy)/*+C*mat(est._R())*C.T()*/; |
| 696 | R = Psi_vbR/(nu-2); |
| 697 | |
| 698 | diffx = est._mu() - xpred; |
| 699 | Psi_vbQ = phi*PsiQ + (1-phi)*PsiQ0+ outer_product(diffx,diffx) /*+mat(est._R()) /*+ A*mat(P0)*A.T()*/; |
| 700 | Q = Psi_vbQ/(nu-2); |
| 701 | } |
| 702 | PsiQ = Psi_vbQ; |
| 703 | PsiR = Psi_vbR; |
| 704 | // cout <<"==" << endl << Psi << endl << diff << endl << P0 << endl << ":" << Q; |
| 705 | log_level.store(logQ, vec(Q._M()._data(),4)); |
| 706 | log_level.store(logR, vec(R._M()._data(),4)); |
| 707 | |
| 708 | } |
| 709 | }; |
| 710 | UIREGISTER(EKFvbR); |
| 711 | |