Changeset 225
- Timestamp:
- 01/12/09 13:28:18 (16 years ago)
- Files:
-
- 6 modified
Legend:
- Unmodified
- Added
- Removed
-
bdm/estim/libKF.cpp
r215 r225 202 202 203 203 // mat Sttm = _P->to_mat(); 204 // cout << preA <<endl; 205 // cout << "_mu:" << _mu <<endl; 204 206 205 207 // if ( !qr ( preA,Q,postA ) ) it_warning ( "QR in kalman unstable!" ); 206 208 if ( !qr ( preA,postA ) ) {it_warning ( "QR in kalman unstable!" );} 209 207 210 208 211 ( _Ry._Ch() ) =postA ( 0,dimy-1, 0,dimy-1 ); … … 229 232 // cout << "P:" <<_P.to_mat() <<endl; 230 233 // cout << "Ry:" <<_Ry.to_mat() <<endl; 231 // cout << "_K:" <<_K <<endl; 234 // cout << "_mu:" <<_mu <<endl; 235 // cout << "dt:" <<dt <<endl; 232 236 233 237 if ( evalll==true ) { //likelihood of observation y -
bdm/estim/libPF.h
r211 r225 52 52 void set_est ( const epdf &epdf0 ); 53 53 void bayes ( const vec &dt ); 54 //!access function 55 vec* __w(){return &_w;} 54 56 }; 55 57 … … 131 133 } 132 134 135 //!Access function 136 BM* _BM(int i){return Bms[i];} 133 137 //SimStr: 134 138 double SSAT; -
bdm/stat/libBM.h
r223 r225 220 220 public: 221 221 222 //! Returns the required moment of the epdf223 // virtual fnc moment ( const int order = 1 );224 222 //! Returns a sample from the density conditioned on \c cond, \f$x \sim epdf(rv|cond)\f$. \param cond is numeric value of \c rv \param ll is a return value of log-likelihood of the sample. 225 223 virtual vec samplecond ( const vec &cond, double &ll ) { … … 256 254 //!access function 257 255 epdf& _epdf() {return *ep;} 256 //!access function 257 epdf* _e() {return ep;} 258 258 }; 259 259 -
bdm/stat/libEF.cpp
r214 r225 59 59 - 0.5*xdim* ( m*log2 + 0.5* ( xdim-1 ) *log2pi ) - lg; 60 60 61 it_assert_debug (((-nkG-nkW)>-Inf) && ((-nkG-nkW)<Inf), "ARX improper");61 it_assert_debug ( ( ( -nkG-nkW ) >-Inf ) && ( ( -nkG-nkW ) <Inf ), "ARX improper" ); 62 62 return -nkG-nkW; 63 63 } … … 78 78 m ( end ) = D ( 0 ) / ( nu -nPsi -2*xdim -2 ); 79 79 return m; 80 } else { 80 } 81 else { 81 82 mat M; 82 83 mat R; 83 mean_mat (M,R);84 return cvectorize ( concat_vertical(M,R));85 } 86 87 } 88 void egiw::mean_mat (mat &M, mat&R) const {84 mean_mat ( M,R ); 85 return cvectorize ( concat_vertical ( M,R ) ); 86 } 87 88 } 89 void egiw::mean_mat ( mat &M, mat&R ) const { 89 90 const mat &L= V._L(); 90 91 const vec &D= V._D(); 91 92 int end = L.rows()-1; 92 93 93 94 ldmat ldR ( L ( 0,xdim-1,0,xdim-1 ), D ( 0,xdim-1 ) / ( nu -nPsi -2*xdim -2 ) ); //exp val of R 94 95 mat iLsub=ltuinv ( L ( xdim,end,xdim,end ) ); 95 96 96 97 // set mean value 97 98 mat Lpsi = L ( xdim,end,0,xdim-1 ); … … 105 106 106 107 for ( i=0; i<rv.count(); i++ ) { 107 GamRNG.setup ( alpha ( i ),beta ( i ) ); 108 if ( beta ( i ) >std::numeric_limits<double>::epsilon() ) { 109 GamRNG.setup ( alpha ( i ),beta ( i ) ); 110 } 111 else { 112 GamRNG.setup ( alpha ( i ),std::numeric_limits<double>::epsilon() ); 113 } 108 114 #pragma omp critical 109 115 smp ( i ) = GamRNG(); … … 136 142 } 137 143 double tmp=res-lognc();; 138 it_assert_debug (std::isfinite(tmp),"Infinite value");144 it_assert_debug ( std::isfinite ( tmp ),"Infinite value" ); 139 145 return tmp; 140 146 } … … 152 158 153 159 //MGamma 154 mgamma::mgamma ( const RV &rv,const RV &rvc ) : mEF ( rv,rvc ), epdf ( rv ) {vec* tmp; epdf._param ( tmp,_beta );};155 156 160 void mgamma::set_parameters ( double k0 ) { 157 161 k=k0; … … 255 259 } 256 260 257 void eEmp::set_samples ( 261 void eEmp::set_samples ( const epdf* epdf0 ) { 258 262 //it_assert_debug(rv==epdf0->rv(),"Wrong epdf0"); 259 263 w=1; -
bdm/stat/libEF.h
r215 r225 317 317 public : 318 318 //! Default constructor 319 egamma ( const RV &rv ) :eEF ( rv ) {};319 egamma ( const RV &rv ) :eEF ( rv ), alpha(rv.count()), beta(rv.count()) {}; 320 320 //! Sets parameters 321 321 void set_parameters ( const vec &a, const vec &b ) {alpha=a,beta=b;}; … … 328 328 void _param ( vec* &a, vec* &b ) {a=αb=β}; 329 329 vec mean() const {vec pom ( alpha ); pom/=beta; return pom;} 330 }; 331 332 /*! 333 \brief Inverse-Gamma posterior density 334 335 Multivariate inverse-Gamma density as product of independent univariate densities. 336 \f[ 337 f(x|\alpha,\beta) = \prod f(x_i|\alpha_i,\beta_i) 338 \f] 339 340 Inverse Gamma can be converted to Gamma using \[ 341 x\sim iG(a,b) => 1/x\sim G(a,1/b) 342 \] 343 This relation is used in sampling. 344 */ 345 346 class eigamma : public eEF { 347 protected: 348 //! Vector \f$\alpha\f$ 349 vec* alpha; 350 //! Vector \f$\beta\f$ (in fact it is 1/beta as used in definition of iG) 351 vec* beta; 352 //!internal egamma 353 egamma eg; 354 public : 355 //! Default constructor 356 eigamma ( const RV &rv ) :eEF ( rv ), eg(rv) {eg._param(alpha,beta);}; 357 //! Sets parameters 358 void set_parameters ( const vec &a, const vec &b ) {*alpha=a,*beta=b;}; 359 vec sample() const {return 1.0/eg.sample();}; 360 //! TODO: is it used anywhere? 361 // mat sample ( int N ) const; 362 double evallog ( const vec &val ) const {return eg.evallog(val);}; 363 double lognc () const {return eg.lognc();}; 364 //! Returns poiter to alpha and beta. Potentially dangerous: use with care! 365 void _param ( vec* &a, vec* &b ) {a=alpha;b=beta;}; 366 vec mean() const {return elem_div(*beta,*alpha-1);} 330 367 }; 331 368 /* … … 471 508 public: 472 509 //! Constructor 473 mgamma ( const RV &rv,const RV &rvc ) ;510 mgamma ( const RV &rv,const RV &rvc ): mEF ( rv,rvc ), epdf ( rv ) {vec* tmp; epdf._param ( tmp,_beta );ep=&epdf;}; 474 511 //! Set value of \c k 475 512 void set_parameters ( double k ); 476 513 void condition ( const vec &val ) {*_beta=k/val;}; 514 }; 515 516 /*! 517 \brief Inverse-Gamma random walk 518 519 Mean value, \f$\mu\f$, of this density is given by \c rvc . 520 Standard deviation of the random walk is proportional to one \f$k\f$-th the mean. 521 This is achieved by setting \f$\alpha=\mu/k+2\f$ and \f$\beta=\mu(\alpha-1)\f$. 522 523 The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$. 524 */ 525 class migamma : public mEF { 526 protected: 527 //! Internal epdf that arise by conditioning on \c rvc 528 eigamma epdf; 529 //! Constant \f$k\f$ 530 double k; 531 //! cache of epdf.beta 532 vec* _beta; 533 //! chaceh of epdf.alpha 534 vec* _alpha; 535 536 public: 537 //! Constructor 538 migamma ( const RV &rv,const RV &rvc ): mEF ( rv,rvc ), epdf ( rv ) {epdf._param ( _alpha,_beta );ep=&epdf;}; 539 //! Set value of \c k 540 void set_parameters ( double k0 ){k=k0;*_alpha=1.0/(k*k)+2;}; 541 void condition ( const vec &val ) { 542 *_beta=elem_mult(val,(*_alpha-1)); 543 }; 477 544 }; 478 545 … … 506 573 }; 507 574 575 576 /*! 577 \brief Inverse-Gamma random walk around a fixed point 578 579 Mean value, \f$\mu\f$, of this density is given by a geometric combination of \c rvc and given fixed point, \f$p\f$. \f$l\f$ is the coefficient of the geometric combimation 580 \f[ \mu = \mu_{t-1} ^{l} p^{1-l}\f] 581 582 ==== Check == vv = 583 Standard deviation of the random walk is proportional to one \f$k\f$-th the mean. 584 This is achieved by setting \f$\alpha=k\f$ and \f$\beta=k/\mu\f$. 585 586 The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$. 587 */ 588 class migamma_fix : public migamma { 589 protected: 590 //! parameter l 591 double l; 592 //! reference vector 593 vec refl; 594 public: 595 //! Constructor 596 migamma_fix ( const RV &rv,const RV &rvc ) : migamma ( rv,rvc ),refl ( rv.count() ) {}; 597 //! Set value of \c k 598 void set_parameters ( double k0 , vec ref0, double l0 ) { 599 migamma::set_parameters ( k0 ); 600 refl=pow ( ref0,1.0-l0 );l=l0; 601 }; 602 603 void condition ( const vec &val ) {vec mean=elem_mult ( refl,pow ( val,l ) ); migamma::condition(mean);}; 604 }; 508 605 //! Switch between various resampling methods. 509 606 enum RESAMPLING_METHOD { MULTINOMIAL = 0, STRATIFIED = 1, SYSTEMATIC = 3 }; -
matlab/pmsm/mpf_test_disp.m
r221 r225 1 1 close all 2 itload('../ mpf_test.it')2 itload('../../mpf_test.it') 3 3 figure(1); 4 4 for i =1:4 … … 22 22 plot(xthM(:,i)') 23 23 if i<4 24 set(gca,'XTick',[]);24 %set(gca,'XTick',[]); 25 25 else 26 26 xlabel('sample [t]');