- Main Page
- Related Pages
- Modules
- Classes
- Files

- File List
- File Members

# libEF.h

Go to the documentation of this file.

    00001
    00013 #ifndef EF_H
    00014 #define EF_H
    00015
    00016
    00017 #include "libBM.h"
    00018 #include "../math/chmat.h"
    00019 //#include "../user_info.h"
    00020
    00021 namespace bdm
    00022 {
    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 ( ) :epdf ( ) {};
    00045                         virtual double lognc() const =0;
    00047                         virtual void dupdate ( mat &v ) {it_error ( "Not implemented" );};
    00049                         virtual double evallog_nn ( const vec &val ) const{it_error ( "Not implemented" );return 0.0;};
    00051                         virtual double evallog ( const vec &val ) const {double tmp;tmp= evallog_nn ( val )-lognc();it_assert_debug ( std::isfinite ( tmp ),"Infinite value" ); return tmp;}
    00053                         virtual vec evallog ( const mat &Val ) const
    00054                         {
    00055                                 vec x ( Val.cols() );
    00056                                 for ( int i=0;i<Val.cols();i++ ) {x ( i ) =evallog_nn ( Val.get_col ( i ) ) ;}
    00057                                 return x-lognc();
    00058                         }
    00060                         virtual void pow ( double p ) {it_error ( "Not implemented" );};
    00061         };
    00062
    00069         class mEF : public mpdf
    00070         {
    00071
    00072                 public:
    00074                         mEF ( ) :mpdf ( ) {};
    00075         };
    00076
    00078         class BMEF : public BM
    00079         {
    00080                 protected:
    00082                         double frg;
    00084                         double last_lognc;
    00085                 public:
    00087                         BMEF ( double frg0=1.0 ) :BM ( ), frg ( frg0 ) {}
    00089                         BMEF ( const BMEF &B ) :BM ( B ), frg ( B.frg ), last_lognc ( B.last_lognc ) {}
    00091                         virtual void set_statistics ( const BMEF* BM0 ) {it_error ( "Not implemented" );};
    00093                         virtual void bayes ( const vec &data, const double w ) {};
    00094                         //original Bayes
    00095                         void bayes ( const vec &dt );
    00097                         virtual void flatten ( const BMEF * B ) {it_error ( "Not implemented" );}
    00099 //      virtual void flatten ( double nu0 ) {it_error ( "Not implemented" );}
    00100
    00101                         BMEF* _copy_ () const {it_error ( "function _copy_ not implemented for this BM" ); return NULL;};
    00102         };
    00103
    00104         template<class sq_T>
    00105         class mlnorm;
    00106
    00112         template<class sq_T>
    00113         class enorm : public eEF
    00114         {
    00115                 protected:
    00117                         vec mu;
    00119                         sq_T R;
    00120                 public:
    00123
    00124                         enorm ( ) :eEF ( ), mu ( ),R ( ) {};
    00125                         enorm ( const vec &mu,const sq_T &R ) {set_parameters ( mu,R );}
    00126                         void set_parameters ( const vec &mu,const sq_T &R );
    00127                         void from_setting(const Setting &root);
    00129
    00132
    00134                         void dupdate ( mat &v,double nu=1.0 );
    00135
    00136                         vec sample() const;
    00137                         mat sample ( int N ) const;
    00138                         double evallog_nn ( const vec &val ) const;
    00139                         double lognc () const;
    00140                         vec mean() const {return mu;}
    00141                         vec variance() const {return diag ( R.to_mat() );}
    00142 //      mlnorm<sq_T>* condition ( const RV &rvn ) const ; <=========== fails to cmpile. Why?
    00143                         mpdf* condition ( const RV &rvn ) const ;
    00144                         enorm<sq_T>* marginal ( const RV &rv ) const;
    00145 //                      epdf* marginal ( const RV &rv ) const;
    00147
    00150
    00151                         vec& _mu() {return mu;}
    00152                         void set_mu ( const vec mu0 ) { mu=mu0;}
    00153                         sq_T& _R() {return R;}
    00154                         const sq_T& _R() const {return R;}
    00156
    00157         };
    00158
    00165         class egiw : public eEF
    00166         {
    00167                 protected:
    00169                         ldmat V;
    00171                         double nu;
    00173                         int dimx;
    00175                         int nPsi;
    00176                 public:
    00179                         egiw() :eEF() {};
    00180                         egiw ( int dimx0, ldmat V0, double nu0=-1.0 ) :eEF() {set_parameters ( dimx0,V0, nu0 );};
    00181
    00182                         void set_parameters ( int dimx0, ldmat V0, double nu0=-1.0 )
    00183                         {
    00184                                 dimx=dimx0;
    00185                                 nPsi = V0.rows()-dimx;
    00186                                 dim = dimx* ( dimx+nPsi ); // size(R) + size(Theta)
    00187
    00188                                 V=V0;
    00189                                 if ( nu0<0 )
    00190                                 {
    00191                                         nu = 0.1 +nPsi +2*dimx +2; // +2 assures finite expected value of R
    00192                                         // terms before that are sufficient for finite normalization
    00193                                 }
    00194                                 else
    00195                                 {
    00196                                         nu=nu0;
    00197                                 }
    00198                         }
    00200
    00201                         vec sample() const;
    00202                         vec mean() const;
    00203                         vec variance() const;
    00204
    00206                         vec est_theta() const;
    00207
    00209                         ldmat est_theta_cov() const;
    00210
    00211                         void mean_mat ( mat &M, mat&R ) const;
    00213                         double evallog_nn ( const vec &val ) const;
    00214                         double lognc () const;
    00215                         void pow ( double p ) {V*=p;nu*=p;};
    00216
    00219
    00220                         ldmat& _V() {return V;}
    00221                         const ldmat& _V() const {return V;}
    00222                         double& _nu()  {return nu;}
    00223                         const double& _nu() const {return nu;}
    00225         };
    00226
    00235         class eDirich: public eEF
    00236         {
    00237                 protected:
    00239                         vec beta;
    00240                 public:
    00243
    00244                         eDirich () : eEF ( ) {};
    00245                         eDirich ( const eDirich &D0 ) : eEF () {set_parameters ( D0.beta );};
    00246                         eDirich ( const vec &beta0 ) {set_parameters ( beta0 );};
    00247                         void set_parameters ( const vec &beta0 )
    00248                         {
    00249                                 beta= beta0;
    00250                                 dim = beta.length();
    00251                         }
    00253
    00254                         vec sample() const {it_error ( "Not implemented" );return vec_1 ( 0.0 );};
    00255                         vec mean() const {return beta/sum(beta);};
    00256                         vec variance() const {double gamma =sum(beta); return elem_mult ( beta, ( beta+1 ) ) / ( gamma* ( gamma+1 ) );}
    00258                         double evallog_nn ( const vec &val ) const
    00259                         {
    00260                                 double tmp; tmp= ( beta-1 ) *log ( val );               it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
    00261                                 return tmp;
    00262                         };
    00263                         double lognc () const
    00264                         {
    00265                                 double tmp;
    00266                                 double gam=sum ( beta );
    00267                                 double lgb=0.0;
    00268                                 for ( int i=0;i<beta.length();i++ ) {lgb+=lgamma ( beta ( i ) );}
    00269                                 tmp= lgb-lgamma ( gam );
    00270                                 it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
    00271                                 return tmp;
    00272                         };
    00274                         vec& _beta()  {return beta;}
    00276         };
    00277
    00279         class multiBM : public BMEF
    00280         {
    00281                 protected:
    00283                         eDirich est;
    00285                         vec &beta;
    00286                 public:
    00288                         multiBM ( ) : BMEF ( ),est ( ),beta ( est._beta() )
    00289                         {
    00290                                 if ( beta.length() >0 ) {last_lognc=est.lognc();}
    00291                                 else{last_lognc=0.0;}
    00292                         }
    00294                         multiBM ( const multiBM &B ) : BMEF ( B ),est ( B.est ),beta ( est._beta() ) {}
    00296                         void set_statistics ( const BM* mB0 ) {const multiBM* mB=dynamic_cast<const multiBM*> ( mB0 ); beta=mB->beta;}
    00297                         void bayes ( const vec &dt )
    00298                         {
    00299                                 if ( frg<1.0 ) {beta*=frg;last_lognc=est.lognc();}
    00300                                 beta+=dt;
    00301                                 if ( evalll ) {ll=est.lognc()-last_lognc;}
    00302                         }
    00303                         double logpred ( const vec &dt ) const
    00304                         {
    00305                                 eDirich pred ( est );
    00306                                 vec &beta = pred._beta();
    00307
    00308                                 double lll;
    00309                                 if ( frg<1.0 )
    00310                                         {beta*=frg;lll=pred.lognc();}
    00311                                 else
    00312                                         if ( evalll ) {lll=last_lognc;}
    00313                                         else{lll=pred.lognc();}
    00314
    00315                                 beta+=dt;
    00316                                 return pred.lognc()-lll;
    00317                         }
    00318                         void flatten ( const BMEF* B )
    00319                         {
    00320                                 const multiBM* E=dynamic_cast<const multiBM*> ( B );
    00321                                 // sum(beta) should be equal to sum(B.beta)
    00322                                 const vec &Eb=E->beta;//const_cast<multiBM*> ( E )->_beta();
    00323                                 beta*= ( sum ( Eb ) /sum ( beta ) );
    00324                                 if ( evalll ) {last_lognc=est.lognc();}
    00325                         }
    00326                         const epdf& posterior() const {return est;};
    00327                         const eDirich* _e() const {return &est;};
    00328                         void set_parameters ( const vec &beta0 )
    00329                         {
    00330                                 est.set_parameters ( beta0 );
    00331                                 if ( evalll ) {last_lognc=est.lognc();}
    00332                         }
    00333         };
    00334
    00344         class egamma : public eEF
    00345         {
    00346                 protected:
    00348                         vec alpha;
    00350                         vec beta;
    00351                 public :
    00354                         egamma ( ) :eEF ( ), alpha ( 0 ), beta ( 0 ) {};
    00355                         egamma ( const vec &a, const vec &b ) {set_parameters ( a, b );};
    00356                         void set_parameters ( const vec &a, const vec &b ) {alpha=a,beta=b;dim = alpha.length();};
    00358
    00359                         vec sample() const;
    00361 //      mat sample ( int N ) const;
    00362                         double evallog ( const vec &val ) const;
    00363                         double lognc () const;
    00365                         vec& _alpha() {return alpha;}
    00366                         vec& _beta() {return beta;}
    00367                         vec mean() const {return elem_div ( alpha,beta );}
    00368                         vec variance() const {return elem_div ( alpha,elem_mult ( beta,beta ) ); }
    00369         };
    00370
    00387         class eigamma : public egamma
    00388         {
    00389                 protected:
    00390                 public :
    00395
    00396                         vec sample() const {return 1.0/egamma::sample();};
    00398                         vec mean() const {return elem_div ( beta,alpha-1 );}
    00399                         vec variance() const {vec mea=mean(); return elem_div ( elem_mult ( mea,mea ),alpha-2 );}
    00400         };
    00401         /*
    00403         class emix : public epdf {
    00404         protected:
    00405                 int n;
    00406                 vec &w;
    00407                 Array<epdf*> Coms;
    00408         public:
    00410                 emix ( const RV &rv, vec &w0): epdf(rv), n(w0.length()), w(w0), Coms(n) {};
    00411                 void set_parameters( int &i, double wi, epdf* ep){w(i)=wi;Coms(i)=ep;}
    00412                 vec mean(){vec pom; for(int i=0;i<n;i++){pom+=Coms(i)->mean()*w(i);} return pom;};
    00413                 vec sample() {it_error ( "Not implemented" );return 0;}
    00414         };
    00415         */
    00416
    00418
    00419         class euni: public epdf
    00420         {
    00421                 protected:
    00423                         vec low;
    00425                         vec high;
    00427                         vec distance;
    00429                         double nk;
    00431                         double lnk;
    00432                 public:
    00435                         euni ( ) :epdf ( ) {}
    00436                         euni ( const vec &low0, const vec &high0 ) {set_parameters ( low0,high0 );}
    00437                         void set_parameters ( const vec &low0, const vec &high0 )
    00438                         {
    00439                                 distance = high0-low0;
    00440                                 it_assert_debug ( min ( distance ) >0.0,"bad support" );
    00441                                 low = low0;
    00442                                 high = high0;
    00443                                 nk = prod ( 1.0/distance );
    00444                                 lnk = log ( nk );
    00445                                 dim = low.length();
    00446                         }
    00448
    00449                         double eval ( const vec &val ) const  {return nk;}
    00450                         double evallog ( const vec &val ) const  {return lnk;}
    00451                         vec sample() const
    00452                         {
    00453                                 vec smp ( dim );
    00454 #pragma omp critical
    00455                                 UniRNG.sample_vector ( dim ,smp );
    00456                                 return low+elem_mult ( distance,smp );
    00457                         }
    00459                         vec mean() const {return ( high-low ) /2.0;}
    00460                         vec variance() const {return ( pow ( high,2 ) +pow ( low,2 ) +elem_mult ( high,low ) ) /3.0;}
    00461         };
    00462
    00463
    00469         template<class sq_T>
    00470         class mlnorm : public mEF
    00471         {
    00472                 protected:
    00474                         enorm<sq_T> epdf;
    00475                         mat A;
    00476                         vec mu_const;
    00477                         vec& _mu; //cached epdf.mu;
    00478                 public:
    00481                         mlnorm ( ) :mEF (),epdf ( ),A ( ),_mu ( epdf._mu() ) {ep =&epdf; };
    00482                         mlnorm ( const  mat &A, const vec &mu0, const sq_T &R ) :epdf ( ),_mu ( epdf._mu() )
    00483                         {
    00484                                 ep =&epdf; set_parameters ( A,mu0,R );
    00485                         };
    00487                         void set_parameters ( const  mat &A, const vec &mu0, const sq_T &R );
    00490                         void condition ( const vec &cond );
    00491
    00493                         vec& _mu_const() {return mu_const;}
    00495                         mat& _A() {return A;}
    00497                         mat _R() {return epdf._R().to_mat();}
    00498
    00499                         template<class sq_M>
    00500                         friend std::ostream &operator<< ( std::ostream &os,  mlnorm<sq_M> &ml );
    00501         };
    00502
    00504         template<class sq_T>
    00505         class mgnorm : public mEF
    00506         {
    00507                 protected:
    00509                         enorm<sq_T> epdf;
    00510                         vec &mu;
    00511                         fnc* g;
    00512                 public:
    00514                         mgnorm() :mu ( epdf._mu() ) {ep=&epdf;}
    00516                         void set_parameters ( fnc* g0, const sq_T &R0 ) {g=g0; epdf.set_parameters ( zeros ( g->dimension() ), R0 );}
    00517                         void condition ( const vec &cond ) {mu=g->eval ( cond );};
    00518
    00519
    00547                         void from_setting( const Setting &root )
    00548                         {
    00549                                 fnc* g = UI::build<fnc>( root, "g" );
    00550
    00551                                 mat R;
    00552                                 if ( root.exists( "dR" ) )
    00553                                 {
    00554                                         vec dR;
    00555                                         UI::get( dR, root, "dR" );
    00556                                         R=diag(dR);
    00557                                 }
    00558                                 else
    00559                                         UI::get( R, root, "R");
    00560
    00561                                 set_parameters(g,R);
    00562                         }
    00563
    00564                         /*void mgnorm::to_setting( Setting &root ) const
    00565                         {
    00566                                 Transport::to_setting( root );
    00567
    00568                                 Setting &kilometers_setting = root.add("kilometers", Setting::TypeInt );
    00569                                 kilometers_setting = kilometers;
    00570
    00571                                 UI::save( passengers, root, "passengers" );
    00572                         }*/
    00573
    00574         };
    00575
    00576         UIREGISTER(mgnorm<chmat>);
    00577
    00578
    00586         class mlstudent : public mlnorm<ldmat>
    00587         {
    00588                 protected:
    00589                         ldmat Lambda;
    00590                         ldmat &_R;
    00591                         ldmat Re;
    00592                 public:
    00593                         mlstudent ( ) :mlnorm<ldmat> (),
    00594                                         Lambda (),      _R ( epdf._R() ) {}
    00595                         void set_parameters ( const mat &A0, const vec &mu0, const ldmat &R0, const ldmat& Lambda0 )
    00596                         {
    00597                                 it_assert_debug ( A0.rows() ==mu0.length(),"" );
    00598                                 it_assert_debug ( R0.rows() ==A0.rows(),"" );
    00599
    00600                                 epdf.set_parameters ( mu0,Lambda ); //
    00601                                 A = A0;
    00602                                 mu_const = mu0;
    00603                                 Re=R0;
    00604                                 Lambda = Lambda0;
    00605                         }
    00606                         void condition ( const vec &cond )
    00607                         {
    00608                                 _mu = A*cond + mu_const;
    00609                                 double zeta;
    00610                                 //ugly hack!
    00611                                 if ( ( cond.length() +1 ) ==Lambda.rows() )
    00612                                 {
    00613                                         zeta = Lambda.invqform ( concat ( cond, vec_1 ( 1.0 ) ) );
    00614                                 }
    00615                                 else
    00616                                 {
    00617                                         zeta = Lambda.invqform ( cond );
    00618                                 }
    00619                                 _R = Re;
    00620                                 _R*= ( 1+zeta );// / ( nu ); << nu is in Re!!!!!!
    00621                         };
    00622
    00623         };
    00633         class mgamma : public mEF
    00634         {
    00635                 protected:
    00637                         egamma epdf;
    00639                         double k;
    00641                         vec &_beta;
    00642
    00643                 public:
    00645                         mgamma ( ) : mEF ( ), epdf (), _beta ( epdf._beta() ) {ep=&epdf;};
    00647                         void set_parameters ( double k, const vec &beta0 );
    00648                         void condition ( const vec &val ) {_beta=k/val;};
    00649         };
    00650
    00660         class migamma : public mEF
    00661         {
    00662                 protected:
    00664                         eigamma epdf;
    00666                         double k;
    00668                         vec &_alpha;
    00670                         vec &_beta;
    00671
    00672                 public:
    00675                         migamma ( ) : mEF (), epdf ( ), _alpha ( epdf._alpha() ), _beta ( epdf._beta() ) {ep=&epdf;};
    00676                         migamma ( const migamma &m ) : mEF (), epdf ( m.epdf ), _alpha ( epdf._alpha() ), _beta ( epdf._beta() ) {ep=&epdf;};
    00678
    00680                         void set_parameters ( int len, double k0 )
    00681                         {
    00682                                 k=k0;
    00683                                 epdf.set_parameters ( ( 1.0/ ( k*k ) +2.0 ) *ones ( len ) /*alpha*/, ones ( len ) /*beta*/ );
    00684                                 dimc = dimension();
    00685                         };
    00686                         void condition ( const vec &val )
    00687                         {
    00688                                 _beta=elem_mult ( val, ( _alpha-1.0 ) );
    00689                         };
    00690         };
    00691
    00692
    00704         class mgamma_fix : public mgamma
    00705         {
    00706                 protected:
    00708                         double l;
    00710                         vec refl;
    00711                 public:
    00713                         mgamma_fix ( ) : mgamma ( ),refl () {};
    00715                         void set_parameters ( double k0 , vec ref0, double l0 )
    00716                         {
    00717                                 mgamma::set_parameters ( k0, ref0 );
    00718                                 refl=pow ( ref0,1.0-l0 );l=l0;
    00719                                 dimc=dimension();
    00720                         };
    00721
    00722                         void condition ( const vec &val ) {vec mean=elem_mult ( refl,pow ( val,l ) ); _beta=k/mean;};
    00723         };
    00724
    00725
    00738         class migamma_ref : public migamma
    00739         {
    00740                 protected:
    00742                         double l;
    00744                         vec refl;
    00745                 public:
    00747                         migamma_ref ( ) : migamma (),refl ( ) {};
    00749                         void set_parameters ( double k0 , vec ref0, double l0 )
    00750                         {
    00751                                 migamma::set_parameters ( ref0.length(), k0 );
    00752                                 refl=pow ( ref0,1.0-l0 );
    00753                                 l=l0;
    00754                                 dimc = dimension();
    00755                         };
    00756
    00757                         void condition ( const vec &val )
    00758                         {
    00759                                 vec mean=elem_mult ( refl,pow ( val,l ) );
    00760                                 migamma::condition ( mean );
    00761                         };
    00762
    00783                         void from_setting( const Setting &root );
    00784
    00785                         // TODO dodelat void to_setting( Setting &root ) const;
    00786         };
    00787
    00788
    00789         UIREGISTER(migamma_ref);
    00790
    00800         class elognorm: public enorm<ldmat>
    00801         {
    00802                 public:
    00803                         vec sample() const {return exp ( enorm<ldmat>::sample() );};
    00804                         vec mean() const {vec var=enorm<ldmat>::variance();return exp ( mu - 0.5*var );};
    00805
    00806         };
    00807
    00819         class mlognorm : public mpdf
    00820         {
    00821                 protected:
    00822                         elognorm eno;
    00824                         double sig2;
    00826                         vec &mu;
    00827                 public:
    00829                         mlognorm ( ) : eno (), mu ( eno._mu() ) {ep=&eno;};
    00831                         void set_parameters ( int size, double k )
    00832                         {
    00833                                 sig2 = 0.5*log ( k*k+1 );
    00834                                 eno.set_parameters ( zeros ( size ),2*sig2*eye ( size ) );
    00835
    00836                                 dimc = size;
    00837                         };
    00838
    00839                         void condition ( const vec &val )
    00840                         {
    00841                                 mu=log ( val )-sig2;//elem_mult ( refl,pow ( val,l ) );
    00842                         };
    00843
    00862                         void from_setting( const Setting &root );
    00863
    00864                         // TODO dodelat void to_setting( Setting &root ) const;
    00865
    00866         };
    00867
    00868         UIREGISTER(mlognorm);
    00869
    00873         class eWishartCh : public epdf
    00874         {
    00875                 protected:
    00877                         chmat Y;
    00879                         int p;
    00881                         double delta;
    00882                 public:
    00883                         void set_parameters ( const mat &Y0, const double delta0 ) {Y=chmat ( Y0 );delta=delta0; p=Y.rows(); dim = p*p; }
    00884                         mat sample_mat() const
    00885                         {
    00886                                 mat X=zeros ( p,p );
    00887
    00888                                 //sample diagonal
    00889                                 for ( int i=0;i<p;i++ )
    00890                                 {
    00891                                         GamRNG.setup ( 0.5* ( delta-i ) , 0.5 ); // no +1 !! index if from 0
    00892 #pragma omp critical
    00893                                         X ( i,i ) =sqrt ( GamRNG() );
    00894                                 }
    00895                                 //do the rest
    00896                                 for ( int i=0;i<p;i++ )
    00897                                 {
    00898                                         for ( int j=i+1;j<p;j++ )
    00899                                         {
    00900 #pragma omp critical
    00901                                                 X ( i,j ) =NorRNG.sample();
    00902                                         }
    00903                                 }
    00904                                 return X*Y._Ch();// return upper triangular part of the decomposition
    00905                         }
    00906                         vec sample () const
    00907                         {
    00908                                 return vec ( sample_mat()._data(),p*p );
    00909                         }
    00911                         void setY ( const mat &Ch0 ) {copy_vector ( dim,Ch0._data(), Y._Ch()._data() );}
    00913                         void _setY ( const vec &ch0 ) {copy_vector ( dim, ch0._data(), Y._Ch()._data() ); }
    00915                         const chmat& getY()const {return Y;}
    00916         };
    00917
    00918         class eiWishartCh: public epdf
    00919         {
    00920                 protected:
    00921                         eWishartCh W;
    00922                         int p;
    00923                         double delta;
    00924                 public:
    00925                         void set_parameters ( const mat &Y0, const double delta0) {
    00926                                 delta = delta0;
    00927                                 W.set_parameters ( inv ( Y0 ),delta0 );
    00928                                 dim = W.dimension(); p=Y0.rows();
    00929                         }
    00930                         vec sample() const {mat iCh; iCh=inv ( W.sample_mat() ); return vec ( iCh._data(),dim );}
    00931                         void _setY ( const vec &y0 )
    00932                         {
    00933                                 mat Ch ( p,p );
    00934                                 mat iCh ( p,p );
    00935                                 copy_vector ( dim, y0._data(), Ch._data() );
    00936
    00937                                 iCh=inv ( Ch );
    00938                                 W.setY ( iCh );
    00939                         }
    00940                         virtual double evallog ( const vec &val ) const {
    00941                                 chmat X(p);
    00942                                 const chmat& Y=W.getY();
    00943
    00944                                 copy_vector(p*p,val._data(),X._Ch()._data());
    00945                                 chmat iX(p);X.inv(iX);
    00946                                 // compute
    00947 //                              \frac{ |\Psi|^{m/2}|X|^{-(m+p+1)/2}e^{-tr(\Psi X^{-1})/2} }{ 2^{mp/2}\Gamma_p(m/2)},
    00948                                 mat M=Y.to_mat()*iX.to_mat();
    00949
    00950                                 double log1 = 0.5*p*(2*Y.logdet())-0.5*(delta+p+1)*(2*X.logdet())-0.5*trace(M);
    00951                                 //Fixme! Multivariate gamma omitted!! it is ok for sampling, but not otherwise!!
    00952
    00953 /*                              if (0) {
    00954                                         mat XX=X.to_mat();
    00955                                         mat YY=Y.to_mat();
    00956
    00957                                         double log2 = 0.5*p*log(det(YY))-0.5*(delta+p+1)*log(det(XX))-0.5*trace(YY*inv(XX));
    00958                                         cout << log1 << "," << log2 << endl;
    00959                                 }*/
    00960                                 return log1;
    00961                         };
    00962
    00963         };
    00964
    00965         class rwiWishartCh : public mpdf
    00966         {
    00967                 protected:
    00968                         eiWishartCh eiW;
    00970                         double sqd;
    00971                         //reference point for diagonal
    00972                         vec refl;
    00973                         double l;
    00974                         int p;
    00975                 public:
    00976                         void set_parameters ( int p0, double k, vec ref0, double l0  )
    00977                         {
    00978                                 p=p0;
    00979                                 double delta = 2/(k*k)+p+3;
    00980                                 sqd=sqrt ( delta-p-1 );
    00981                                 l=l0;
    00982                                 refl=pow(ref0,1-l);
    00983
    00984                                 eiW.set_parameters ( eye ( p ),delta );
    00985                                 ep=&eiW;
    00986                                 dimc=eiW.dimension();
    00987                         }
    00988                         void condition ( const vec &c ) {
    00989                                 vec z=c;
    00990                                 int ri=0;
    00991                                 for(int i=0;i<p*p;i+=(p+1)){//trace diagonal element
    00992                                         z(i) = pow(z(i),l)*refl(ri);
    00993                                         ri++;
    00994                                 }
    00995
    00996                                 eiW._setY ( sqd*z );
    00997                         }
    00998         };
    00999
    01001         enum RESAMPLING_METHOD { MULTINOMIAL = 0, STRATIFIED = 1, SYSTEMATIC = 3 };
    01007         class eEmp: public epdf
    01008         {
    01009                 protected :
    01011                         int n;
    01013                         vec w;
    01015                         Array<vec> samples;
    01016                 public:
    01019                         eEmp ( ) :epdf ( ),w ( ),samples ( ) {};
    01021                         eEmp ( const eEmp &e ) : epdf ( e ), w ( e.w ), samples ( e.samples ) {};
    01023
    01025                         void set_statistics ( const vec &w0, const epdf* pdf0 );
    01027                         void set_statistics ( const epdf* pdf0 , int n ) {set_statistics ( ones ( n ) /n,pdf0 );};
    01029                         void set_samples ( const epdf* pdf0 );
    01031                         void set_parameters ( int n0, bool copy=true ) {n=n0; w.set_size ( n0,copy );samples.set_size ( n0,copy );};
    01033                         vec& _w()  {return w;};
    01035                         const vec& _w() const {return w;};
    01037                         Array<vec>& _samples() {return samples;};
    01039                         const Array<vec>& _samples() const {return samples;};
    01041                         ivec resample ( RESAMPLING_METHOD method=SYSTEMATIC );
    01043                         vec sample() const {it_error ( "Not implemented" );return 0;}
    01045                         double evallog ( const vec &val ) const {it_error ( "Not implemented" );return 0.0;}
    01046                         vec mean() const
    01047                         {
    01048                                 vec pom=zeros ( dim );
    01049                                 for ( int i=0;i<n;i++ ) {pom+=samples ( i ) *w ( i );}
    01050                                 return pom;
    01051                         }
    01052                         vec variance() const
    01053                         {
    01054                                 vec pom=zeros ( dim );
    01055                                 for ( int i=0;i<n;i++ ) {pom+=pow ( samples ( i ),2 ) *w ( i );}
    01056                                 return pom-pow ( mean(),2 );
    01057                         }
    01059                         void qbounds ( vec &lb, vec &ub, double perc=0.95 ) const
    01060                         {
    01061                                 // lb in inf so than it will be pushed below;
    01062                                 lb.set_size ( dim );
    01063                                 ub.set_size ( dim );
    01064                                 lb = std::numeric_limits<double>::infinity();
    01065                                 ub = -std::numeric_limits<double>::infinity();
    01066                                 int j;
    01067                                 for ( int i=0;i<n;i++ )
    01068                                 {
    01069                                         for ( j=0;j<dim; j++ )
    01070                                         {
    01071                                                 if ( samples ( i ) ( j ) <lb ( j ) ) {lb ( j ) =samples ( i ) ( j );}
    01072                                                 if ( samples ( i ) ( j ) >ub ( j ) ) {ub ( j ) =samples ( i ) ( j );}
    01073                                         }
    01074                                 }
    01075                         }
    01076         };
    01077
    01078
    01080
    01081         template<class sq_T>
    01082         void enorm<sq_T>::set_parameters ( const vec &mu0, const sq_T &R0 )
    01083         {
    01084 //Fixme test dimensions of mu0 and R0;
    01085                 mu = mu0;
    01086                 R = R0;
    01087                 dim = mu0.length();
    01088         };
    01089
    01090         template<class sq_T>
    01091         void enorm<sq_T>::from_setting(const Setting &root){
    01092                 vec mu;
    01093                 UI::get(mu,root,"mu");
    01094                 mat R;
    01095                 UI::get(R,root,"R");
    01096                 set_parameters(mu,R);
    01097
    01098                 RV* r = UI::build<RV>(root,"rv");
    01099                 set_rv(*r);
    01100                 delete r;
    01101         }
    01102
    01103         template<class sq_T>
    01104         void enorm<sq_T>::dupdate ( mat &v, double nu )
    01105         {
    01106                 //
    01107         };
    01108
    01109 // template<class sq_T>
    01110 // void enorm<sq_T>::tupdate ( double phi, mat &vbar, double nubar ) {
    01111 //      //
    01112 // };
    01113
    01114         template<class sq_T>
    01115         vec enorm<sq_T>::sample() const
    01116         {
    01117                 vec x ( dim );
    01118 #pragma omp critical
    01119                 NorRNG.sample_vector ( dim,x );
    01120                 vec smp = R.sqrt_mult ( x );
    01121
    01122                 smp += mu;
    01123                 return smp;
    01124         };
    01125
    01126         template<class sq_T>
    01127         mat enorm<sq_T>::sample ( int N ) const
    01128         {
    01129                 mat X ( dim,N );
    01130                 vec x ( dim );
    01131                 vec pom;
    01132                 int i;
    01133
    01134                 for ( i=0;i<N;i++ )
    01135                 {
    01136 #pragma omp critical
    01137                         NorRNG.sample_vector ( dim,x );
    01138                         pom = R.sqrt_mult ( x );
    01139                         pom +=mu;
    01140                         X.set_col ( i, pom );
    01141                 }
    01142
    01143                 return X;
    01144         };
    01145
    01146 // template<class sq_T>
    01147 // double enorm<sq_T>::eval ( const vec &val ) const {
    01148 //      double pdfl,e;
    01149 //      pdfl = evallog ( val );
    01150 //      e = exp ( pdfl );
    01151 //      return e;
    01152 // };
    01153
    01154         template<class sq_T>
    01155         double enorm<sq_T>::evallog_nn ( const vec &val ) const
    01156         {
    01157                 // 1.83787706640935 = log(2pi)
    01158                 double tmp=-0.5* ( R.invqform ( mu-val ) );// - lognc();
    01159                 return  tmp;
    01160         };
    01161
    01162         template<class sq_T>
    01163         inline double enorm<sq_T>::lognc () const
    01164         {
    01165                 // 1.83787706640935 = log(2pi)
    01166                 double tmp=0.5* ( R.cols() * 1.83787706640935 +R.logdet() );
    01167                 return tmp;
    01168         };
    01169
    01170         template<class sq_T>
    01171         void mlnorm<sq_T>::set_parameters ( const mat &A0, const vec &mu0, const sq_T &R0 )
    01172         {
    01173                 it_assert_debug ( A0.rows() ==mu0.length(),"" );
    01174                 it_assert_debug ( A0.rows() ==R0.rows(),"" );
    01175
    01176                 epdf.set_parameters ( zeros ( A0.rows() ),R0 );
    01177                 A = A0;
    01178                 mu_const = mu0;
    01179                 dimc=A0.cols();
    01180         }
    01181
    01182 // template<class sq_T>
    01183 // vec mlnorm<sq_T>::samplecond (const  vec &cond, double &lik ) {
    01184 //      this->condition ( cond );
    01185 //      vec smp = epdf.sample();
    01186 //      lik = epdf.eval ( smp );
    01187 //      return smp;
    01188 // }
    01189
    01190 // template<class sq_T>
    01191 // mat mlnorm<sq_T>::samplecond (const vec &cond, vec &lik, int n ) {
    01192 //      int i;
    01193 //      int dim = rv.count();
    01194 //      mat Smp ( dim,n );
    01195 //      vec smp ( dim );
    01196 //      this->condition ( cond );
    01197 //
    01198 //      for ( i=0; i<n; i++ ) {
    01199 //              smp = epdf.sample();
    01200 //              lik ( i ) = epdf.eval ( smp );
    01201 //              Smp.set_col ( i ,smp );
    01202 //      }
    01203 //
    01204 //      return Smp;
    01205 // }
    01206
    01207         template<class sq_T>
    01208         void mlnorm<sq_T>::condition ( const vec &cond )
    01209         {
    01210                 _mu = A*cond + mu_const;
    01211 //R is already assigned;
    01212         }
    01213
    01214         template<class sq_T>
    01215         enorm<sq_T>* enorm<sq_T>::marginal ( const RV &rvn ) const
    01216         {
    01217                 it_assert_debug ( isnamed(), "rv description is not assigned" );
    01218                 ivec irvn = rvn.dataind ( rv );
    01219
    01220                 sq_T Rn ( R,irvn ); //select rows and columns of R
    01221
    01222                 enorm<sq_T>* tmp = new enorm<sq_T>;
    01223                 tmp->set_rv ( rvn );
    01224                 tmp->set_parameters ( mu ( irvn ), Rn );
    01225                 return tmp;
    01226         }
    01227
    01228         template<class sq_T>
    01229         mpdf* enorm<sq_T>::condition ( const RV &rvn ) const
    01230         {
    01231
    01232                 it_assert_debug ( isnamed(),"rvs are not assigned" );
    01233
    01234                 RV rvc = rv.subt ( rvn );
    01235                 it_assert_debug ( ( rvc._dsize() +rvn._dsize() ==rv._dsize() ),"wrong rvn" );
    01236                 //Permutation vector of the new R
    01237                 ivec irvn = rvn.dataind ( rv );
    01238                 ivec irvc = rvc.dataind ( rv );
    01239                 ivec perm=concat ( irvn , irvc );
    01240                 sq_T Rn ( R,perm );
    01241
    01242                 //fixme - could this be done in general for all sq_T?
    01243                 mat S=Rn.to_mat();
    01244                 //fixme
    01245                 int n=rvn._dsize()-1;
    01246                 int end=R.rows()-1;
    01247                 mat S11 = S.get ( 0,n, 0, n );
    01248                 mat S12 = S.get ( 0, n , rvn._dsize(), end );
    01249                 mat S22 = S.get ( rvn._dsize(), end, rvn._dsize(), end );
    01250
    01251                 vec mu1 = mu ( irvn );
    01252                 vec mu2 = mu ( irvc );
    01253                 mat A=S12*inv ( S22 );
    01254                 sq_T R_n ( S11 - A *S12.T() );
    01255
    01256                 mlnorm<sq_T>* tmp=new mlnorm<sq_T> ( );
    01257                 tmp->set_rv ( rvn ); tmp->set_rvc ( rvc );
    01258                 tmp->set_parameters ( A,mu1-A*mu2,R_n );
    01259                 return tmp;
    01260         }
    01261
    01263
    01264         template<class sq_T>
    01265         std::ostream &operator<< ( std::ostream &os,  mlnorm<sq_T> &ml )
    01266         {
    01267                 os << "A:"<< ml.A<<endl;
    01268                 os << "mu:"<< ml.mu_const<<endl;
    01269                 os << "R:" << ml.epdf._R().to_mat() <<endl;
    01270                 return os;
    01271         };
    01272
    01273 }
    01274 #endif //EF_H