- 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         sq_T _iR;
    00079         bool cached;
    00081         int dim;
    00082 public:
    00083 //      enorm() :eEF() {};
    00085         enorm ( RV &rv );
    00087         void set_parameters ( const vec &mu,const sq_T &R );
    00089         void tupdate ( double phi, mat &vbar, double nubar );
    00091         void dupdate ( mat &v,double nu=1.0 );
    00092
    00093         vec sample() const;
    00095         mat sample ( int N ) const;
    00096         double eval ( const vec &val ) const ;
    00097         double evalpdflog ( const vec &val ) const;
    00098         vec mean()const {return mu;}
    00099
    00100 //Access methods
    00102         vec* _mu() {return &mu;}
    00103
    00105         void _R ( sq_T* &pR, sq_T* &piR ) {
    00106                 pR=&R;
    00107                 piR=&_iR;
    00108         }
    00109
    00111         void _cached ( bool what ) {cached=what;}
    00112 };
    00113
    00123 class egamma : public eEF {
    00124 protected:
    00126         vec alpha;
    00128         vec beta;
    00129 public :
    00131         egamma ( const RV &rv ) :eEF ( rv ) {};
    00133         void set_parameters ( const vec &a, const vec &b ) {alpha=a,beta=b;};
    00134         vec sample() const;
    00136         mat sample ( int N ) const;
    00137         double evalpdflog ( const vec &val ) const;
    00139         void _param ( vec* &a, vec* &b ) {a=&alpha;b=&beta;};
    00140         vec mean()const {vec pom(alpha); pom/=beta; return pom;}
    00141 };
    00142 /*
    00144 class emix : public epdf {
    00145 protected:
    00146         int n;
    00147         vec &w;
    00148         Array<epdf*> Coms;
    00149 public:
    00151         emix ( const RV &rv, vec &w0): epdf(rv), n(w0.length()), w(w0), Coms(n) {};
    00152         void set_parameters( int &i, double wi, epdf* ep){w(i)=wi;Coms(i)=ep;}
    00153         vec mean(){vec pom; for(int i=0;i<n;i++){pom+=Coms(i)->mean()*w(i);} return pom;};
    00154         vec sample() {it_error ( "Not implemented" );return 0;}
    00155 };
    00156 */
    00157
    00159
    00160 class euni: public epdf {
    00161 protected:
    00163         vec low;
    00165         vec high;
    00167         vec distance;
    00169         double nk;
    00171         double lnk;
    00172 public:
    00174         euni ( const RV rv ) :epdf ( rv ) {}
    00175         double eval ( const vec &val ) const  {return nk;}
    00176         double evalpdflog ( const vec &val ) const  {return lnk;}
    00177         vec sample() const {
    00178                 vec smp ( rv.count() ); UniRNG.sample_vector ( rv.count(),smp );
    00179                 return low+distance*smp;
    00180         }
    00182         void set_parameters ( const vec &low0, const vec &high0 ) {
    00183                 distance = high0-low0;
    00184                 it_assert_debug ( min ( distance ) >0.0,"bad support" );
    00185                 low = low0;
    00186                 high = high0;
    00187                 nk = prod ( 1.0/distance );
    00188                 lnk = log ( nk );
    00189         }
    00190         vec mean() const {vec pom=high; pom-=low; pom/=2.0; return pom;}
    00191 };
    00192
    00193
    00199 template<class sq_T>
    00200 class mlnorm : public mEF {
    00202         enorm<sq_T> epdf;
    00203         vec* _mu; //cached epdf.mu;
    00204         mat A;
    00205 public:
    00207         mlnorm ( RV &rv,RV &rvc );
    00209         void set_parameters ( const  mat &A, const sq_T &R );
    00211         vec samplecond ( vec &cond, double &lik );
    00213         mat samplecond ( vec &cond, vec &lik, int n );
    00215         void condition ( vec &cond );
    00216 };
    00217
    00227 class mgamma : public mEF {
    00229         egamma epdf;
    00231         double k;
    00233         vec* _beta;
    00234
    00235 public:
    00237         mgamma ( const RV &rv,const RV &rvc );
    00239         void set_parameters ( double k );
    00241         vec samplecond ( vec &cond, double &lik );
    00243         mat samplecond ( vec &cond, vec &lik, int n );
    00244         void condition ( const vec &val ) {*_beta=k/val;};
    00245 };
    00246
    00248 enum RESAMPLING_METHOD { MULTINOMIAL = 0, STRATIFIED = 1, SYSTEMATIC = 3 };
    00254 class eEmp: public epdf {
    00255 protected :
    00257         int n;
    00259         vec w;
    00261         Array<vec> samples;
    00262 public:
    00264         eEmp ( const RV &rv0 ,int n0) :epdf ( rv0 ),n(n0),w(n),samples(n) {};
    00266         void set_parameters ( const vec &w0, epdf* pdf0 );
    00268         vec& _w()  {return w;};
    00270         Array<vec>& _samples() {return samples;};
    00272         ivec resample ( RESAMPLING_METHOD method = SYSTEMATIC );
    00274         vec sample() const {it_error ( "Not implemented" );return 0;}
    00276         double evalpdflog(const vec &val) const {it_error ( "Not implemented" );return 0.0;}
    00277         vec mean()const {vec pom=zeros(rv.count());
    00278                 for (int i=0;i<n;i++){pom+=samples(i)*w(i);}
    00279                 return pom;
    00280         }
    00281 };
    00282
    00283
    00285
    00286 template<class sq_T>
    00287 enorm<sq_T>::enorm ( RV &rv ) :eEF(rv), mu ( rv.count() ),R ( rv.count() ),_iR ( rv.count() ),cached ( false ),dim ( rv.count() ) {};
    00288
    00289 template<class sq_T>
    00290 void enorm<sq_T>::set_parameters ( const vec &mu0, const sq_T &R0 ) {
    00291 //Fixme test dimensions of mu0 and R0;
    00292         mu = mu0;
    00293         R = R0;
    00294         if ( _iR.rows() !=R.rows() ) _iR=R; // memory allocation!
    00295         R.inv ( _iR ); //update cache
    00296         cached=true;
    00297 };
    00298
    00299 template<class sq_T>
    00300 void enorm<sq_T>::dupdate ( mat &v, double nu ) {
    00301         //
    00302 };
    00303
    00304 template<class sq_T>
    00305 void enorm<sq_T>::tupdate ( double phi, mat &vbar, double nubar ) {
    00306         //
    00307 };
    00308
    00309 template<class sq_T>
    00310 vec enorm<sq_T>::sample() const {
    00311         vec x ( dim );
    00312         NorRNG.sample_vector ( dim,x );
    00313         vec smp = R.sqrt_mult ( x );
    00314
    00315         smp += mu;
    00316         return smp;
    00317 };
    00318
    00319 template<class sq_T>
    00320 mat enorm<sq_T>::sample ( int N )const {
    00321         mat X ( dim,N );
    00322         vec x ( dim );
    00323         vec pom;
    00324         int i;
    00325
    00326         for ( i=0;i<N;i++ ) {
    00327                 NorRNG.sample_vector ( dim,x );
    00328                 pom = R.sqrt_mult ( x );
    00329                 pom +=mu;
    00330                 X.set_col ( i, pom );
    00331         }
    00332
    00333         return X;
    00334 };
    00335
    00336 template<class sq_T>
    00337 double enorm<sq_T>::eval ( const vec &val ) const {
    00338         double pdfl,e;
    00339         pdfl = evalpdflog ( val );
    00340         e = exp ( pdfl );
    00341         return e;
    00342 };
    00343
    00344 template<class sq_T>
    00345 double enorm<sq_T>::evalpdflog ( const vec &val ) const {
    00346         if ( !cached ) {it_error("this should not happen, see cached");}
    00347
    00348         // 1.83787706640935 = log(2pi)
    00349         return -0.5* ( R.cols() * 1.83787706640935 +R.logdet() +_iR.qform ( mu-val ) );
    00350 };
    00351
    00352
    00353 template<class sq_T>
    00354 mlnorm<sq_T>::mlnorm ( RV &rv0,RV &rvc0 ) :mEF ( rv0,rvc0 ),epdf ( rv ),A ( rv0.count(),rv0.count() ) {
    00355         _mu = epdf._mu();
    00356 }
    00357
    00358 template<class sq_T>
    00359 void mlnorm<sq_T>::set_parameters ( const mat &A0, const sq_T &R0 ) {
    00360         epdf.set_parameters ( zeros ( rv.count() ),R0 );
    00361         A = A0;
    00362 }
    00363
    00364 template<class sq_T>
    00365 vec mlnorm<sq_T>::samplecond ( vec &cond, double &lik ) {
    00366         this->condition ( cond );
    00367         vec smp = epdf.sample();
    00368         lik = epdf.eval ( smp );
    00369         return smp;
    00370 }
    00371
    00372 template<class sq_T>
    00373 mat mlnorm<sq_T>::samplecond ( vec &cond, vec &lik, int n ) {
    00374         int i;
    00375         int dim = rv.count();
    00376         mat Smp ( dim,n );
    00377         vec smp ( dim );
    00378         this->condition ( cond );
    00379
    00380         for ( i=0; i<n; i++ ) {
    00381                 smp = epdf.sample();
    00382                 lik ( i ) = epdf.eval ( smp );
    00383                 Smp.set_col ( i ,smp );
    00384         }
    00385
    00386         return Smp;
    00387 }
    00388
    00389 template<class sq_T>
    00390 void mlnorm<sq_T>::condition ( vec &cond ) {
    00391         *_mu = A*cond;
    00392 //R is already assigned;
    00393 }
    00394
    00396
    00397
    00398 #endif //EF_H