- Main Page
- Related Pages
- Classes
- Files

# work/git/mixpp/bdm/stat/libEF.h

Go to the documentation of this file.

    00001
    00013 #ifndef EF_H
    00014 #define EF_H
    00015
    00016 #include <itpp/itbase.h>
    00017 #include "../math/libDC.h"
    00018 #include "libBM.h"
    00019 #include "../itpp_ext.h"
    00020 //#include <std>
    00021
    00022 using namespace itpp;
    00023
    00024
    00026 extern Uniform_RNG UniRNG;
    00028 extern Normal_RNG NorRNG;
    00030 extern Gamma_RNG GamRNG;
    00031
    00038 class eEF : public epdf {
    00039 public:
    00040 //      eEF() :epdf() {};
    00042         eEF ( const RV &rv ) :epdf ( rv ) {};
    00044         virtual double lognc() const =0;
    00046         virtual void dupdate ( mat &v ) {it_error ( "Not implemneted" );};
    00048         virtual double evalpdflog_nn ( const vec &val ) const{it_error ( "Not implemneted" );return 0.0;};
    00050         virtual double evalpdflog ( const vec &val ) const {return evalpdflog_nn ( val )-lognc();}
    00052         virtual vec evalpdflog ( const mat &Val ) const {
    00053                 vec x ( Val.cols() );
    00054                 for ( int i=0;i<Val.cols();i++ ) {x ( i ) =evalpdflog_nn ( Val.get_col ( i ) ) ;}
    00055                 return x-lognc();
    00056         }
    00058         virtual void pow ( double p ) {it_error ( "Not implemented" );};
    00059 };
    00060
    00067 class mEF : public mpdf {
    00068
    00069 public:
    00071         mEF ( const RV &rv0, const RV &rvc0 ) :mpdf ( rv0,rvc0 ) {};
    00072 };
    00073
    00075 class BMEF : public BM {
    00076 protected:
    00078         double frg;
    00080         double last_lognc;
    00081 public:
    00083         BMEF ( const RV &rv, double frg0=1.0 ) :BM ( rv ), frg ( frg0 ) {}
    00085         BMEF ( const BMEF &B ) :BM ( B ), frg ( B.frg ), last_lognc ( B.last_lognc ) {}
    00087         virtual void set_statistics ( const BMEF* BM0 ) {it_error ( "Not implemented" );};
    00089         virtual void bayes ( const vec &data, const double w ) {};
    00090         //original Bayes
    00091         void bayes ( const vec &dt );
    00093         virtual void flatten ( BMEF * B) {it_error ( "Not implemented" );}
    00094 };
    00095
    00101 template<class sq_T>
    00102
    00103 class enorm : public eEF {
    00104 protected:
    00106         vec mu;
    00108         sq_T R;
    00110         int dim;
    00111 public:
    00112 //      enorm() :eEF() {};
    00114         enorm ( RV &rv );
    00116         void set_parameters ( const vec &mu,const sq_T &R );
    00118         void tupdate ( double phi, mat &vbar, double nubar );
    00120         void dupdate ( mat &v,double nu=1.0 );
    00121
    00122         vec sample() const;
    00124         mat sample ( int N ) const;
    00125         double eval ( const vec &val ) const ;
    00126         double evalpdflog ( const vec &val ) const;
    00127         double lognc () const;
    00128         vec mean() const {return mu;}
    00129
    00130 //Access methods
    00132         vec& _mu() {return mu;}
    00133
    00135         void set_mu ( const vec mu0 ) { mu=mu0;}
    00136
    00138         sq_T& _R() {return R;}
    00139
    00141         mat getR () {return R.to_mat();}
    00142 };
    00143
    00150 class egiw : public eEF {
    00151 protected:
    00153         ldmat V;
    00155         double nu;
    00157         int xdim;
    00159         int nPsi;
    00160 public:
    00162         egiw ( RV rv, mat V0, double nu0 ) : eEF ( rv ), V ( V0 ), nu ( nu0 ) {
    00163                 xdim = rv.count() /V.rows();
    00164                 it_assert_debug ( rv.count() ==xdim*V.rows(),"Incompatible V0." );
    00165                 nPsi = V.rows()-xdim;
    00166         }
    00168         egiw ( RV rv, ldmat V0, double nu0 ) : eEF ( rv ), V ( V0 ), nu ( nu0 ) {
    00169                 xdim = rv.count() /V.rows();
    00170                 it_assert_debug ( rv.count() ==xdim*V.rows(),"Incompatible V0." );
    00171                 nPsi = V.rows()-xdim;
    00172         }
    00173
    00174         vec sample() const;
    00175         vec mean() const;
    00176         void mean_mat ( mat &M, mat&R ) const;
    00178         double evalpdflog_nn ( const vec &val ) const;
    00179         double lognc () const;
    00180
    00181         //Access
    00183         ldmat& _V() {return V;}
    00185         double& _nu() {return nu;}
    00186         void pow ( double p );
    00187 };
    00188
    00197 class eDirich: public eEF {
    00198 protected:
    00200         vec beta;
    00201 public:
    00203         eDirich ( const RV &rv, const vec &beta0 ) : eEF ( rv ),beta ( beta0 ) {it_assert_debug ( rv.count() ==beta.length(),"Incompatible statistics" ); };
    00205         eDirich ( const eDirich &D0 ) : eEF ( D0.rv ),beta ( D0.beta ) {};
    00206         vec sample() const {it_error ( "Not implemented" );return vec_1 ( 0.0 );};
    00207         vec mean() const {return beta/sum ( beta );};
    00209         double evalpdflog_nn ( const vec &val ) const {return ( beta-1 ) *log ( val );};
    00210         double lognc () const {
    00211                 double gam=sum ( beta );
    00212                 double lgb=0.0;
    00213                 for ( int i=0;i<beta.length();i++ ) {lgb+=lgamma ( beta ( i ) );}
    00214                 return lgb-lgamma ( gam );
    00215         };
    00217         vec& _beta() {return beta;}
    00218 };
    00219
    00221 class multiBM : public BMEF {
    00222 protected:
    00224         eDirich est;
    00225         vec &beta;
    00226 public:
    00228         multiBM ( const RV &rv, const vec beta0 ) : BMEF ( rv ),est ( rv,beta0 ),beta ( est._beta() ) {last_lognc=est.lognc();}
    00230         multiBM ( const multiBM &B ) : BMEF ( B ),est ( rv,B.beta ),beta ( est._beta() ) {}
    00231
    00232         void set_statistics ( const BM* mB0 ) {const multiBM* mB=dynamic_cast<const multiBM*> ( mB0 ); beta=mB->beta;}
    00233         void bayes ( const vec &dt ) {
    00234                 if ( frg<1.0 ) {beta*=frg;last_lognc=est.lognc();}
    00235                 beta+=dt;
    00236                 if ( evalll ) {ll=est.lognc()-last_lognc;}
    00237         }
    00238         double logpred ( const vec &dt ) const {
    00239                 eDirich pred ( est );
    00240                 vec &beta = pred._beta();
    00241
    00242                 double lll;
    00243                 if ( frg<1.0 )
    00244                         {beta*=frg;lll=pred.lognc();}
    00245                 else
    00246                         if ( evalll ) {lll=last_lognc;}
    00247                         else{lll=pred.lognc();}
    00248
    00249                 beta+=dt;
    00250                 return pred.lognc()-lll;
    00251         }
    00252         void flatten (BMEF* B ) {
    00253                 eDirich* E=dynamic_cast<eDirich*>(B);
    00254                 // sum(beta) should be equal to sum(B.beta)
    00255                 const vec &Eb=E->_beta();
    00256                 est.pow ( sum(beta)/sum(Eb) );
    00257                 if(evalll){last_lognc=est.lognc();}
    00258         }
    00259         const epdf& _epdf() const {return est;};
    00261 };
    00262
    00272 class egamma : public eEF {
    00273 protected:
    00275         vec alpha;
    00277         vec beta;
    00278 public :
    00280         egamma ( const RV &rv ) :eEF ( rv ) {};
    00282         void set_parameters ( const vec &a, const vec &b ) {alpha=a,beta=b;};
    00283         vec sample() const;
    00285 //      mat sample ( int N ) const;
    00286         double evalpdflog ( const vec &val ) const;
    00287         double lognc () const;
    00289         void _param ( vec* &a, vec* &b ) {a=&alpha;b=&beta;};
    00290         vec mean() const {vec pom ( alpha ); pom/=beta; return pom;}
    00291 };
    00292 /*
    00294 class emix : public epdf {
    00295 protected:
    00296         int n;
    00297         vec &w;
    00298         Array<epdf*> Coms;
    00299 public:
    00301         emix ( const RV &rv, vec &w0): epdf(rv), n(w0.length()), w(w0), Coms(n) {};
    00302         void set_parameters( int &i, double wi, epdf* ep){w(i)=wi;Coms(i)=ep;}
    00303         vec mean(){vec pom; for(int i=0;i<n;i++){pom+=Coms(i)->mean()*w(i);} return pom;};
    00304         vec sample() {it_error ( "Not implemented" );return 0;}
    00305 };
    00306 */
    00307
    00309
    00310 class euni: public epdf {
    00311 protected:
    00313         vec low;
    00315         vec high;
    00317         vec distance;
    00319         double nk;
    00321         double lnk;
    00322 public:
    00324         euni ( const RV rv ) :epdf ( rv ) {}
    00325         double eval ( const vec &val ) const  {return nk;}
    00326         double evalpdflog ( const vec &val ) const  {return lnk;}
    00327         vec sample() const {
    00328                 vec smp ( rv.count() );
    00329 #pragma omp critical
    00330                 UniRNG.sample_vector ( rv.count(),smp );
    00331                 return low+elem_mult ( distance,smp );
    00332         }
    00334         void set_parameters ( const vec &low0, const vec &high0 ) {
    00335                 distance = high0-low0;
    00336                 it_assert_debug ( min ( distance ) >0.0,"bad support" );
    00337                 low = low0;
    00338                 high = high0;
    00339                 nk = prod ( 1.0/distance );
    00340                 lnk = log ( nk );
    00341         }
    00342         vec mean() const {vec pom=high; pom-=low; pom/=2.0; return pom;}
    00343 };
    00344
    00345
    00351 template<class sq_T>
    00352 class mlnorm : public mEF {
    00354         enorm<sq_T> epdf;
    00355         mat A;
    00356         vec& _mu; //cached epdf.mu;
    00357 public:
    00359         mlnorm ( RV &rv,RV &rvc );
    00361         void set_parameters ( const  mat &A, const sq_T &R );
    00363         vec samplecond ( vec &cond, double &lik );
    00365         mat samplecond ( vec &cond, vec &lik, int n );
    00367         void condition ( vec &cond );
    00368 };
    00369
    00379 class mgamma : public mEF {
    00380 protected:
    00382         egamma epdf;
    00384         double k;
    00386         vec* _beta;
    00387
    00388 public:
    00390         mgamma ( const RV &rv,const RV &rvc );
    00392         void set_parameters ( double k );
    00393         void condition ( const vec &val ) {*_beta=k/val;};
    00394 };
    00395
    00407 class mgamma_fix : public mgamma {
    00408 protected:
    00410         double l;
    00412         vec refl;
    00413 public:
    00415         mgamma_fix ( const RV &rv,const RV &rvc ) : mgamma ( rv,rvc ),refl ( rv.count() ) {};
    00417         void set_parameters ( double k0 , vec ref0, double l0 ) {
    00418                 mgamma::set_parameters ( k0 );
    00419                 refl=pow ( ref0,1.0-l0 );l=l0;
    00420         };
    00421
    00422         void condition ( const vec &val ) {vec mean=elem_mult ( refl,pow ( val,l ) ); *_beta=k/mean;};
    00423 };
    00424
    00426 enum RESAMPLING_METHOD { MULTINOMIAL = 0, STRATIFIED = 1, SYSTEMATIC = 3 };
    00432 class eEmp: public epdf {
    00433 protected :
    00435         int n;
    00437         vec w;
    00439         Array<vec> samples;
    00440 public:
    00442         eEmp ( const RV &rv0 ,int n0 ) :epdf ( rv0 ),n ( n0 ),w ( n ),samples ( n ) {};
    00444         void set_parameters ( const vec &w0, epdf* pdf0 );
    00446         vec& _w()  {return w;};
    00448         Array<vec>& _samples() {return samples;};
    00450         ivec resample ( RESAMPLING_METHOD method = SYSTEMATIC );
    00452         vec sample() const {it_error ( "Not implemented" );return 0;}
    00454         double evalpdflog ( const vec &val ) const {it_error ( "Not implemented" );return 0.0;}
    00455         vec mean() const {
    00456                 vec pom=zeros ( rv.count() );
    00457                 for ( int i=0;i<n;i++ ) {pom+=samples ( i ) *w ( i );}
    00458                 return pom;
    00459         }
    00460 };
    00461
    00462
    00464
    00465 template<class sq_T>
    00466 enorm<sq_T>::enorm ( RV &rv ) :eEF ( rv ), mu ( rv.count() ),R ( rv.count() ),dim ( rv.count() ) {};
    00467
    00468 template<class sq_T>
    00469 void enorm<sq_T>::set_parameters ( const vec &mu0, const sq_T &R0 ) {
    00470 //Fixme test dimensions of mu0 and R0;
    00471         mu = mu0;
    00472         R = R0;
    00473 };
    00474
    00475 template<class sq_T>
    00476 void enorm<sq_T>::dupdate ( mat &v, double nu ) {
    00477         //
    00478 };
    00479
    00480 template<class sq_T>
    00481 void enorm<sq_T>::tupdate ( double phi, mat &vbar, double nubar ) {
    00482         //
    00483 };
    00484
    00485 template<class sq_T>
    00486 vec enorm<sq_T>::sample() const {
    00487         vec x ( dim );
    00488         NorRNG.sample_vector ( dim,x );
    00489         vec smp = R.sqrt_mult ( x );
    00490
    00491         smp += mu;
    00492         return smp;
    00493 };
    00494
    00495 template<class sq_T>
    00496 mat enorm<sq_T>::sample ( int N ) const {
    00497         mat X ( dim,N );
    00498         vec x ( dim );
    00499         vec pom;
    00500         int i;
    00501
    00502         for ( i=0;i<N;i++ ) {
    00503                 NorRNG.sample_vector ( dim,x );
    00504                 pom = R.sqrt_mult ( x );
    00505                 pom +=mu;
    00506                 X.set_col ( i, pom );
    00507         }
    00508
    00509         return X;
    00510 };
    00511
    00512 template<class sq_T>
    00513 double enorm<sq_T>::eval ( const vec &val ) const {
    00514         double pdfl,e;
    00515         pdfl = evalpdflog ( val );
    00516         e = exp ( pdfl );
    00517         return e;
    00518 };
    00519
    00520 template<class sq_T>
    00521 double enorm<sq_T>::evalpdflog ( const vec &val ) const {
    00522         // 1.83787706640935 = log(2pi)
    00523         return  -0.5* ( +R.invqform ( mu-val ) ) - lognc();
    00524 };
    00525
    00526 template<class sq_T>
    00527 inline double enorm<sq_T>::lognc () const {
    00528         // 1.83787706640935 = log(2pi)
    00529         return -0.5* ( R.cols() * 1.83787706640935 +R.logdet() );
    00530 };
    00531
    00532 template<class sq_T>
    00533 mlnorm<sq_T>::mlnorm ( RV &rv0,RV &rvc0 ) :mEF ( rv0,rvc0 ),epdf ( rv0 ),A ( rv0.count(),rv0.count() ),_mu ( epdf._mu() ) {
    00534         ep =&epdf;
    00535 }
    00536
    00537 template<class sq_T>
    00538 void mlnorm<sq_T>::set_parameters ( const mat &A0, const sq_T &R0 ) {
    00539         epdf.set_parameters ( zeros ( rv.count() ),R0 );
    00540         A = A0;
    00541 }
    00542
    00543 template<class sq_T>
    00544 vec mlnorm<sq_T>::samplecond ( vec &cond, double &lik ) {
    00545         this->condition ( cond );
    00546         vec smp = epdf.sample();
    00547         lik = epdf.eval ( smp );
    00548         return smp;
    00549 }
    00550
    00551 template<class sq_T>
    00552 mat mlnorm<sq_T>::samplecond ( vec &cond, vec &lik, int n ) {
    00553         int i;
    00554         int dim = rv.count();
    00555         mat Smp ( dim,n );
    00556         vec smp ( dim );
    00557         this->condition ( cond );
    00558
    00559         for ( i=0; i<n; i++ ) {
    00560                 smp = epdf.sample();
    00561                 lik ( i ) = epdf.eval ( smp );
    00562                 Smp.set_col ( i ,smp );
    00563         }
    00564
    00565         return Smp;
    00566 }
    00567
    00568 template<class sq_T>
    00569 void mlnorm<sq_T>::condition ( vec &cond ) {
    00570         _mu = A*cond;
    00571 //R is already assigned;
    00572 }
    00573
    00575
    00576
    00577 #endif //EF_H