Changeset 178 for bdm/stat/libEF.h
- Timestamp:
- 10/15/08 19:04:09 (16 years ago)
- Files:
-
- 1 modified
Legend:
- Unmodified
- Added
- Removed
-
bdm/stat/libEF.h
r176 r178 44 44 virtual double lognc() const =0; 45 45 //!TODO decide if it is really needed 46 virtual void dupdate ( mat &v ) {it_error ( "Not implem neted" );};46 virtual void dupdate ( mat &v ) {it_error ( "Not implemented" );}; 47 47 //!Evaluate normalized log-probability 48 virtual double evalpdflog_nn ( const vec &val ) const{it_error ( "Not implem neted" );return 0.0;};48 virtual double evalpdflog_nn ( const vec &val ) const{it_error ( "Not implemented" );return 0.0;}; 49 49 //!Evaluate normalized log-probability 50 50 virtual double evalpdflog ( const vec &val ) const {return evalpdflog_nn ( val )-lognc();} … … 90 90 //original Bayes 91 91 void bayes ( const vec &dt ); 92 //!Flatten the posterior 93 virtual void flatten ( BMEF * B ) {it_error ( "Not implemented" );} 94 }; 92 //!Flatten the posterior according to the given BMEF (of the same type!) 93 virtual void flatten ( const BMEF * B ) {it_error ( "Not implemented" );} 94 //!Flatten the posterior as if to keep nu0 data 95 virtual void flatten ( double nu0 ) {it_error ( "Not implemented" );} 96 }; 97 98 template<class sq_T> 99 class mlnorm; 95 100 96 101 /*! … … 100 105 */ 101 106 template<class sq_T> 102 103 107 class enorm : public eEF { 104 108 protected: … … 110 114 int dim; 111 115 public: 112 // enorm() :eEF() {};113 116 //!Default constructor 114 enorm ( RV &rv );117 enorm ( const RV &rv ); 115 118 //! Set mean value \c mu and covariance \c R 116 119 void set_parameters ( const vec &mu,const sq_T &R ); 117 // ! tupdate in exponential form (not really handy)118 void tupdate ( double phi, mat &vbar, double nubar );120 ////! tupdate in exponential form (not really handy) 121 //void tupdate ( double phi, mat &vbar, double nubar ); 119 122 //! dupdate in exponential form (not really handy) 120 123 void dupdate ( mat &v,double nu=1.0 ); … … 124 127 mat sample ( int N ) const; 125 128 double eval ( const vec &val ) const ; 126 double evalpdflog ( const vec &val ) const;129 double evalpdflog_nn ( const vec &val ) const; 127 130 double lognc () const; 128 131 vec mean() const {return mu;} 129 132 mlnorm<sq_T>* condition ( const RV &rvn ); 133 enorm<sq_T>* marginal ( const RV &rv ); 130 134 //Access methods 131 135 //! returns a pointer to the internal mean value. Use with Care! … … 139 143 140 144 //! access method 141 mat getR () {return R.to_mat();}145 // mat getR () {return R.to_mat();} 142 146 }; 143 147 … … 215 219 }; 216 220 //!access function 217 vec& _beta() {return beta;}221 vec& _beta() {return beta;} 218 222 //!Set internal parameters 219 void set_parameters (const vec &beta0){220 if (beta0.length()!=beta.length()){221 it_assert_debug (rv.length()==1,"Undefined");222 rv.set_size (0,beta0.length());223 void set_parameters ( const vec &beta0 ) { 224 if ( beta0.length() !=beta.length() ) { 225 it_assert_debug ( rv.length() ==1,"Undefined" ); 226 rv.set_size ( 0,beta0.length() ); 223 227 } 224 228 beta= beta0; … … 258 262 return pred.lognc()-lll; 259 263 } 260 void flatten ( BMEF* B ) {261 eDirich* E=dynamic_cast<eDirich*> ( B );264 void flatten ( const BMEF* B ) { 265 const eDirich* E=dynamic_cast<const eDirich*> ( B ); 262 266 // sum(beta) should be equal to sum(B.beta) 263 const vec &Eb= E->_beta();267 const vec &Eb=const_cast<eDirich*> ( E )->_beta(); 264 268 est.pow ( sum ( beta ) /sum ( Eb ) ); 265 269 if ( evalll ) {last_lognc=est.lognc();} … … 267 271 const epdf& _epdf() const {return est;}; 268 272 void set_parameters ( const vec &beta0 ) { 269 est.set_parameters (beta0);273 est.set_parameters ( beta0 ); 270 274 rv = est._rv(); 271 if (evalll){last_lognc=est.lognc();}275 if ( evalll ) {last_lognc=est.lognc();} 272 276 } 273 277 }; … … 359 363 \brief Normal distributed linear function with linear function of mean value; 360 364 361 Mean value \f$mu=A*rvc \f$.365 Mean value \f$mu=A*rvc+mu_0\f$. 362 366 */ 363 367 template<class sq_T> … … 366 370 enorm<sq_T> epdf; 367 371 mat A; 372 vec mu_const; 368 373 vec& _mu; //cached epdf.mu; 369 374 public: 370 375 //! Constructor 371 mlnorm ( RV &rv,RV &rvc );376 mlnorm (const RV &rv, const RV &rvc ); 372 377 //! Set \c A and \c R 373 void set_parameters ( const mat &A, const sq_T &R );378 void set_parameters ( const mat &A, const vec &mu0, const sq_T &R ); 374 379 //!Generate one sample of the posterior 375 vec samplecond ( vec &cond, double &lik );380 vec samplecond (const vec &cond, double &lik ); 376 381 //!Generate matrix of samples of the posterior 377 mat samplecond ( vec &cond, vec &lik, int n );382 mat samplecond (const vec &cond, vec &lik, int n ); 378 383 //! Set value of \c rvc . Result of this operation is stored in \c epdf use function \c _ep to access it. 379 void condition ( vec &cond );384 void condition (const vec &cond ); 380 385 }; 381 386 … … 453 458 //! Default constructor 454 459 eEmp ( const RV &rv0 ,int n0 ) :epdf ( rv0 ),n ( n0 ),w ( n ),samples ( n ) {}; 460 //! Set samples and weights 461 void set_parameters ( const vec &w0, const epdf* pdf0 ); 455 462 //! Set sample 456 void set_ parameters ( const vec &w0,epdf* pdf0 );463 void set_samples ( const epdf* pdf0 ); 457 464 //! Potentially dangerous, use with care. 458 465 vec& _w() {return w;}; … … 476 483 477 484 template<class sq_T> 478 enorm<sq_T>::enorm ( RV &rv ) :eEF ( rv ), mu ( rv.count() ),R ( rv.count() ),dim ( rv.count() ) {};485 enorm<sq_T>::enorm ( const RV &rv ) :eEF ( rv ), mu ( rv.count() ),R ( rv.count() ),dim ( rv.count() ) {}; 479 486 480 487 template<class sq_T> … … 490 497 }; 491 498 492 template<class sq_T>493 void enorm<sq_T>::tupdate ( double phi, mat &vbar, double nubar ) {494 //495 };499 // template<class sq_T> 500 // void enorm<sq_T>::tupdate ( double phi, mat &vbar, double nubar ) { 501 // // 502 // }; 496 503 497 504 template<class sq_T> … … 531 538 532 539 template<class sq_T> 533 double enorm<sq_T>::evalpdflog ( const vec &val ) const {540 double enorm<sq_T>::evalpdflog_nn ( const vec &val ) const { 534 541 // 1.83787706640935 = log(2pi) 535 return -0.5* ( +R.invqform ( mu-val ) )- lognc();542 return -0.5* ( R.invqform ( mu-val ) );// - lognc(); 536 543 }; 537 544 … … 539 546 inline double enorm<sq_T>::lognc () const { 540 547 // 1.83787706640935 = log(2pi) 541 return -0.5* ( R.cols() * 1.83787706640935 +R.logdet() );542 }; 543 544 template<class sq_T> 545 mlnorm<sq_T>::mlnorm ( RV &rv0,RV &rvc0 ) :mEF ( rv0,rvc0 ),epdf ( rv0 ),A ( rv0.count(),rv0.count() ),_mu ( epdf._mu() ) {548 return 0.5* ( R.cols() * 1.83787706640935 +R.logdet() ); 549 }; 550 551 template<class sq_T> 552 mlnorm<sq_T>::mlnorm ( const RV &rv0, const RV &rvc0 ) :mEF ( rv0,rvc0 ),epdf ( rv0 ),A ( rv0.count(),rv0.count() ),_mu ( epdf._mu() ) { 546 553 ep =&epdf; 547 554 } 548 555 549 556 template<class sq_T> 550 void mlnorm<sq_T>::set_parameters ( const mat &A0, const sq_T &R0 ) {557 void mlnorm<sq_T>::set_parameters ( const mat &A0, const vec &mu0, const sq_T &R0 ) { 551 558 epdf.set_parameters ( zeros ( rv.count() ),R0 ); 552 559 A = A0; 560 mu_const = mu0; 553 561 } 554 562 555 563 template<class sq_T> 556 vec mlnorm<sq_T>::samplecond ( vec &cond, double &lik ) {564 vec mlnorm<sq_T>::samplecond (const vec &cond, double &lik ) { 557 565 this->condition ( cond ); 558 566 vec smp = epdf.sample(); … … 562 570 563 571 template<class sq_T> 564 mat mlnorm<sq_T>::samplecond ( vec &cond, vec &lik, int n ) {572 mat mlnorm<sq_T>::samplecond (const vec &cond, vec &lik, int n ) { 565 573 int i; 566 574 int dim = rv.count(); … … 579 587 580 588 template<class sq_T> 581 void mlnorm<sq_T>::condition ( vec &cond ) {582 _mu = A*cond ;589 void mlnorm<sq_T>::condition (const vec &cond ) { 590 _mu = A*cond + mu_const; 583 591 //R is already assigned; 584 592 } 585 593 594 template<class sq_T> 595 enorm<sq_T>* enorm<sq_T>::marginal ( const RV &rvn ) { 596 ivec irvn = rvn.dataind ( rv ); 597 598 sq_T Rn ( R,irvn ); 599 enorm<sq_T>* tmp = new enorm<sq_T>( rvn ); 600 tmp->set_parameters ( mu ( irvn ), Rn ); 601 return tmp; 602 } 603 604 template<class sq_T> 605 mlnorm<sq_T>* enorm<sq_T>::condition ( const RV &rvn ) { 606 607 RV rvc = rv.subt ( rvn ); 608 it_assert_debug ( ( rvc.count() +rvn.count() ==rv.count() ),"wrong rvn" ); 609 //Permutation vector of the new R 610 ivec irvn = rvn.dataind ( rv ); 611 ivec irvc = rvc.dataind ( rv ); 612 ivec perm=concat ( irvn , irvc ); 613 sq_T Rn ( R,perm ); 614 615 //fixme - could this be done in general for all sq_T? 616 mat S=R.to_mat(); 617 //fixme 618 int n=rvn.count()-1; 619 int end=R.rows()-1; 620 mat S11 = S.get ( 0,n, 0, n ); 621 mat S12 = S.get ( rvn.count(), end, 0, n ); 622 mat S22 = S.get ( rvn.count(), end, rvn.count(), end ); 623 624 vec mu1 = mu ( irvn ); 625 vec mu2 = mu ( irvc ); 626 mat A=S12*inv ( S22 ); 627 sq_T R_n ( S11 - A *S12.T() ); 628 629 mlnorm<sq_T>* tmp=new mlnorm<sq_T> ( rvn,rvc ); 630 631 tmp->set_parameters ( A,mu1-A*mu2,R_n ); 632 return tmp; 633 } 634 586 635 /////////// 587 636