00001
00013 #ifndef EF_H
00014 #define EF_H
00015
00016
00017 #include "libBM.h"
00018 #include "../math/chmat.h"
00019
00020
00021 namespace bdm {
00022
00023
00025 extern Uniform_RNG UniRNG;
00027 extern Normal_RNG NorRNG;
00029 extern Gamma_RNG GamRNG;
00030
00037 class eEF : public epdf {
00038 public:
00039
00041 eEF ( ) :epdf ( ) {};
00043 virtual double lognc() const =0;
00045 virtual void dupdate ( mat &v ) {it_error ( "Not implemented" );};
00047 virtual double evallog_nn ( const vec &val ) const{it_error ( "Not implemented" );return 0.0;};
00049 virtual double evallog ( const vec &val ) const {double tmp;tmp= evallog_nn ( val )-lognc();it_assert_debug ( std::isfinite ( tmp ),"Infinite value" ); return tmp;}
00051 virtual vec evallog ( const mat &Val ) const {
00052 vec x ( Val.cols() );
00053 for ( int i=0;i<Val.cols();i++ ) {x ( i ) =evallog_nn ( Val.get_col ( i ) ) ;}
00054 return x-lognc();
00055 }
00057 virtual void pow ( double p ) {it_error ( "Not implemented" );};
00058 };
00059
00066 class mEF : public mpdf {
00067
00068 public:
00070 mEF ( ) :mpdf ( ) {};
00071 };
00072
00074 class BMEF : public BM {
00075 protected:
00077 double frg;
00079 double last_lognc;
00080 public:
00082 BMEF ( double frg0=1.0 ) :BM ( ), frg ( frg0 ) {}
00084 BMEF ( const BMEF &B ) :BM ( B ), frg ( B.frg ), last_lognc ( B.last_lognc ) {}
00086 virtual void set_statistics ( const BMEF* BM0 ) {it_error ( "Not implemented" );};
00088 virtual void bayes ( const vec &data, const double w ) {};
00089
00090 void bayes ( const vec &dt );
00092 virtual void flatten ( const BMEF * B ) {it_error ( "Not implemented" );}
00094
00095
00096 BMEF* _copy_ ( bool changerv=false ) {it_error ( "function _copy_ not implemented for this BM" ); return NULL;};
00097 };
00098
00099 template<class sq_T>
00100 class mlnorm;
00101
00107 template<class sq_T>
00108 class enorm : public eEF {
00109 protected:
00111 vec mu;
00113 sq_T R;
00114 public:
00117
00118 enorm ( );
00119 enorm (const vec &mu,const sq_T &R ){set_parameters(mu,R);}
00120 void set_parameters ( const vec &mu,const sq_T &R );
00122
00125
00127 void dupdate ( mat &v,double nu=1.0 );
00128
00129 vec sample() const;
00130 mat sample ( int N ) const;
00131 double evallog_nn ( const vec &val ) const;
00132 double lognc () const;
00133 vec mean() const {return mu;}
00134 vec variance() const {return diag ( R.to_mat() );}
00135
00136 mpdf* condition ( const RV &rvn ) const ;
00137
00138 epdf* marginal ( const RV &rv ) const;
00140
00143
00144 vec& _mu() {return mu;}
00145 void set_mu ( const vec mu0 ) { mu=mu0;}
00146 sq_T& _R() {return R;}
00147 const sq_T& _R() const {return R;}
00149
00150 };
00151
00158 class egiw : public eEF {
00159 protected:
00161 ldmat V;
00163 double nu;
00165 int dimx;
00167 int nPsi;
00168 public:
00171 egiw() :eEF() {};
00172 egiw(int dimx0, ldmat V0, double nu0=-1.0 ) :eEF() {set_parameters(dimx0,V0, nu0);};
00173
00174 void set_parameters (int dimx0, ldmat V0, double nu0=-1.0 ) {
00175 dimx=dimx0;
00176 dim = dimx*(dimx+nPsi);
00177 nPsi = V.rows()-dimx;
00178
00179 V=V0;
00180 if ( nu0<0 ) {
00181 nu = 0.1 +nPsi +2*xdim +2;
00182
00183 }
00184 else {
00185 nu=nu0;
00186 }
00187 }
00189
00190 vec sample() const;
00191 vec mean() const;
00192 vec variance() const;
00193 void mean_mat ( mat &M, mat&R ) const;
00195 double evallog_nn ( const vec &val ) const;
00196 double lognc () const;
00197 void pow ( double p ) {V*=p;nu*=p;};
00198
00201
00202 ldmat& _V() {return V;}
00203 const ldmat& _V() const {return V;}
00204 double& _nu() {return nu;}
00205 const double& _nu() const {return nu;}
00207 };
00208
00217 class eDirich: public eEF {
00218 protected:
00220 vec beta;
00222 double gamma;
00223 public:
00225 eDirich () : eEF ( ) {};
00227 eDirich ( const eDirich &D0 ) : eEF () {set_parameters ( D0.beta );};
00228 vec sample() const {it_error ( "Not implemented" );return vec_1 ( 0.0 );};
00229 vec mean() const {return beta/gamma;};
00230 vec variance() const {return elem_mult ( beta, ( beta+1 ) ) / ( gamma* ( gamma+1 ) );}
00232 double evallog_nn ( const vec &val ) const {
00233 double tmp; tmp= ( beta-1 ) *log ( val ); it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
00234 return tmp;
00235 };
00236 double lognc () const {
00237 double tmp;
00238 double gam=sum ( beta );
00239 double lgb=0.0;
00240 for ( int i=0;i<beta.length();i++ ) {lgb+=lgamma ( beta ( i ) );}
00241 tmp= lgb-lgamma ( gam );
00242 it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
00243 return tmp;
00244 };
00246 vec& _beta() {return beta;}
00248 void set_parameters ( const vec &beta0 ) {
00249 beta= beta0;
00250 dim = beta.length();
00251 gamma = sum ( beta );
00252 }
00253 };
00254
00256 class multiBM : public BMEF {
00257 protected:
00259 eDirich est;
00261 vec β
00262 public:
00264 multiBM ( ) : BMEF ( ),est ( ),beta ( est._beta() ) {if ( beta.length() >0 ) {last_lognc=est.lognc();}else{last_lognc=0.0;}}
00266 multiBM ( const multiBM &B ) : BMEF ( B ),est ( B.est ),beta ( est._beta() ) {}
00268 void set_statistics ( const BM* mB0 ) {const multiBM* mB=dynamic_cast<const multiBM*> ( mB0 ); beta=mB->beta;}
00269 void bayes ( const vec &dt ) {
00270 if ( frg<1.0 ) {beta*=frg;last_lognc=est.lognc();}
00271 beta+=dt;
00272 if ( evalll ) {ll=est.lognc()-last_lognc;}
00273 }
00274 double logpred ( const vec &dt ) const {
00275 eDirich pred ( est );
00276 vec &beta = pred._beta();
00277
00278 double lll;
00279 if ( frg<1.0 )
00280 {beta*=frg;lll=pred.lognc();}
00281 else
00282 if ( evalll ) {lll=last_lognc;}
00283 else{lll=pred.lognc();}
00284
00285 beta+=dt;
00286 return pred.lognc()-lll;
00287 }
00288 void flatten ( const BMEF* B ) {
00289 const multiBM* E=dynamic_cast<const multiBM*> ( B );
00290
00291 const vec &Eb=E->beta;
00292 beta*= ( sum ( Eb ) /sum ( beta ) );
00293 if ( evalll ) {last_lognc=est.lognc();}
00294 }
00295 const epdf& _epdf() const {return est;};
00296 const eDirich* _e() const {return &est;};
00297 void set_parameters ( const vec &beta0 ) {
00298 est.set_parameters ( beta0 );
00299 if ( evalll ) {last_lognc=est.lognc();}
00300 }
00301 };
00302
00312 class egamma : public eEF {
00313 protected:
00315 vec alpha;
00317 vec beta;
00318 public :
00320 egamma ( ) :eEF ( ), alpha(), beta() {};
00322 void set_parameters ( const vec &a, const vec &b ) {alpha=a,beta=b;dim = alpha.length();};
00323 vec sample() const;
00325
00326 double evallog ( const vec &val ) const;
00327 double lognc () const;
00329 vec& _alpha() {return alpha;}
00330 vec& _beta() {return beta;}
00331 vec mean() const {return elem_div ( alpha,beta );}
00332 vec variance() const {return elem_div ( alpha,elem_mult ( beta,beta ) ); }
00333 };
00334
00349 class eigamma : public eEF {
00350 protected:
00352 egamma eg;
00354 vec α
00356 vec β
00357 public :
00359 eigamma ( ) :eEF ( ), eg(),alpha ( eg._alpha() ), beta ( eg._beta() ) {};
00361 void set_parameters ( const vec &a, const vec &b ) {alpha=a,beta=b;};
00362 vec sample() const {return 1.0/eg.sample();};
00364
00365 double evallog ( const vec &val ) const {return eg.evallog ( val );};
00366 double lognc () const {return eg.lognc();};
00368 vec mean() const {return elem_div ( beta,alpha-1 );}
00369 vec variance() const {vec mea=mean(); return elem_div ( elem_mult ( mea,mea ),alpha-2 );}
00370 vec& _alpha() {return alpha;}
00371 vec& _beta() {return beta;}
00372 };
00373
00375
00376
00377
00378
00379
00380
00382
00383
00384
00385
00386
00387
00388
00390
00391 class euni: public epdf {
00392 protected:
00394 vec low;
00396 vec high;
00398 vec distance;
00400 double nk;
00402 double lnk;
00403 public:
00405 euni ( ) :epdf ( ) {}
00406 double eval ( const vec &val ) const {return nk;}
00407 double evallog ( const vec &val ) const {return lnk;}
00408 vec sample() const {
00409 vec smp ( dim );
00410 #pragma omp critical
00411 UniRNG.sample_vector ( dim ,smp );
00412 return low+elem_mult ( distance,smp );
00413 }
00415 void set_parameters ( const vec &low0, const vec &high0 ) {
00416 distance = high0-low0;
00417 it_assert_debug ( min ( distance ) >0.0,"bad support" );
00418 low = low0;
00419 high = high0;
00420 nk = prod ( 1.0/distance );
00421 lnk = log ( nk );
00422 dim = low.length();
00423 }
00424 vec mean() const {return ( high-low ) /2.0;}
00425 vec variance() const {return ( pow ( high,2 ) +pow ( low,2 ) +elem_mult ( high,low ) ) /3.0;}
00426 };
00427
00428
00434 template<class sq_T>
00435 class mlnorm : public mEF {
00436 protected:
00438 enorm<sq_T> epdf;
00439 mat A;
00440 vec mu_const;
00441 vec& _mu;
00442 public:
00444 mlnorm ( );
00446 void set_parameters ( const mat &A, const vec &mu0, const sq_T &R );
00447
00448
00449
00450
00452 void condition ( const vec &cond );
00453
00455 vec& _mu_const() {return mu_const;}
00457 mat& _A() {return A;}
00459 mat _R() {return epdf._R().to_mat();}
00460
00461 template<class sq_M>
00462 friend std::ostream &operator<< ( std::ostream &os, mlnorm<sq_M> &ml );
00463 };
00464
00472 class mlstudent : public mlnorm<ldmat> {
00473 protected:
00474 ldmat Lambda;
00475 ldmat &_R;
00476 ldmat Re;
00477 public:
00478 mlstudent ( ) :mlnorm<ldmat> (),
00479 Lambda (), _R ( epdf._R() ) {}
00480 void set_parameters ( const mat &A0, const vec &mu0, const ldmat &R0, const ldmat& Lambda0 ) {
00481 it_assert_debug ( A0.rows() ==mu0.length(),"" );
00482 it_assert_debug ( R0.rows() ==A0.rows(),"" );
00483
00484 epdf.set_parameters ( mu0,Lambda );
00485 A = A0;
00486 mu_const = mu0;
00487 Re=R0;
00488 Lambda = Lambda0;
00489 }
00490 void condition ( const vec &cond ) {
00491 _mu = A*cond + mu_const;
00492 double zeta;
00493
00494 if ( ( cond.length() +1 ) ==Lambda.rows() ) {
00495 zeta = Lambda.invqform ( concat ( cond, vec_1 ( 1.0 ) ) );
00496 }
00497 else {
00498 zeta = Lambda.invqform ( cond );
00499 }
00500 _R = Re;
00501 _R*= ( 1+zeta );
00502 };
00503
00504 };
00514 class mgamma : public mEF {
00515 protected:
00517 egamma epdf;
00519 double k;
00521 vec &_beta;
00522
00523 public:
00525 mgamma ( ) : mEF ( ), epdf (), _beta ( epdf._beta() ) {ep=&epdf;};
00527 void set_parameters ( double k );
00528 void condition ( const vec &val ) {_beta=k/val;};
00529 };
00530
00540 class migamma : public mEF {
00541 protected:
00543 eigamma epdf;
00545 double k;
00547 vec &_alpha;
00549 vec &_beta;
00550
00551 public:
00554 migamma ( ) : mEF (), epdf ( ), _alpha ( epdf._alpha() ), _beta ( epdf._beta() ) {ep=&epdf;};
00555 migamma ( const migamma &m ) : mEF (), epdf ( m.epdf ), _alpha ( epdf._alpha() ), _beta ( epdf._beta() ) {ep=&epdf;};
00557
00559 void set_parameters ( int len, double k0 ) {
00560 k=k0;
00561 epdf.set_parameters ( ( 1.0/ ( k*k ) +2.0 ) *ones ( len ) , ones ( len ) );
00562 };
00563 void condition ( const vec &val ) {
00564 _beta=elem_mult ( val, ( _alpha-1.0 ) );
00565 };
00566 };
00567
00579 class mgamma_fix : public mgamma {
00580 protected:
00582 double l;
00584 vec refl;
00585 public:
00587 mgamma_fix ( ) : mgamma ( ),refl () {};
00589 void set_parameters ( double k0 , vec ref0, double l0 ) {
00590 mgamma::set_parameters ( k0 );
00591 refl=pow ( ref0,1.0-l0 );l=l0;
00592 };
00593
00594 void condition ( const vec &val ) {vec mean=elem_mult ( refl,pow ( val,l ) ); _beta=k/mean;};
00595 };
00596
00597
00610 class migamma_fix : public migamma {
00611 protected:
00613 double l;
00615 vec refl;
00616 public:
00618 migamma_fix ( ) : migamma (),refl ( ) {};
00620 void set_parameters ( double k0 , vec ref0, double l0 ) {
00621 migamma::set_parameters (ref0.length(), k0 );
00622 refl=pow ( ref0,1.0-l0 );
00623 l=l0;
00624 };
00625
00626 void condition ( const vec &val ) {
00627 vec mean=elem_mult ( refl,pow ( val,l ) );
00628 migamma::condition ( mean );
00629 };
00630 };
00632 enum RESAMPLING_METHOD { MULTINOMIAL = 0, STRATIFIED = 1, SYSTEMATIC = 3 };
00638 class eEmp: public epdf {
00639 protected :
00641 int n;
00643 vec w;
00645 Array<vec> samples;
00646 public:
00648 eEmp ( int n0 ) :epdf ( ),n ( n0 ),w ( n ),samples ( n ) {};
00650 void set_parameters ( const vec &w0, const epdf* pdf0 );
00652 void set_samples ( const epdf* pdf0 );
00654 void set_n ( int n0, bool copy=true ) {w.set_size ( n0,copy );samples.set_size ( n0,copy );};
00656 vec& _w() {return w;};
00658 const vec& _w() const {return w;};
00660 Array<vec>& _samples() {return samples;};
00662 const Array<vec>& _samples() const {return samples;};
00664 ivec resample ( RESAMPLING_METHOD method = SYSTEMATIC );
00666 vec sample() const {it_error ( "Not implemented" );return 0;}
00668 double evallog ( const vec &val ) const {it_error ( "Not implemented" );return 0.0;}
00669 vec mean() const {
00670 vec pom=zeros ( dim );
00671 for ( int i=0;i<n;i++ ) {pom+=samples ( i ) *w ( i );}
00672 return pom;
00673 }
00674 vec variance() const {
00675 vec pom=zeros ( dim );
00676 for ( int i=0;i<n;i++ ) {pom+=pow ( samples ( i ),2 ) *w ( i );}
00677 return pom-pow ( mean(),2 );
00678 }
00679 };
00680
00681
00683
00684 template<class sq_T>
00685 enorm<sq_T>::enorm ( ) :eEF ( ), mu ( ),R ( ) {};
00686
00687 template<class sq_T>
00688 void enorm<sq_T>::set_parameters ( const vec &mu0, const sq_T &R0 ) {
00689
00690 mu = mu0;
00691 R = R0;
00692 dim = mu0.length();
00693 };
00694
00695 template<class sq_T>
00696 void enorm<sq_T>::dupdate ( mat &v, double nu ) {
00697
00698 };
00699
00700
00701
00702
00703
00704
00705 template<class sq_T>
00706 vec enorm<sq_T>::sample() const {
00707 vec x ( dim );
00708 #pragma omp critical
00709 NorRNG.sample_vector ( dim,x );
00710 vec smp = R.sqrt_mult ( x );
00711
00712 smp += mu;
00713 return smp;
00714 };
00715
00716 template<class sq_T>
00717 mat enorm<sq_T>::sample ( int N ) const {
00718 mat X ( dim,N );
00719 vec x ( dim );
00720 vec pom;
00721 int i;
00722
00723 for ( i=0;i<N;i++ ) {
00724 #pragma omp critical
00725 NorRNG.sample_vector ( dim,x );
00726 pom = R.sqrt_mult ( x );
00727 pom +=mu;
00728 X.set_col ( i, pom );
00729 }
00730
00731 return X;
00732 };
00733
00734
00735
00736
00737
00738
00739
00740
00741
00742 template<class sq_T>
00743 double enorm<sq_T>::evallog_nn ( const vec &val ) const {
00744
00745 double tmp=-0.5* ( R.invqform ( mu-val ) );
00746 return tmp;
00747 };
00748
00749 template<class sq_T>
00750 inline double enorm<sq_T>::lognc () const {
00751
00752 double tmp=0.5* ( R.cols() * 1.83787706640935 +R.logdet() );
00753 return tmp;
00754 };
00755
00756 template<class sq_T>
00757 mlnorm<sq_T>::mlnorm ( ) :mEF (),epdf ( ),A ( ),_mu ( epdf._mu() ) {
00758 ep =&epdf;
00759 }
00760
00761 template<class sq_T>
00762 void mlnorm<sq_T>::set_parameters ( const mat &A0, const vec &mu0, const sq_T &R0 ) {
00763 it_assert_debug ( A0.rows() ==mu0.length(),"" );
00764 it_assert_debug ( A0.rows() ==R0.rows(),"" );
00765
00766 epdf.set_parameters ( zeros ( A.rows() ),R0 );
00767 A = A0;
00768 mu_const = mu0;
00769 }
00770
00771
00772
00773
00774
00775
00776
00777
00778
00779
00780
00781
00782
00783
00784
00785
00786
00787
00788
00789
00790
00791
00792
00793
00794
00795
00796 template<class sq_T>
00797 void mlnorm<sq_T>::condition ( const vec &cond ) {
00798 _mu = A*cond + mu_const;
00799
00800 }
00801
00802 template<class sq_T>
00803 epdf* enorm<sq_T>::marginal ( const RV &rvn ) const {
00804 ivec irvn = rvn.dataind ( rv );
00805
00806 sq_T Rn ( R,irvn );
00807 enorm<sq_T>* tmp = new enorm<sq_T> ( );
00808 tmp->set_rv ( rvn );
00809 tmp->set_parameters ( mu ( irvn ), Rn );
00810 return tmp;
00811 }
00812
00813 template<class sq_T>
00814 mpdf* enorm<sq_T>::condition ( const RV &rvn ) const {
00815
00816 it_assert_debug ( isnamed(),"" );
00817
00818 RV rvc = rv.subt ( rvn );
00819 it_assert_debug ( ( rvc._dsize() +rvn._dsize() ==rv._dsize() ),"wrong rvn" );
00820
00821 ivec irvn = rvn.dataind ( rv );
00822 ivec irvc = rvc.dataind ( rv );
00823 ivec perm=concat ( irvn , irvc );
00824 sq_T Rn ( R,perm );
00825
00826
00827 mat S=Rn.to_mat();
00828
00829 int n=rvn._dsize()-1;
00830 int end=R.rows()-1;
00831 mat S11 = S.get ( 0,n, 0, n );
00832 mat S12 = S.get ( 0, n , rvn._dsize(), end );
00833 mat S22 = S.get ( rvn._dsize(), end, rvn._dsize(), end );
00834
00835 vec mu1 = mu ( irvn );
00836 vec mu2 = mu ( irvc );
00837 mat A=S12*inv ( S22 );
00838 sq_T R_n ( S11 - A *S12.T() );
00839
00840 mlnorm<sq_T>* tmp=new mlnorm<sq_T> ( );
00841
00842 tmp->set_parameters ( A,mu1-A*mu2,R_n );
00843 return tmp;
00844 }
00845
00847
00848 template<class sq_T>
00849 std::ostream &operator<< ( std::ostream &os, mlnorm<sq_T> &ml ) {
00850 os << "A:"<< ml.A<<endl;
00851 os << "mu:"<< ml.mu_const<<endl;
00852 os << "R:" << ml.epdf._R().to_mat() <<endl;
00853 return os;
00854 };
00855
00856 }
00857 #endif //EF_H