- 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 implemented" );};
    00048         virtual double evalpdflog_nn ( const vec &val ) const{it_error ( "Not implemented" );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 ( const BMEF * B ) {it_error ( "Not implemented" );}
    00095         virtual void flatten ( double nu0 ) {it_error ( "Not implemented" );}
    00096 };
    00097
    00098 template<class sq_T>
    00099 class mlnorm;
    00100
    00106 template<class sq_T>
    00107 class enorm : public eEF {
    00108 protected:
    00110         vec mu;
    00112         sq_T R;
    00114         int dim;
    00115 public:
    00117         enorm ( const RV &rv );
    00119         void set_parameters ( const vec &mu,const sq_T &R );
    00121         //void tupdate ( double phi, mat &vbar, double nubar );
    00123         void dupdate ( mat &v,double nu=1.0 );
    00124
    00125         vec sample() const;
    00127         mat sample ( int N ) const;
    00128         double eval ( const vec &val ) const ;
    00129         double evalpdflog_nn ( const vec &val ) const;
    00130         double lognc () const;
    00131         vec mean() const {return mu;}
    00132         mlnorm<sq_T>* condition ( const RV &rvn );
    00133         enorm<sq_T>* marginal ( const RV &rv );
    00134 //Access methods
    00136         vec& _mu() {return mu;}
    00137
    00139         void set_mu ( const vec mu0 ) { mu=mu0;}
    00140
    00142         sq_T& _R() {return R;}
    00143
    00145 //      mat getR () {return R.to_mat();}
    00146 };
    00147
    00154 class egiw : public eEF {
    00155 protected:
    00157         ldmat V;
    00159         double nu;
    00161         int xdim;
    00163         int nPsi;
    00164 public:
    00166         egiw ( RV rv, mat V0, double nu0 ) : eEF ( rv ), V ( V0 ), nu ( nu0 ) {
    00167                 xdim = rv.count() /V.rows();
    00168                 it_assert_debug ( rv.count() ==xdim*V.rows(),"Incompatible V0." );
    00169                 nPsi = V.rows()-xdim;
    00170         }
    00172         egiw ( RV rv, ldmat V0, double nu0 ) : eEF ( rv ), V ( V0 ), nu ( nu0 ) {
    00173                 xdim = rv.count() /V.rows();
    00174                 it_assert_debug ( rv.count() ==xdim*V.rows(),"Incompatible V0." );
    00175                 nPsi = V.rows()-xdim;
    00176         }
    00177
    00178         vec sample() const;
    00179         vec mean() const;
    00180         void mean_mat ( mat &M, mat&R ) const;
    00182         double evalpdflog_nn ( const vec &val ) const;
    00183         double lognc () const;
    00184
    00185         //Access
    00187         ldmat& _V() {return V;}
    00189         double& _nu() {return nu;}
    00190         void pow ( double p );
    00191 };
    00192
    00201 class eDirich: public eEF {
    00202 protected:
    00204         vec beta;
    00205 public:
    00207         eDirich ( const RV &rv, const vec &beta0 ) : eEF ( rv ),beta ( beta0 ) {it_assert_debug ( rv.count() ==beta.length(),"Incompatible statistics" ); };
    00209         eDirich ( const eDirich &D0 ) : eEF ( D0.rv ),beta ( D0.beta ) {};
    00210         vec sample() const {it_error ( "Not implemented" );return vec_1 ( 0.0 );};
    00211         vec mean() const {return beta/sum ( beta );};
    00213         double evalpdflog_nn ( const vec &val ) const {return ( beta-1 ) *log ( val );};
    00214         double lognc () const {
    00215                 double gam=sum ( beta );
    00216                 double lgb=0.0;
    00217                 for ( int i=0;i<beta.length();i++ ) {lgb+=lgamma ( beta ( i ) );}
    00218                 return lgb-lgamma ( gam );
    00219         };
    00221         vec& _beta()  {return beta;}
    00223         void set_parameters ( const vec &beta0 ) {
    00224                 if ( beta0.length() !=beta.length() ) {
    00225                         it_assert_debug ( rv.length() ==1,"Undefined" );
    00226                         rv.set_size ( 0,beta0.length() );
    00227                 }
    00228                 beta= beta0;
    00229         }
    00230 };
    00231
    00233 class multiBM : public BMEF {
    00234 protected:
    00236         eDirich est;
    00238         vec &beta;
    00239 public:
    00241         multiBM ( const RV &rv, const vec beta0 ) : BMEF ( rv ),est ( rv,beta0 ),beta ( est._beta() ) {last_lognc=est.lognc();}
    00243         multiBM ( const multiBM &B ) : BMEF ( B ),est ( rv,B.beta ),beta ( est._beta() ) {}
    00245         void set_statistics ( const BM* mB0 ) {const multiBM* mB=dynamic_cast<const multiBM*> ( mB0 ); beta=mB->beta;}
    00246         void bayes ( const vec &dt ) {
    00247                 if ( frg<1.0 ) {beta*=frg;last_lognc=est.lognc();}
    00248                 beta+=dt;
    00249                 if ( evalll ) {ll=est.lognc()-last_lognc;}
    00250         }
    00251         double logpred ( const vec &dt ) const {
    00252                 eDirich pred ( est );
    00253                 vec &beta = pred._beta();
    00254
    00255                 double lll;
    00256                 if ( frg<1.0 )
    00257                         {beta*=frg;lll=pred.lognc();}
    00258                 else
    00259                         if ( evalll ) {lll=last_lognc;}
    00260                         else{lll=pred.lognc();}
    00261
    00262                 beta+=dt;
    00263                 return pred.lognc()-lll;
    00264         }
    00265         void flatten ( const BMEF* B ) {
    00266                 const eDirich* E=dynamic_cast<const eDirich*> ( B );
    00267                 // sum(beta) should be equal to sum(B.beta)
    00268                 const vec &Eb=const_cast<eDirich*> ( E )->_beta();
    00269                 est.pow ( sum ( beta ) /sum ( Eb ) );
    00270                 if ( evalll ) {last_lognc=est.lognc();}
    00271         }
    00272         const epdf& _epdf() const {return est;};
    00273         void set_parameters ( const vec &beta0 ) {
    00274                 est.set_parameters ( beta0 );
    00275                 rv = est._rv();
    00276                 if ( evalll ) {last_lognc=est.lognc();}
    00277         }
    00278 };
    00279
    00289 class egamma : public eEF {
    00290 protected:
    00292         vec alpha;
    00294         vec beta;
    00295 public :
    00297         egamma ( const RV &rv ) :eEF ( rv ) {};
    00299         void set_parameters ( const vec &a, const vec &b ) {alpha=a,beta=b;};
    00300         vec sample() const;
    00302 //      mat sample ( int N ) const;
    00303         double evalpdflog ( const vec &val ) const;
    00304         double lognc () const;
    00306         void _param ( vec* &a, vec* &b ) {a=&alpha;b=&beta;};
    00307         vec mean() const {vec pom ( alpha ); pom/=beta; return pom;}
    00308 };
    00309 /*
    00311 class emix : public epdf {
    00312 protected:
    00313         int n;
    00314         vec &w;
    00315         Array<epdf*> Coms;
    00316 public:
    00318         emix ( const RV &rv, vec &w0): epdf(rv), n(w0.length()), w(w0), Coms(n) {};
    00319         void set_parameters( int &i, double wi, epdf* ep){w(i)=wi;Coms(i)=ep;}
    00320         vec mean(){vec pom; for(int i=0;i<n;i++){pom+=Coms(i)->mean()*w(i);} return pom;};
    00321         vec sample() {it_error ( "Not implemented" );return 0;}
    00322 };
    00323 */
    00324
    00326
    00327 class euni: public epdf {
    00328 protected:
    00330         vec low;
    00332         vec high;
    00334         vec distance;
    00336         double nk;
    00338         double lnk;
    00339 public:
    00341         euni ( const RV rv ) :epdf ( rv ) {}
    00342         double eval ( const vec &val ) const  {return nk;}
    00343         double evalpdflog ( const vec &val ) const  {return lnk;}
    00344         vec sample() const {
    00345                 vec smp ( rv.count() );
    00346 #pragma omp critical
    00347                 UniRNG.sample_vector ( rv.count(),smp );
    00348                 return low+elem_mult ( distance,smp );
    00349         }
    00351         void set_parameters ( const vec &low0, const vec &high0 ) {
    00352                 distance = high0-low0;
    00353                 it_assert_debug ( min ( distance ) >0.0,"bad support" );
    00354                 low = low0;
    00355                 high = high0;
    00356                 nk = prod ( 1.0/distance );
    00357                 lnk = log ( nk );
    00358         }
    00359         vec mean() const {vec pom=high; pom-=low; pom/=2.0; return pom;}
    00360 };
    00361
    00362
    00368 template<class sq_T>
    00369 class mlnorm : public mEF {
    00371         enorm<sq_T> epdf;
    00372         mat A;
    00373         vec mu_const;
    00374         vec& _mu; //cached epdf.mu;
    00375 public:
    00377         mlnorm (const RV &rv, const RV &rvc );
    00379         void set_parameters ( const  mat &A, const vec &mu0, const sq_T &R );
    00381         vec samplecond (const vec &cond, double &lik );
    00383         mat samplecond (const vec &cond, vec &lik, int n );
    00385         void condition (const vec &cond );
    00386 };
    00387
    00397 class mgamma : public mEF {
    00398 protected:
    00400         egamma epdf;
    00402         double k;
    00404         vec* _beta;
    00405
    00406 public:
    00408         mgamma ( const RV &rv,const RV &rvc );
    00410         void set_parameters ( double k );
    00411         void condition ( const vec &val ) {*_beta=k/val;};
    00412 };
    00413
    00425 class mgamma_fix : public mgamma {
    00426 protected:
    00428         double l;
    00430         vec refl;
    00431 public:
    00433         mgamma_fix ( const RV &rv,const RV &rvc ) : mgamma ( rv,rvc ),refl ( rv.count() ) {};
    00435         void set_parameters ( double k0 , vec ref0, double l0 ) {
    00436                 mgamma::set_parameters ( k0 );
    00437                 refl=pow ( ref0,1.0-l0 );l=l0;
    00438         };
    00439
    00440         void condition ( const vec &val ) {vec mean=elem_mult ( refl,pow ( val,l ) ); *_beta=k/mean;};
    00441 };
    00442
    00444 enum RESAMPLING_METHOD { MULTINOMIAL = 0, STRATIFIED = 1, SYSTEMATIC = 3 };
    00450 class eEmp: public epdf {
    00451 protected :
    00453         int n;
    00455         vec w;
    00457         Array<vec> samples;
    00458 public:
    00460         eEmp ( const RV &rv0 ,int n0 ) :epdf ( rv0 ),n ( n0 ),w ( n ),samples ( n ) {};
    00462         void set_parameters ( const vec &w0, const epdf* pdf0 );
    00464         void set_samples ( const epdf* pdf0 );
    00466         vec& _w()  {return w;};
    00468         Array<vec>& _samples() {return samples;};
    00470         ivec resample ( RESAMPLING_METHOD method = SYSTEMATIC );
    00472         vec sample() const {it_error ( "Not implemented" );return 0;}
    00474         double evalpdflog ( const vec &val ) const {it_error ( "Not implemented" );return 0.0;}
    00475         vec mean() const {
    00476                 vec pom=zeros ( rv.count() );
    00477                 for ( int i=0;i<n;i++ ) {pom+=samples ( i ) *w ( i );}
    00478                 return pom;
    00479         }
    00480 };
    00481
    00482
    00484
    00485 template<class sq_T>
    00486 enorm<sq_T>::enorm ( const RV &rv ) :eEF ( rv ), mu ( rv.count() ),R ( rv.count() ),dim ( rv.count() ) {};
    00487
    00488 template<class sq_T>
    00489 void enorm<sq_T>::set_parameters ( const vec &mu0, const sq_T &R0 ) {
    00490 //Fixme test dimensions of mu0 and R0;
    00491         mu = mu0;
    00492         R = R0;
    00493 };
    00494
    00495 template<class sq_T>
    00496 void enorm<sq_T>::dupdate ( mat &v, double nu ) {
    00497         //
    00498 };
    00499
    00500 // template<class sq_T>
    00501 // void enorm<sq_T>::tupdate ( double phi, mat &vbar, double nubar ) {
    00502 //      //
    00503 // };
    00504
    00505 template<class sq_T>
    00506 vec enorm<sq_T>::sample() const {
    00507         vec x ( dim );
    00508         NorRNG.sample_vector ( dim,x );
    00509         vec smp = R.sqrt_mult ( x );
    00510
    00511         smp += mu;
    00512         return smp;
    00513 };
    00514
    00515 template<class sq_T>
    00516 mat enorm<sq_T>::sample ( int N ) const {
    00517         mat X ( dim,N );
    00518         vec x ( dim );
    00519         vec pom;
    00520         int i;
    00521
    00522         for ( i=0;i<N;i++ ) {
    00523                 NorRNG.sample_vector ( dim,x );
    00524                 pom = R.sqrt_mult ( x );
    00525                 pom +=mu;
    00526                 X.set_col ( i, pom );
    00527         }
    00528
    00529         return X;
    00530 };
    00531
    00532 template<class sq_T>
    00533 double enorm<sq_T>::eval ( const vec &val ) const {
    00534         double pdfl,e;
    00535         pdfl = evalpdflog ( val );
    00536         e = exp ( pdfl );
    00537         return e;
    00538 };
    00539
    00540 template<class sq_T>
    00541 double enorm<sq_T>::evalpdflog_nn ( const vec &val ) const {
    00542         // 1.83787706640935 = log(2pi)
    00543         return  -0.5* ( R.invqform ( mu-val ) );// - lognc();
    00544 };
    00545
    00546 template<class sq_T>
    00547 inline double enorm<sq_T>::lognc () const {
    00548         // 1.83787706640935 = log(2pi)
    00549         return 0.5* ( R.cols() * 1.83787706640935 +R.logdet() );
    00550 };
    00551
    00552 template<class sq_T>
    00553 mlnorm<sq_T>::mlnorm ( const RV &rv0, const RV &rvc0 ) :mEF ( rv0,rvc0 ),epdf ( rv0 ),A ( rv0.count(),rv0.count() ),_mu ( epdf._mu() ) {
    00554         ep =&epdf;
    00555 }
    00556
    00557 template<class sq_T>
    00558 void mlnorm<sq_T>::set_parameters ( const mat &A0, const vec &mu0, const sq_T &R0 ) {
    00559         epdf.set_parameters ( zeros ( rv.count() ),R0 );
    00560         A = A0;
    00561         mu_const = mu0;
    00562 }
    00563
    00564 template<class sq_T>
    00565 vec mlnorm<sq_T>::samplecond (const  vec &cond, double &lik ) {
    00566         this->condition ( cond );
    00567         vec smp = epdf.sample();
    00568         lik = epdf.eval ( smp );
    00569         return smp;
    00570 }
    00571
    00572 template<class sq_T>
    00573 mat mlnorm<sq_T>::samplecond (const vec &cond, vec &lik, int n ) {
    00574         int i;
    00575         int dim = rv.count();
    00576         mat Smp ( dim,n );
    00577         vec smp ( dim );
    00578         this->condition ( cond );
    00579
    00580         for ( i=0; i<n; i++ ) {
    00581                 smp = epdf.sample();
    00582                 lik ( i ) = epdf.eval ( smp );
    00583                 Smp.set_col ( i ,smp );
    00584         }
    00585
    00586         return Smp;
    00587 }
    00588
    00589 template<class sq_T>
    00590 void mlnorm<sq_T>::condition (const vec &cond ) {
    00591         _mu = A*cond + mu_const;
    00592 //R is already assigned;
    00593 }
    00594
    00595 template<class sq_T>
    00596 enorm<sq_T>* enorm<sq_T>::marginal ( const RV &rvn ) {
    00597         ivec irvn = rvn.dataind ( rv );
    00598
    00599         sq_T Rn ( R,irvn );
    00600         enorm<sq_T>* tmp = new enorm<sq_T>( rvn );
    00601         tmp->set_parameters ( mu ( irvn ), Rn );
    00602         return tmp;
    00603 }
    00604
    00605 template<class sq_T>
    00606 mlnorm<sq_T>* enorm<sq_T>::condition ( const RV &rvn ) {
    00607
    00608         RV rvc = rv.subt ( rvn );
    00609         it_assert_debug ( ( rvc.count() +rvn.count() ==rv.count() ),"wrong rvn" );
    00610         //Permutation vector of the new R
    00611         ivec irvn = rvn.dataind ( rv );
    00612         ivec irvc = rvc.dataind ( rv );
    00613         ivec perm=concat ( irvn , irvc );
    00614         sq_T Rn ( R,perm );
    00615
    00616         //fixme - could this be done in general for all sq_T?
    00617         mat S=R.to_mat();
    00618         //fixme
    00619         int n=rvn.count()-1;
    00620         int end=R.rows()-1;
    00621         mat S11 = S.get ( 0,n, 0, n );
    00622         mat S12 = S.get ( rvn.count(), end, 0, n );
    00623         mat S22 = S.get ( rvn.count(), end, rvn.count(), end );
    00624
    00625         vec mu1 = mu ( irvn );
    00626         vec mu2 = mu ( irvc );
    00627         mat A=S12*inv ( S22 );
    00628         sq_T R_n ( S11 - A *S12.T() );
    00629
    00630         mlnorm<sq_T>* tmp=new mlnorm<sq_T> ( rvn,rvc );
    00631
    00632         tmp->set_parameters ( A,mu1-A*mu2,R_n );
    00633         return tmp;
    00634 }
    00635
    00637
    00638
    00639 #endif //EF_H