- Main Page
- Related Pages
- Classes
- Files

# work/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 tupdate ( double phi, mat &vbar, double nubar ) {};
    00048         virtual void dupdate ( mat &v,double nu=1.0 ) {};
    00049 };
    00050
    00057 class mEF : public mpdf {
    00058
    00059 public:
    00061         mEF ( const RV &rv0, const RV &rvc0 ) :mpdf ( rv0,rvc0 ) {};
    00062 };
    00063
    00069 template<class sq_T>
    00070
    00071 class enorm : public eEF {
    00072 protected:
    00074         vec mu;
    00076         sq_T R;
    00078         int dim;
    00079 public:
    00080 //      enorm() :eEF() {};
    00082         enorm ( RV &rv );
    00084         void set_parameters ( const vec &mu,const sq_T &R );
    00086         void tupdate ( double phi, mat &vbar, double nubar );
    00088         void dupdate ( mat &v,double nu=1.0 );
    00089
    00090         vec sample() const;
    00092         mat sample ( int N ) const;
    00093         double eval ( const vec &val ) const ;
    00094         double evalpdflog ( const vec &val ) const;
    00095         double lognc () const;
    00096         vec mean() const {return mu;}
    00097
    00098 //Access methods
    00100         vec& _mu() {return mu;}
    00101
    00103         void set_mu(const vec mu0) { mu=mu0;}
    00104
    00106         sq_T& _R() {return R;}
    00107
    00109         mat getR () {return R.to_mat();}
    00110 };
    00111
    00117 class egiw : public eEF {
    00118 protected:
    00120         ldmat V;
    00122         double nu;
    00123 public:
    00125         egiw(RV rv, mat V0, double nu0): eEF(rv), V(V0), nu(nu0) {
    00126                 it_assert_debug(rv.count()==V.rows(),"Incompatible V0.");
    00127         }
    00128
    00129         vec sample() const;
    00130         vec mean() const;
    00131         double evalpdflog ( const vec &val ) const;
    00132         double lognc () const;
    00133
    00134         //Access
    00136         ldmat& _V() {return V;}
    00138         double& _nu() {return nu;}
    00139
    00140 };
    00141
    00151 class egamma : public eEF {
    00152 protected:
    00154         vec alpha;
    00156         vec beta;
    00157 public :
    00159         egamma ( const RV &rv ) :eEF ( rv ) {};
    00161         void set_parameters ( const vec &a, const vec &b ) {alpha=a,beta=b;};
    00162         vec sample() const;
    00164 //      mat sample ( int N ) const;
    00165         double evalpdflog ( const vec &val ) const;
    00166         double lognc () const;
    00168         void _param ( vec* &a, vec* &b ) {a=&alpha;b=&beta;};
    00169         vec mean() const {vec pom ( alpha ); pom/=beta; return pom;}
    00170 };
    00171 /*
    00173 class emix : public epdf {
    00174 protected:
    00175         int n;
    00176         vec &w;
    00177         Array<epdf*> Coms;
    00178 public:
    00180         emix ( const RV &rv, vec &w0): epdf(rv), n(w0.length()), w(w0), Coms(n) {};
    00181         void set_parameters( int &i, double wi, epdf* ep){w(i)=wi;Coms(i)=ep;}
    00182         vec mean(){vec pom; for(int i=0;i<n;i++){pom+=Coms(i)->mean()*w(i);} return pom;};
    00183         vec sample() {it_error ( "Not implemented" );return 0;}
    00184 };
    00185 */
    00186
    00188
    00189 class euni: public epdf {
    00190 protected:
    00192         vec low;
    00194         vec high;
    00196         vec distance;
    00198         double nk;
    00200         double lnk;
    00201 public:
    00203         euni ( const RV rv ) :epdf ( rv ) {}
    00204         double eval ( const vec &val ) const  {return nk;}
    00205         double evalpdflog ( const vec &val ) const  {return lnk;}
    00206         vec sample() const {
    00207                 vec smp ( rv.count() );
    00208                 #pragma omp critical
    00209                 UniRNG.sample_vector ( rv.count(),smp );
    00210                 return low+elem_mult(distance,smp);
    00211         }
    00213         void set_parameters ( const vec &low0, const vec &high0 ) {
    00214                 distance = high0-low0;
    00215                 it_assert_debug ( min ( distance ) >0.0,"bad support" );
    00216                 low = low0;
    00217                 high = high0;
    00218                 nk = prod ( 1.0/distance );
    00219                 lnk = log ( nk );
    00220         }
    00221         vec mean() const {vec pom=high; pom-=low; pom/=2.0; return pom;}
    00222 };
    00223
    00224
    00230 template<class sq_T>
    00231 class mlnorm : public mEF {
    00233         enorm<sq_T> epdf;
    00234         mat A;
    00235         vec& _mu; //cached epdf.mu;
    00236 public:
    00238         mlnorm ( RV &rv,RV &rvc );
    00240         void set_parameters ( const  mat &A, const sq_T &R );
    00242         vec samplecond ( vec &cond, double &lik );
    00244         mat samplecond ( vec &cond, vec &lik, int n );
    00246         void condition ( vec &cond );
    00247 };
    00248
    00258 class mgamma : public mEF {
    00259 protected:
    00261         egamma epdf;
    00263         double k;
    00265         vec* _beta;
    00266
    00267 public:
    00269         mgamma ( const RV &rv,const RV &rvc );
    00271         void set_parameters ( double k );
    00273         vec samplecond ( vec &cond, double &lik );
    00275         mat samplecond ( vec &cond, vec &lik, int n );
    00276         void condition ( const vec &val ) {*_beta=k/val;};
    00277 };
    00278
    00290 class mgamma_fix : public mgamma {
    00291 protected:
    00293         double l;
    00295         vec refl;
    00296 public:
    00298         mgamma_fix ( const RV &rv,const RV &rvc ) : mgamma ( rv,rvc ),refl ( rv.count() ) {};
    00300         void set_parameters ( double k0 , vec ref0, double l0 ) {
    00301                 mgamma::set_parameters ( k0 );
    00302                 refl=pow ( ref0,1.0-l0 );l=l0;
    00303         };
    00304
    00305         void condition ( const vec &val ) {vec mean=elem_mult ( refl,pow ( val,l ) ); *_beta=k/mean;};
    00306 };
    00307
    00309 enum RESAMPLING_METHOD { MULTINOMIAL = 0, STRATIFIED = 1, SYSTEMATIC = 3 };
    00315 class eEmp: public epdf {
    00316 protected :
    00318         int n;
    00320         vec w;
    00322         Array<vec> samples;
    00323 public:
    00325         eEmp ( const RV &rv0 ,int n0 ) :epdf ( rv0 ),n ( n0 ),w ( n ),samples ( n ) {};
    00327         void set_parameters ( const vec &w0, epdf* pdf0 );
    00329         vec& _w()  {return w;};
    00331         Array<vec>& _samples() {return samples;};
    00333         ivec resample ( RESAMPLING_METHOD method = SYSTEMATIC );
    00335         vec sample() const {it_error ( "Not implemented" );return 0;}
    00337         double evalpdflog ( const vec &val ) const {it_error ( "Not implemented" );return 0.0;}
    00338         vec mean() const {
    00339                 vec pom=zeros ( rv.count() );
    00340                 for ( int i=0;i<n;i++ ) {pom+=samples ( i ) *w ( i );}
    00341                 return pom;
    00342         }
    00343 };
    00344
    00345
    00347
    00348 template<class sq_T>
    00349 enorm<sq_T>::enorm ( RV &rv ) :eEF ( rv ), mu ( rv.count() ),R ( rv.count() ),dim ( rv.count() ) {};
    00350
    00351 template<class sq_T>
    00352 void enorm<sq_T>::set_parameters ( const vec &mu0, const sq_T &R0 ) {
    00353 //Fixme test dimensions of mu0 and R0;
    00354         mu = mu0;
    00355         R = R0;
    00356 };
    00357
    00358 template<class sq_T>
    00359 void enorm<sq_T>::dupdate ( mat &v, double nu ) {
    00360         //
    00361 };
    00362
    00363 template<class sq_T>
    00364 void enorm<sq_T>::tupdate ( double phi, mat &vbar, double nubar ) {
    00365         //
    00366 };
    00367
    00368 template<class sq_T>
    00369 vec enorm<sq_T>::sample() const {
    00370         vec x ( dim );
    00371         NorRNG.sample_vector ( dim,x );
    00372         vec smp = R.sqrt_mult ( x );
    00373
    00374         smp += mu;
    00375         return smp;
    00376 };
    00377
    00378 template<class sq_T>
    00379 mat enorm<sq_T>::sample ( int N ) const {
    00380         mat X ( dim,N );
    00381         vec x ( dim );
    00382         vec pom;
    00383         int i;
    00384
    00385         for ( i=0;i<N;i++ ) {
    00386                 NorRNG.sample_vector ( dim,x );
    00387                 pom = R.sqrt_mult ( x );
    00388                 pom +=mu;
    00389                 X.set_col ( i, pom );
    00390         }
    00391
    00392         return X;
    00393 };
    00394
    00395 template<class sq_T>
    00396 double enorm<sq_T>::eval ( const vec &val ) const {
    00397         double pdfl,e;
    00398         pdfl = evalpdflog ( val );
    00399         e = exp ( pdfl );
    00400         return e;
    00401 };
    00402
    00403 template<class sq_T>
    00404 double enorm<sq_T>::evalpdflog ( const vec &val ) const {
    00405         // 1.83787706640935 = log(2pi)
    00406         return  -0.5* (  +R.invqform ( mu-val ) ) - lognc();
    00407 };
    00408
    00409 template<class sq_T>
    00410 inline double enorm<sq_T>::lognc () const {
    00411         // 1.83787706640935 = log(2pi)
    00412         return -0.5* ( R.cols() * 1.83787706640935 +R.logdet());
    00413 };
    00414
    00415 template<class sq_T>
    00416 mlnorm<sq_T>::mlnorm ( RV &rv0,RV &rvc0 ) :mEF ( rv0,rvc0 ),epdf ( rv0 ),A ( rv0.count(),rv0.count() ),_mu(epdf._mu()) { ep =&epdf;
    00417 }
    00418
    00419 template<class sq_T>
    00420 void mlnorm<sq_T>::set_parameters ( const mat &A0, const sq_T &R0 ) {
    00421         epdf.set_parameters ( zeros ( rv.count() ),R0 );
    00422         A = A0;
    00423 }
    00424
    00425 template<class sq_T>
    00426 vec mlnorm<sq_T>::samplecond ( vec &cond, double &lik ) {
    00427         this->condition ( cond );
    00428         vec smp = epdf.sample();
    00429         lik = epdf.eval ( smp );
    00430         return smp;
    00431 }
    00432
    00433 template<class sq_T>
    00434 mat mlnorm<sq_T>::samplecond ( vec &cond, vec &lik, int n ) {
    00435         int i;
    00436         int dim = rv.count();
    00437         mat Smp ( dim,n );
    00438         vec smp ( dim );
    00439         this->condition ( cond );
    00440
    00441         for ( i=0; i<n; i++ ) {
    00442                 smp = epdf.sample();
    00443                 lik ( i ) = epdf.eval ( smp );
    00444                 Smp.set_col ( i ,smp );
    00445         }
    00446
    00447         return Smp;
    00448 }
    00449
    00450 template<class sq_T>
    00451 void mlnorm<sq_T>::condition ( vec &cond ) {
    00452         _mu = A*cond;
    00453 //R is already assigned;
    00454 }
    00455
    00457
    00458
    00459 #endif //EF_H