394 | | if (log_level[logweights]) { |
395 | | L.add_vector( log_level, logweights, RV ( particles.length()), prefix); |
396 | | } |
397 | | if (log_level[logmeans]) { |
| 394 | if (log_level[logweights]) { |
| 395 | L.add_vector( log_level, logweights, RV ( particles.length()), prefix); |
| 396 | } |
| 397 | if (log_level[logNeff]) { |
| 398 | L.add_vector( log_level, logNeff, RV ( 1), prefix); |
| 399 | } |
| 400 | if (log_level[logmeans]) { |
550 | | //vec vt=pred_v->sample(); |
551 | | vec vt = bm->samplepred(); |
| 553 | vec xt; |
| 554 | vec xthat; |
| 555 | //vec vt=pred_v->sample(); |
| 556 | // vec vt = bm->samplepred(); |
| 557 | if (1){ |
| 558 | //vec vt=pred_v->sample(); |
| 559 | int dim_x = bm->dimensiony(); |
| 560 | int dim_y = dimensiony(); |
| 561 | |
| 562 | vec post_mean =bm->posterior().mean(); |
| 563 | mat SigV; |
| 564 | if (post_mean.length()==dim_x){ |
| 565 | SigV=diag(vec(post_mean._data(), dim_x)); |
| 566 | } else { |
| 567 | SigV=mat(post_mean._data(), dim_x, dim_x); |
| 568 | } |
| 569 | |
| 570 | vec &xtm=est_emp.point; |
| 571 | vec g_args(g->dimensionc()); |
| 572 | x2g.filldown(xtm,g_args); |
| 573 | cond2g.filldown(cond,g_args); |
| 574 | xthat = g->eval(g_args); |
| 575 | |
| 576 | mat IC = eye(dim_x+dim_y); |
| 577 | IC.set_submatrix(0,0,eye(dim_x)); |
| 578 | mat SigJo=eye(dim_x + dim_y); |
| 579 | SigJo.set_submatrix(0,0,SigV); |
| 580 | mlnorm<ldmat>* mlfy=dynamic_cast<mlnorm<ldmat>* >(fy.get()); |
| 581 | |
| 582 | SigJo.set_submatrix(dim_x,dim_x, mlfy->_R()); |
| 583 | IC.set_submatrix(dim_x,0,-mlfy->_A()); |
| 584 | |
| 585 | mat iIC=inv(IC); |
| 586 | enorm<chmat> Jo; |
| 587 | Jo._mu()=iIC*concat(xthat, zeros(dim_y)); |
| 588 | Jo._R()=iIC*SigJo*iIC.T(); |
| 589 | Jo.set_rv(concat(bm->_yrv(),yrv)); |
| 590 | Jo.validate(); |
| 591 | shared_ptr<pdf> Cond=Jo.condition(bm->_yrv()); |
| 592 | //new sample |
| 593 | |
| 594 | xt = Cond->samplecond(dt);// + vt; |
| 595 | } else{ |
| 596 | |
| 856 | class NoiseParticleXYprop: public NoiseParticleXY{ |
| 857 | public: |
| 858 | BM* _copy() const { |
| 859 | return new NoiseParticleXYprop(*this); |
| 860 | }; |
| 861 | void bayes(const vec &dt, const vec &cond) { |
| 862 | int dim_x = bmx->dimensiony(); |
| 863 | int dim_y = bmy->dimensiony(); |
| 864 | |
| 865 | //vec vt=pred_v->sample(); |
| 866 | mat SigV=mat(bmx->posterior().mean()._data(), dim_x, dim_x); |
| 867 | mat SigW=mat(bmy->posterior().mean()._data(), dim_y, dim_y); |
| 868 | |
| 869 | vec &xtm=est_emp.point; |
| 870 | vec g_args(g->dimensionc()); |
| 871 | x2g.filldown(xtm,g_args); |
| 872 | cond2g.filldown(cond,g_args); |
| 873 | vec xthat = g->eval(g_args); |
| 874 | |
| 875 | mat IC = eye(dim_x+dim_y); |
| 876 | IC.set_submatrix(0,0,eye(dim_x)); |
| 877 | mat SigJo=eye(dim_x + dim_y); |
| 878 | SigJo.set_submatrix(0,0,SigV); |
| 879 | SigJo.set_submatrix(dim_x,dim_x, 10*SigW); |
| 880 | |
| 881 | diffbifn* hb=dynamic_cast<bilinfn*>(h.get()); |
| 882 | mat C=zeros(dim_y,dim_x); |
| 883 | if (hb){ |
| 884 | hb->dfdx_cond(xtm, cond, C,true); |
| 885 | IC.set_submatrix(dim_x,0,-C); |
| 886 | } |
| 887 | |
| 888 | mat iIC=inv(IC); |
| 889 | enorm<chmat> Jo; |
| 890 | Jo._mu()=iIC*concat(xthat, zeros(dim_y)); |
| 891 | Jo._R()=iIC*SigJo*iIC.T(); |
| 892 | Jo.set_rv(concat(bmx->_yrv(), bmy->_yrv())); |
| 893 | Jo.validate(); |
| 894 | shared_ptr<pdf> Cond=Jo.condition(bmx->_yrv()); |
| 895 | //new sample |
| 896 | |
| 897 | vec xt = Cond->samplecond(dt);// + vt; |
| 898 | |
| 899 | est_emp.point=xt; |
| 900 | |
| 901 | // the vector [v_t] updates bm, |
| 902 | bmx->bayes(xt-xthat); |
| 903 | |
| 904 | // residue of observation |
| 905 | vec h_args(h->dimensionc()); |
| 906 | x2h.filldown(xt,h_args); |
| 907 | cond2h.filldown(cond,h_args); |
| 908 | |
| 909 | vec z_y =h->eval(h_args)-dt; |
| 910 | // ARX *abm = dynamic_cast<ARX*>(bmy.get()); |
| 911 | /* double ll2; |
| 912 | i f (abm){ //ARX* |
| 913 | shared_ptr<epdf> pr_y(abm->epredictor_student(empty_vec)); |
| 914 | ll2=pr_y->evallog(z_y); |
| 915 | } else{ |
| 916 | shared_ptr<epdf> pr_y(bmy->epredictor(empty_vec)); |
| 917 | ll2=pr_y->evallog(z_y); |
| 918 | }*/ |
| 919 | |
| 920 | bmy->bayes(z_y); |
| 921 | // test _lls |
| 922 | ll= bmy->_ll(); |
| 923 | } |
| 924 | }; |
| 925 | UIREGISTER(NoiseParticleXYprop); |
| 926 | |