- Main Page
- Classes
- Files
- Related Pages

# 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
    00040 public:
    00041 //      eEF() :epdf() {};
    00043         eEF ( const RV &rv ) :epdf ( rv ) {};
    00045         virtual void tupdate ( double phi, mat &vbar, double nubar ) {};
    00047         virtual void dupdate ( mat &v,double nu=1.0 ) {};
    00048 };
    00049
    00056 class mEF : public mpdf {
    00057
    00058 public:
    00060         mEF ( const RV &rv0, const RV &rvc0 ) :mpdf ( rv0,rvc0 ) {};
    00061 };
    00062
    00068 template<class sq_T>
    00069
    00070 class enorm : public eEF {
    00071 protected:
    00073         vec mu;
    00075         sq_T R;
    00077         int dim;
    00078 public:
    00079 //      enorm() :eEF() {};
    00081         enorm ( RV &rv );
    00083         void set_parameters ( const vec &mu,const sq_T &R );
    00085         void tupdate ( double phi, mat &vbar, double nubar );
    00087         void dupdate ( mat &v,double nu=1.0 );
    00088
    00089         vec sample() const;
    00091         mat sample ( int N ) const;
    00092         double eval ( const vec &val ) const ;
    00093         double evalpdflog ( const vec &val ) const;
    00094         vec mean() const {return mu;}
    00095
    00096 //Access methods
    00098         vec& _mu() {return mu;}
    00099
    00101         void set_mu(const vec mu0) { mu=mu0;}
    00102
    00104         sq_T& _R() {return R;}
    00105
    00107         mat getR () {return R.to_mat();}
    00108 };
    00109
    00119 class egamma : public eEF {
    00120 protected:
    00122         vec alpha;
    00124         vec beta;
    00125 public :
    00127         egamma ( const RV &rv ) :eEF ( rv ) {};
    00129         void set_parameters ( const vec &a, const vec &b ) {alpha=a,beta=b;};
    00130         vec sample() const;
    00132         mat sample ( int N ) const;
    00133         double evalpdflog ( const vec &val ) const;
    00135         void _param ( vec* &a, vec* &b ) {a=&alpha;b=&beta;};
    00136         vec mean() const {vec pom ( alpha ); pom/=beta; return pom;}
    00137 };
    00138 /*
    00140 class emix : public epdf {
    00141 protected:
    00142         int n;
    00143         vec &w;
    00144         Array<epdf*> Coms;
    00145 public:
    00147         emix ( const RV &rv, vec &w0): epdf(rv), n(w0.length()), w(w0), Coms(n) {};
    00148         void set_parameters( int &i, double wi, epdf* ep){w(i)=wi;Coms(i)=ep;}
    00149         vec mean(){vec pom; for(int i=0;i<n;i++){pom+=Coms(i)->mean()*w(i);} return pom;};
    00150         vec sample() {it_error ( "Not implemented" );return 0;}
    00151 };
    00152 */
    00153
    00155
    00156 class euni: public epdf {
    00157 protected:
    00159         vec low;
    00161         vec high;
    00163         vec distance;
    00165         double nk;
    00167         double lnk;
    00168 public:
    00170         euni ( const RV rv ) :epdf ( rv ) {}
    00171         double eval ( const vec &val ) const  {return nk;}
    00172         double evalpdflog ( const vec &val ) const  {return lnk;}
    00173         vec sample() const {
    00174                 vec smp ( rv.count() ); UniRNG.sample_vector ( rv.count(),smp );
    00175                 return low+distance*smp;
    00176         }
    00178         void set_parameters ( const vec &low0, const vec &high0 ) {
    00179                 distance = high0-low0;
    00180                 it_assert_debug ( min ( distance ) >0.0,"bad support" );
    00181                 low = low0;
    00182                 high = high0;
    00183                 nk = prod ( 1.0/distance );
    00184                 lnk = log ( nk );
    00185         }
    00186         vec mean() const {vec pom=high; pom-=low; pom/=2.0; return pom;}
    00187 };
    00188
    00189
    00195 template<class sq_T>
    00196 class mlnorm : public mEF {
    00198         enorm<sq_T> epdf;
    00199         vec& _mu; //cached epdf.mu;
    00200         mat A;
    00201 public:
    00203         mlnorm ( RV &rv,RV &rvc );
    00205         void set_parameters ( const  mat &A, const sq_T &R );
    00207         vec samplecond ( vec &cond, double &lik );
    00209         mat samplecond ( vec &cond, vec &lik, int n );
    00211         void condition ( vec &cond );
    00212 };
    00213
    00223 class mgamma : public mEF {
    00224 protected:
    00226         egamma epdf;
    00228         double k;
    00230         vec* _beta;
    00231
    00232 public:
    00234         mgamma ( const RV &rv,const RV &rvc );
    00236         void set_parameters ( double k );
    00238         vec samplecond ( vec &cond, double &lik );
    00240         mat samplecond ( vec &cond, vec &lik, int n );
    00241         void condition ( const vec &val ) {*_beta=k/val;};
    00242 };
    00243
    00255 class mgamma_fix : public mgamma {
    00256 protected:
    00257         double l;
    00258         vec refl;
    00259 public:
    00261         mgamma_fix ( const RV &rv,const RV &rvc ) : mgamma ( rv,rvc ),refl ( rv.count() ) {};
    00263         void set_parameters ( double k0 , vec ref0, double l0 ) {
    00264                 mgamma::set_parameters ( k0 );
    00265                 refl=pow ( ref0,1.0-l0 );l=l0;
    00266         };
    00267
    00268         void condition ( const vec &val ) {vec mean=elem_mult ( refl,pow ( val,l ) ); *_beta=k/mean;};
    00269 };
    00270
    00272 enum RESAMPLING_METHOD { MULTINOMIAL = 0, STRATIFIED = 1, SYSTEMATIC = 3 };
    00278 class eEmp: public epdf {
    00279 protected :
    00281         int n;
    00283         vec w;
    00285         Array<vec> samples;
    00286 public:
    00288         eEmp ( const RV &rv0 ,int n0 ) :epdf ( rv0 ),n ( n0 ),w ( n ),samples ( n ) {};
    00290         void set_parameters ( const vec &w0, epdf* pdf0 );
    00292         vec& _w()  {return w;};
    00294         Array<vec>& _samples() {return samples;};
    00296         ivec resample ( RESAMPLING_METHOD method = SYSTEMATIC );
    00298         vec sample() const {it_error ( "Not implemented" );return 0;}
    00300         double evalpdflog ( const vec &val ) const {it_error ( "Not implemented" );return 0.0;}
    00301         vec mean() const {
    00302                 vec pom=zeros ( rv.count() );
    00303                 for ( int i=0;i<n;i++ ) {pom+=samples ( i ) *w ( i );}
    00304                 return pom;
    00305         }
    00306 };
    00307
    00308
    00310
    00311 template<class sq_T>
    00312 enorm<sq_T>::enorm ( RV &rv ) :eEF ( rv ), mu ( rv.count() ),R ( rv.count() ),dim ( rv.count() ) {};
    00313
    00314 template<class sq_T>
    00315 void enorm<sq_T>::set_parameters ( const vec &mu0, const sq_T &R0 ) {
    00316 //Fixme test dimensions of mu0 and R0;
    00317         mu = mu0;
    00318         R = R0;
    00319 };
    00320
    00321 template<class sq_T>
    00322 void enorm<sq_T>::dupdate ( mat &v, double nu ) {
    00323         //
    00324 };
    00325
    00326 template<class sq_T>
    00327 void enorm<sq_T>::tupdate ( double phi, mat &vbar, double nubar ) {
    00328         //
    00329 };
    00330
    00331 template<class sq_T>
    00332 vec enorm<sq_T>::sample() const {
    00333         vec x ( dim );
    00334         NorRNG.sample_vector ( dim,x );
    00335         vec smp = R.sqrt_mult ( x );
    00336
    00337         smp += mu;
    00338         return smp;
    00339 };
    00340
    00341 template<class sq_T>
    00342 mat enorm<sq_T>::sample ( int N ) const {
    00343         mat X ( dim,N );
    00344         vec x ( dim );
    00345         vec pom;
    00346         int i;
    00347
    00348         for ( i=0;i<N;i++ ) {
    00349                 NorRNG.sample_vector ( dim,x );
    00350                 pom = R.sqrt_mult ( x );
    00351                 pom +=mu;
    00352                 X.set_col ( i, pom );
    00353         }
    00354
    00355         return X;
    00356 };
    00357
    00358 template<class sq_T>
    00359 double enorm<sq_T>::eval ( const vec &val ) const {
    00360         double pdfl,e;
    00361         pdfl = evalpdflog ( val );
    00362         e = exp ( pdfl );
    00363         return e;
    00364 };
    00365
    00366 template<class sq_T>
    00367 double enorm<sq_T>::evalpdflog ( const vec &val ) const {
    00368         // 1.83787706640935 = log(2pi)
    00369         return -0.5* ( R.cols() * 1.83787706640935 +R.logdet() +R.invqform ( mu-val ) );
    00370 };
    00371
    00372
    00373 template<class sq_T>
    00374 mlnorm<sq_T>::mlnorm ( RV &rv0,RV &rvc0 ) :mEF ( rv0,rvc0 ),epdf ( rv ),A ( rv0.count(),rv0.count() ),_mu(epdf._mu()) {
    00375 }
    00376
    00377 template<class sq_T>
    00378 void mlnorm<sq_T>::set_parameters ( const mat &A0, const sq_T &R0 ) {
    00379         epdf.set_parameters ( zeros ( rv.count() ),R0 );
    00380         A = A0;
    00381 }
    00382
    00383 template<class sq_T>
    00384 vec mlnorm<sq_T>::samplecond ( vec &cond, double &lik ) {
    00385         this->condition ( cond );
    00386         vec smp = epdf.sample();
    00387         lik = epdf.eval ( smp );
    00388         return smp;
    00389 }
    00390
    00391 template<class sq_T>
    00392 mat mlnorm<sq_T>::samplecond ( vec &cond, vec &lik, int n ) {
    00393         int i;
    00394         int dim = rv.count();
    00395         mat Smp ( dim,n );
    00396         vec smp ( dim );
    00397         this->condition ( cond );
    00398
    00399         for ( i=0; i<n; i++ ) {
    00400                 smp = epdf.sample();
    00401                 lik ( i ) = epdf.eval ( smp );
    00402                 Smp.set_col ( i ,smp );
    00403         }
    00404
    00405         return Smp;
    00406 }
    00407
    00408 template<class sq_T>
    00409 void mlnorm<sq_T>::condition ( vec &cond ) {
    00410         _mu = A*cond;
    00411 //R is already assigned;
    00412 }
    00413
    00415
    00416
    00417 #endif //EF_H