- Main Page
- Related Pages
- Classes
- Files

- File List
- File Members

# exp_family.h

Go to the documentation of this file.

    00001
    00013 #ifndef EF_H
    00014 #define EF_H
    00015
    00016
    00017 #include "../shared_ptr.h"
    00018 #include "../base/bdmbase.h"
    00019 #include "../math/chmat.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 {
    00052                                 double tmp;
    00053                                 tmp= evallog_nn ( val )-lognc();
    00054 //                              it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
    00055                                 return tmp;}
    00057                         virtual vec evallog ( const mat &Val ) const
    00058                         {
    00059                                 vec x ( Val.cols() );
    00060                                 for ( int i=0;i<Val.cols();i++ ) {x ( i ) =evallog_nn ( Val.get_col ( i ) ) ;}
    00061                                 return x-lognc();
    00062                         }
    00064                         virtual void pow ( double p ) {it_error ( "Not implemented" );};
    00065         };
    00066
    00073         class mEF : public mpdf
    00074         {
    00075
    00076                 public:
    00078                         mEF ( ) :mpdf ( ) {};
    00079         };
    00080
    00082         class BMEF : public BM
    00083         {
    00084                 protected:
    00086                         double frg;
    00088                         double last_lognc;
    00089                 public:
    00091                         BMEF ( double frg0=1.0 ) :BM ( ), frg ( frg0 ) {}
    00093                         BMEF ( const BMEF &B ) :BM ( B ), frg ( B.frg ), last_lognc ( B.last_lognc ) {}
    00095                         virtual void set_statistics ( const BMEF* BM0 ) {it_error ( "Not implemented" );};
    00097                         virtual void bayes ( const vec &data, const double w ) {};
    00098                         //original Bayes
    00099                         void bayes ( const vec &dt );
    00101                         virtual void flatten ( const BMEF * B ) {it_error ( "Not implemented" );}
    00103 //      virtual void flatten ( double nu0 ) {it_error ( "Not implemented" );}
    00104
    00105                         BMEF* _copy_ () const {it_error ( "function _copy_ not implemented for this BM" ); return NULL;};
    00106         };
    00107
    00108         template<class sq_T>
    00109         class mlnorm;
    00110
    00116         template<class sq_T>
    00117         class enorm : public eEF
    00118         {
    00119                 protected:
    00121                         vec mu;
    00123                         sq_T R;
    00124                 public:
    00127
    00128                         enorm ( ) :eEF ( ), mu ( ),R ( ) {};
    00129                         enorm ( const vec &mu,const sq_T &R ) {set_parameters ( mu,R );}
    00130                         void set_parameters ( const vec &mu,const sq_T &R );
    00131                         void from_setting(const Setting &root);
    00132                         void validate() {
    00133                                 it_assert(mu.length()==R.rows(),"parameters mismatch");
    00134                                 dim = mu.length();
    00135                         }
    00137
    00140
    00142                         void dupdate ( mat &v,double nu=1.0 );
    00143
    00144                         vec sample() const;
    00145
    00146                         double evallog_nn ( const vec &val ) const;
    00147                         double lognc () const;
    00148                         vec mean() const {return mu;}
    00149                         vec variance() const {return diag ( R.to_mat() );}
    00150 //      mlnorm<sq_T>* condition ( const RV &rvn ) const ; <=========== fails to cmpile. Why?
    00151                         mpdf* condition ( const RV &rvn ) const ;
    00152                         enorm<sq_T>* marginal ( const RV &rv ) const;
    00153 //                      epdf* marginal ( const RV &rv ) const;
    00155
    00158
    00159                         vec& _mu() {return mu;}
    00160                         void set_mu ( const vec mu0 ) { mu=mu0;}
    00161                         sq_T& _R() {return R;}
    00162                         const sq_T& _R() const {return R;}
    00164
    00165         };
    00166         UIREGISTER(enorm<chmat>);
    00167         UIREGISTER(enorm<ldmat>);
    00168         UIREGISTER(enorm<fsqmat>);
    00169
    00170
    00177         class egiw : public eEF
    00178         {
    00179                 protected:
    00181                         ldmat V;
    00183                         double nu;
    00185                         int dimx;
    00187                         int nPsi;
    00188                 public:
    00191                         egiw() :eEF() {};
    00192                         egiw ( int dimx0, ldmat V0, double nu0=-1.0 ) :eEF() {set_parameters ( dimx0,V0, nu0 );};
    00193
    00194                         void set_parameters ( int dimx0, ldmat V0, double nu0=-1.0 )
    00195                         {
    00196                                 dimx=dimx0;
    00197                                 nPsi = V0.rows()-dimx;
    00198                                 dim = dimx* ( dimx+nPsi ); // size(R) + size(Theta)
    00199
    00200                                 V=V0;
    00201                                 if ( nu0<0 )
    00202                                 {
    00203                                         nu = 0.1 +nPsi +2*dimx +2; // +2 assures finite expected value of R
    00204                                         // terms before that are sufficient for finite normalization
    00205                                 }
    00206                                 else
    00207                                 {
    00208                                         nu=nu0;
    00209                                 }
    00210                         }
    00212
    00213                         vec sample() const;
    00214                         vec mean() const;
    00215                         vec variance() const;
    00216
    00218                         vec est_theta() const;
    00219
    00221                         ldmat est_theta_cov() const;
    00222
    00223                         void mean_mat ( mat &M, mat&R ) const;
    00225                         double evallog_nn ( const vec &val ) const;
    00226                         double lognc () const;
    00227                         void pow ( double p ) {V*=p;nu*=p;};
    00228
    00231
    00232                         ldmat& _V() {return V;}
    00233                         const ldmat& _V() const {return V;}
    00234                         double& _nu()  {return nu;}
    00235                         const double& _nu() const {return nu;}
    00236                         void from_setting(const Setting &set){
    00237                                 UI::get(nu,set,"nu");
    00238                                 UI::get(dimx,set,"dimx");
    00239                                 set.lookupValue("nu",nu);
    00240                                 set.lookupValue("dimx",dimx);
    00241                                 mat V;
    00242                                 UI::get(V,set,"V", UI::compulsory);
    00243                                 set_parameters(dimx, V, nu);
    00244                                 RV* rv=UI::build<RV>(set,"rv", UI::compulsory);
    00245                                 set_rv(*rv);
    00246                                 delete rv;
    00247                         }
    00249         };
    00250         UIREGISTER(egiw);
    00251
    00260         class eDirich: public eEF
    00261         {
    00262                 protected:
    00264                         vec beta;
    00265                 public:
    00268
    00269                         eDirich () : eEF ( ) {};
    00270                         eDirich ( const eDirich &D0 ) : eEF () {set_parameters ( D0.beta );};
    00271                         eDirich ( const vec &beta0 ) {set_parameters ( beta0 );};
    00272                         void set_parameters ( const vec &beta0 )
    00273                         {
    00274                                 beta= beta0;
    00275                                 dim = beta.length();
    00276                         }
    00278
    00279                         vec sample() const {it_error ( "Not implemented" );return vec_1 ( 0.0 );};
    00280                         vec mean() const {return beta/sum(beta);};
    00281                         vec variance() const {double gamma =sum(beta); return elem_mult ( beta, ( beta+1 ) ) / ( gamma* ( gamma+1 ) );}
    00283                         double evallog_nn ( const vec &val ) const
    00284                         {
    00285                                 double tmp; tmp= ( beta-1 ) *log ( val );
    00286 //                              it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
    00287                                 return tmp;
    00288                         };
    00289                         double lognc () const
    00290                         {
    00291                                 double tmp;
    00292                                 double gam=sum ( beta );
    00293                                 double lgb=0.0;
    00294                                 for ( int i=0;i<beta.length();i++ ) {lgb+=lgamma ( beta ( i ) );}
    00295                                 tmp= lgb-lgamma ( gam );
    00296 //                              it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
    00297                                 return tmp;
    00298                         };
    00300                         vec& _beta()  {return beta;}
    00302         };
    00303
    00305         class multiBM : public BMEF
    00306         {
    00307                 protected:
    00309                         eDirich est;
    00311                         vec &beta;
    00312                 public:
    00314                         multiBM ( ) : BMEF ( ),est ( ),beta ( est._beta() )
    00315                         {
    00316                                 if ( beta.length() >0 ) {last_lognc=est.lognc();}
    00317                                 else{last_lognc=0.0;}
    00318                         }
    00320                         multiBM ( const multiBM &B ) : BMEF ( B ),est ( B.est ),beta ( est._beta() ) {}
    00322                         void set_statistics ( const BM* mB0 ) {const multiBM* mB=dynamic_cast<const multiBM*> ( mB0 ); beta=mB->beta;}
    00323                         void bayes ( const vec &dt )
    00324                         {
    00325                                 if ( frg<1.0 ) {beta*=frg;last_lognc=est.lognc();}
    00326                                 beta+=dt;
    00327                                 if ( evalll ) {ll=est.lognc()-last_lognc;}
    00328                         }
    00329                         double logpred ( const vec &dt ) const
    00330                         {
    00331                                 eDirich pred ( est );
    00332                                 vec &beta = pred._beta();
    00333
    00334                                 double lll;
    00335                                 if ( frg<1.0 )
    00336                                         {beta*=frg;lll=pred.lognc();}
    00337                                 else
    00338                                         if ( evalll ) {lll=last_lognc;}
    00339                                         else{lll=pred.lognc();}
    00340
    00341                                 beta+=dt;
    00342                                 return pred.lognc()-lll;
    00343                         }
    00344                         void flatten ( const BMEF* B )
    00345                         {
    00346                                 const multiBM* E=dynamic_cast<const multiBM*> ( B );
    00347                                 // sum(beta) should be equal to sum(B.beta)
    00348                                 const vec &Eb=E->beta;//const_cast<multiBM*> ( E )->_beta();
    00349                                 beta*= ( sum ( Eb ) /sum ( beta ) );
    00350                                 if ( evalll ) {last_lognc=est.lognc();}
    00351                         }
    00352                         const epdf& posterior() const {return est;};
    00353                         const eDirich* _e() const {return &est;};
    00354                         void set_parameters ( const vec &beta0 )
    00355                         {
    00356                                 est.set_parameters ( beta0 );
    00357                                 if ( evalll ) {last_lognc=est.lognc();}
    00358                         }
    00359         };
    00360
    00370         class egamma : public eEF
    00371         {
    00372                 protected:
    00374                         vec alpha;
    00376                         vec beta;
    00377                 public :
    00380                         egamma ( ) :eEF ( ), alpha ( 0 ), beta ( 0 ) {};
    00381                         egamma ( const vec &a, const vec &b ) {set_parameters ( a, b );};
    00382                         void set_parameters ( const vec &a, const vec &b ) {alpha=a,beta=b;dim = alpha.length();};
    00384
    00385                         vec sample() const;
    00387 //      mat sample ( int N ) const;
    00388                         double evallog ( const vec &val ) const;
    00389                         double lognc () const;
    00391                         vec& _alpha() {return alpha;}
    00392                         vec& _beta() {return beta;}
    00393                         vec mean() const {return elem_div ( alpha,beta );}
    00394                         vec variance() const {return elem_div ( alpha,elem_mult ( beta,beta ) ); }
    00395
    00404                         void from_setting(const Setting &set){
    00405                                 epdf::from_setting(set); // reads rv
    00406                                 UI::get(alpha,set,"alpha", UI::compulsory);
    00407                                 UI::get(beta,set,"beta", UI::compulsory);
    00408                                 validate();
    00409                         }
    00410                         void validate(){
    00411                                 it_assert(alpha.length() ==beta.length(), "parameters do not match");
    00412                                 dim =alpha.length();
    00413                         }
    00414         };
    00415 UIREGISTER(egamma);
    00432         class eigamma : public egamma
    00433         {
    00434                 protected:
    00435                 public :
    00440
    00441                         vec sample() const {return 1.0/egamma::sample();};
    00443                         vec mean() const {return elem_div ( beta,alpha-1 );}
    00444                         vec variance() const {vec mea=mean(); return elem_div ( elem_mult ( mea,mea ),alpha-2 );}
    00445         };
    00446         /*
    00448         class emix : public epdf {
    00449         protected:
    00450                 int n;
    00451                 vec &w;
    00452                 Array<epdf*> Coms;
    00453         public:
    00455                 emix ( const RV &rv, vec &w0): epdf(rv), n(w0.length()), w(w0), Coms(n) {};
    00456                 void set_parameters( int &i, double wi, epdf* ep){w(i)=wi;Coms(i)=ep;}
    00457                 vec mean(){vec pom; for(int i=0;i<n;i++){pom+=Coms(i)->mean()*w(i);} return pom;};
    00458                 vec sample() {it_error ( "Not implemented" );return 0;}
    00459         };
    00460         */
    00461
    00463
    00464         class euni: public epdf
    00465         {
    00466                 protected:
    00468                         vec low;
    00470                         vec high;
    00472                         vec distance;
    00474                         double nk;
    00476                         double lnk;
    00477                 public:
    00480                         euni ( ) :epdf ( ) {}
    00481                         euni ( const vec &low0, const vec &high0 ) {set_parameters ( low0,high0 );}
    00482                         void set_parameters ( const vec &low0, const vec &high0 )
    00483                         {
    00484                                 distance = high0-low0;
    00485                                 it_assert_debug ( min ( distance ) >0.0,"bad support" );
    00486                                 low = low0;
    00487                                 high = high0;
    00488                                 nk = prod ( 1.0/distance );
    00489                                 lnk = log ( nk );
    00490                                 dim = low.length();
    00491                         }
    00493
    00494                         double eval ( const vec &val ) const  {return nk;}
    00495                         double evallog ( const vec &val ) const  {
    00496                                 if (any(val<low) && any(val>high)) {return inf;}
    00497                                 else return lnk;
    00498                         }
    00499                         vec sample() const
    00500                         {
    00501                                 vec smp ( dim );
    00502 #pragma omp critical
    00503                                 UniRNG.sample_vector ( dim ,smp );
    00504                                 return low+elem_mult ( distance,smp );
    00505                         }
    00507                         vec mean() const {return ( high-low ) /2.0;}
    00508                         vec variance() const {return ( pow ( high,2 ) +pow ( low,2 ) +elem_mult ( high,low ) ) /3.0;}
    00517                         void from_setting(const Setting &set){
    00518                                 epdf::from_setting(set); // reads rv and rvc
    00519
    00520                                 UI::get(high,set,"high", UI::compulsory);
    00521                                 UI::get(low,set,"low", UI::compulsory);
    00522                         }
    00523         };
    00524
    00525
    00531         template<class sq_T>
    00532         class mlnorm : public mEF
    00533         {
    00534                 protected:
    00536                         shared_ptr<enorm<sq_T> > iepdf;
    00537                         mat A;
    00538                         vec mu_const;
    00539                         vec& _mu; //cached epdf.mu;
    00540                 public:
    00543                         mlnorm():iepdf(new enorm<sq_T>()), _mu(iepdf->_mu()) { set_ep(iepdf); };
    00544                         mlnorm (const mat &A, const vec &mu0, const sq_T &R ) :iepdf(new enorm<sq_T>()), _mu(iepdf->_mu())
    00545                         {
    00546                                 set_ep(iepdf); set_parameters ( A,mu0,R );
    00547                         }
    00548
    00550                         void set_parameters ( const  mat &A, const vec &mu0, const sq_T &R );
    00553                         void condition ( const vec &cond );
    00554
    00556                         vec& _mu_const() {return mu_const;}
    00558                         mat& _A() {return A;}
    00560                         mat _R() { return iepdf->_R().to_mat(); }
    00561
    00562                         template<class sq_M>
    00563                         friend std::ostream &operator<< ( std::ostream &os,  mlnorm<sq_M> &ml );
    00564
    00565                         void from_setting(const Setting &set){
    00566                                 mpdf::from_setting(set);
    00567
    00568                                 UI::get(A,set,"A", UI::compulsory);
    00569                                 UI::get(mu_const,set,"const", UI::compulsory);
    00570                                 mat R0;
    00571                                 UI::get(R0,set,"R", UI::compulsory);
    00572                                 set_parameters(A,mu_const,R0);
    00573                         };
    00574         };
    00575         UIREGISTER(mlnorm<ldmat>);
    00576         UIREGISTER(mlnorm<fsqmat>);
    00577         UIREGISTER(mlnorm<chmat>);
    00578
    00580         template<class sq_T>
    00581         class mgnorm : public mEF
    00582         {
    00583                 protected:
    00585                         shared_ptr<enorm<sq_T> > iepdf;
    00586                         vec &mu;
    00587                         fnc* g;
    00588                 public:
    00590                         mgnorm():iepdf(new enorm<sq_T>()), mu(iepdf->_mu()) { set_ep(iepdf); }
    00592                         void set_parameters ( fnc* g0, const sq_T &R0 ) {g=g0; iepdf->set_parameters ( zeros ( g->dimension() ), R0 );}
    00593                         void condition ( const vec &cond ) {mu=g->eval ( cond );};
    00594
    00595
    00623                         void from_setting( const Setting &set )
    00624                         {
    00625                                 fnc* g = UI::build<fnc>( set, "g", UI::compulsory );
    00626
    00627                                 mat R;
    00628                                 vec dR;
    00629                                 if ( UI::get( dR, set, "dR" ) )
    00630                                         R=diag(dR);
    00631                                 else
    00632                                         UI::get( R, set, "R", UI::compulsory);
    00633
    00634                                 set_parameters(g,R);
    00635                         }
    00636         };
    00637
    00638         UIREGISTER(mgnorm<chmat>);
    00639
    00640
    00648         class mlstudent : public mlnorm<ldmat>
    00649         {
    00650                 protected:
    00651                         ldmat Lambda;
    00652                         ldmat &_R;
    00653                         ldmat Re;
    00654                 public:
    00655                         mlstudent ( ) :mlnorm<ldmat> (),
    00656                                         Lambda (),      _R ( iepdf->_R() ) {}
    00657                         void set_parameters ( const mat &A0, const vec &mu0, const ldmat &R0, const ldmat& Lambda0 )
    00658                         {
    00659                                 it_assert_debug ( A0.rows() ==mu0.length(),"" );
    00660                                 it_assert_debug ( R0.rows() ==A0.rows(),"" );
    00661
    00662                                 iepdf->set_parameters ( mu0,Lambda ); //
    00663                                 A = A0;
    00664                                 mu_const = mu0;
    00665                                 Re=R0;
    00666                                 Lambda = Lambda0;
    00667                         }
    00668                         void condition ( const vec &cond )
    00669                         {
    00670                                 _mu = A*cond + mu_const;
    00671                                 double zeta;
    00672                                 //ugly hack!
    00673                                 if ( ( cond.length() +1 ) ==Lambda.rows() )
    00674                                 {
    00675                                         zeta = Lambda.invqform ( concat ( cond, vec_1 ( 1.0 ) ) );
    00676                                 }
    00677                                 else
    00678                                 {
    00679                                         zeta = Lambda.invqform ( cond );
    00680                                 }
    00681                                 _R = Re;
    00682                                 _R*= ( 1+zeta );// / ( nu ); << nu is in Re!!!!!!
    00683                         };
    00684
    00685         };
    00695         class mgamma : public mEF
    00696         {
    00697                 protected:
    00699                         shared_ptr<egamma> iepdf;
    00700
    00702                         double k;
    00703
    00705                         vec &_beta;
    00706
    00707                 public:
    00709                         mgamma():iepdf(new egamma()), k(0),
    00710                             _beta(iepdf->_beta()) {
    00711                             set_ep(iepdf);
    00712                         }
    00713
    00715                         void set_parameters(double k, const vec &beta0);
    00716
    00717                         void condition ( const vec &val ) {_beta=k/val;};
    00727                         void from_setting(const Setting &set){
    00728                                 mpdf::from_setting(set); // reads rv and rvc
    00729                                 vec betatmp; // ugly but necessary
    00730                                 UI::get(betatmp,set,"beta", UI::compulsory);
    00731                                 UI::get(k,set,"k", UI::compulsory);
    00732                                 set_parameters(k,betatmp);
    00733                         }
    00734         };
    00735         UIREGISTER(mgamma);
    00736
    00746         class migamma : public mEF
    00747         {
    00748                 protected:
    00750                         shared_ptr<eigamma> iepdf;
    00751
    00753                         double k;
    00754
    00756                         vec &_alpha;
    00757
    00759                         vec &_beta;
    00760
    00761                 public:
    00764                         migamma():iepdf(new eigamma()),
    00765                             k(0),
    00766                             _alpha(iepdf->_alpha()),
    00767                             _beta(iepdf->_beta()) {
    00768                             set_ep(iepdf);
    00769                         }
    00770
    00771                         migamma(const migamma &m):iepdf(m.iepdf),
    00772                             k(0),
    00773                             _alpha(iepdf->_alpha()),
    00774                             _beta(iepdf->_beta()) {
    00775                             set_ep(iepdf);
    00776                         }
    00778
    00780                         void set_parameters ( int len, double k0 )
    00781                         {
    00782                                 k=k0;
    00783                                 iepdf->set_parameters ( ( 1.0/ ( k*k ) +2.0 ) *ones ( len ) /*alpha*/, ones ( len ) /*beta*/ );
    00784                                 dimc = dimension();
    00785                         };
    00786                         void condition ( const vec &val )
    00787                         {
    00788                                 _beta=elem_mult ( val, ( _alpha-1.0 ) );
    00789                         };
    00790         };
    00791
    00792
    00804         class mgamma_fix : public mgamma
    00805         {
    00806                 protected:
    00808                         double l;
    00810                         vec refl;
    00811                 public:
    00813                         mgamma_fix ( ) : mgamma ( ),refl () {};
    00815                         void set_parameters ( double k0 , vec ref0, double l0 )
    00816                         {
    00817                                 mgamma::set_parameters ( k0, ref0 );
    00818                                 refl=pow ( ref0,1.0-l0 );l=l0;
    00819                                 dimc=dimension();
    00820                         };
    00821
    00822                         void condition ( const vec &val ) {vec mean=elem_mult ( refl,pow ( val,l ) ); _beta=k/mean;};
    00823         };
    00824
    00825
    00838         class migamma_ref : public migamma
    00839         {
    00840                 protected:
    00842                         double l;
    00844                         vec refl;
    00845                 public:
    00847                         migamma_ref ( ) : migamma (),refl ( ) {};
    00849                         void set_parameters ( double k0 , vec ref0, double l0 )
    00850                         {
    00851                                 migamma::set_parameters ( ref0.length(), k0 );
    00852                                 refl=pow ( ref0,1.0-l0 );
    00853                                 l=l0;
    00854                                 dimc = dimension();
    00855                         };
    00856
    00857                         void condition ( const vec &val )
    00858                         {
    00859                                 vec mean=elem_mult ( refl,pow ( val,l ) );
    00860                                 migamma::condition ( mean );
    00861                         };
    00862
    00883                         void from_setting( const Setting &set );
    00884
    00885                         // TODO dodelat void to_setting( Setting &set ) const;
    00886         };
    00887
    00888
    00889         UIREGISTER(migamma_ref);
    00890
    00900         class elognorm: public enorm<ldmat>
    00901         {
    00902                 public:
    00903                         vec sample() const {return exp ( enorm<ldmat>::sample() );};
    00904                         vec mean() const {vec var=enorm<ldmat>::variance();return exp ( mu - 0.5*var );};
    00905
    00906         };
    00907
    00919         class mlognorm : public mpdf
    00920         {
    00921                 protected:
    00922                         shared_ptr<elognorm> eno;
    00923
    00925                         double sig2;
    00926
    00928                         vec &mu;
    00929                 public:
    00931                         mlognorm():eno(new elognorm()),
    00932                             sig2(0),
    00933                             mu(eno->_mu()) {
    00934                             set_ep(eno);
    00935                         }
    00936
    00938                         void set_parameters ( int size, double k )
    00939                         {
    00940                                 sig2 = 0.5*log ( k*k+1 );
    00941                                 eno->set_parameters ( zeros ( size ),2*sig2*eye ( size ) );
    00942
    00943                                 dimc = size;
    00944                         };
    00945
    00946                         void condition ( const vec &val )
    00947                         {
    00948                                 mu=log ( val )-sig2;//elem_mult ( refl,pow ( val,l ) );
    00949                         };
    00950
    00969                         void from_setting( const Setting &set );
    00970
    00971                         // TODO dodelat void to_setting( Setting &set ) const;
    00972
    00973         };
    00974
    00975         UIREGISTER(mlognorm);
    00976
    00980         class eWishartCh : public epdf
    00981         {
    00982                 protected:
    00984                         chmat Y;
    00986                         int p;
    00988                         double delta;
    00989                 public:
    00990                         void set_parameters ( const mat &Y0, const double delta0 ) {Y=chmat ( Y0 );delta=delta0; p=Y.rows(); dim = p*p; }
    00991                         mat sample_mat() const
    00992                         {
    00993                                 mat X=zeros ( p,p );
    00994
    00995                                 //sample diagonal
    00996                                 for ( int i=0;i<p;i++ )
    00997                                 {
    00998                                         GamRNG.setup ( 0.5* ( delta-i ) , 0.5 ); // no +1 !! index if from 0
    00999 #pragma omp critical
    01000                                         X ( i,i ) =sqrt ( GamRNG() );
    01001                                 }
    01002                                 //do the rest
    01003                                 for ( int i=0;i<p;i++ )
    01004                                 {
    01005                                         for ( int j=i+1;j<p;j++ )
    01006                                         {
    01007 #pragma omp critical
    01008                                                 X ( i,j ) =NorRNG.sample();
    01009                                         }
    01010                                 }
    01011                                 return X*Y._Ch();// return upper triangular part of the decomposition
    01012                         }
    01013                         vec sample () const
    01014                         {
    01015                                 return vec ( sample_mat()._data(),p*p );
    01016                         }
    01018                         void setY ( const mat &Ch0 ) {copy_vector ( dim,Ch0._data(), Y._Ch()._data() );}
    01020                         void _setY ( const vec &ch0 ) {copy_vector ( dim, ch0._data(), Y._Ch()._data() ); }
    01022                         const chmat& getY()const {return Y;}
    01023         };
    01024
    01025         class eiWishartCh: public epdf
    01026         {
    01027                 protected:
    01028                         eWishartCh W;
    01029                         int p;
    01030                         double delta;
    01031                 public:
    01032                         void set_parameters ( const mat &Y0, const double delta0) {
    01033                                 delta = delta0;
    01034                                 W.set_parameters ( inv ( Y0 ),delta0 );
    01035                                 dim = W.dimension(); p=Y0.rows();
    01036                         }
    01037                         vec sample() const {mat iCh; iCh=inv ( W.sample_mat() ); return vec ( iCh._data(),dim );}
    01038                         void _setY ( const vec &y0 )
    01039                         {
    01040                                 mat Ch ( p,p );
    01041                                 mat iCh ( p,p );
    01042                                 copy_vector ( dim, y0._data(), Ch._data() );
    01043
    01044                                 iCh=inv ( Ch );
    01045                                 W.setY ( iCh );
    01046                         }
    01047                         virtual double evallog ( const vec &val ) const {
    01048                                 chmat X(p);
    01049                                 const chmat& Y=W.getY();
    01050
    01051                                 copy_vector(p*p,val._data(),X._Ch()._data());
    01052                                 chmat iX(p);X.inv(iX);
    01053                                 // compute
    01054 //                              \frac{ |\Psi|^{m/2}|X|^{-(m+p+1)/2}e^{-tr(\Psi X^{-1})/2} }{ 2^{mp/2}\Gamma_p(m/2)},
    01055                                 mat M=Y.to_mat()*iX.to_mat();
    01056
    01057                                 double log1 = 0.5*p*(2*Y.logdet())-0.5*(delta+p+1)*(2*X.logdet())-0.5*trace(M);
    01058                                 //Fixme! Multivariate gamma omitted!! it is ok for sampling, but not otherwise!!
    01059
    01060 /*                              if (0) {
    01061                                         mat XX=X.to_mat();
    01062                                         mat YY=Y.to_mat();
    01063
    01064                                         double log2 = 0.5*p*log(det(YY))-0.5*(delta+p+1)*log(det(XX))-0.5*trace(YY*inv(XX));
    01065                                         cout << log1 << "," << log2 << endl;
    01066                                 }*/
    01067                                 return log1;
    01068                         };
    01069
    01070         };
    01071
    01072         class rwiWishartCh : public mpdf
    01073         {
    01074                 protected:
    01075                         shared_ptr<eiWishartCh> eiW;
    01077                         double sqd;
    01078                         //reference point for diagonal
    01079                         vec refl;
    01080                         double l;
    01081                         int p;
    01082
    01083                 public:
    01084                         rwiWishartCh():eiW(new eiWishartCh()),
    01085                             sqd(0), l(0), p(0) {
    01086                             set_ep(eiW);
    01087                         }
    01088
    01089                         void set_parameters ( int p0, double k, vec ref0, double l0  )
    01090                         {
    01091                                 p=p0;
    01092                                 double delta = 2/(k*k)+p+3;
    01093                                 sqd=sqrt ( delta-p-1 );
    01094                                 l=l0;
    01095                                 refl=pow(ref0,1-l);
    01096
    01097                                 eiW->set_parameters ( eye ( p ),delta );
    01098                                 dimc=eiW->dimension();
    01099                         }
    01100                         void condition ( const vec &c ) {
    01101                                 vec z=c;
    01102                                 int ri=0;
    01103                                 for(int i=0;i<p*p;i+=(p+1)){//trace diagonal element
    01104                                         z(i) = pow(z(i),l)*refl(ri);
    01105                                         ri++;
    01106                                 }
    01107
    01108                                 eiW->_setY ( sqd*z );
    01109                         }
    01110         };
    01111
    01113         enum RESAMPLING_METHOD { MULTINOMIAL = 0, STRATIFIED = 1, SYSTEMATIC = 3 };
    01119         class eEmp: public epdf
    01120         {
    01121                 protected :
    01123                         int n;
    01125                         vec w;
    01127                         Array<vec> samples;
    01128                 public:
    01131                         eEmp ( ) :epdf ( ),w ( ),samples ( ) {};
    01133                         eEmp ( const eEmp &e ) : epdf ( e ), w ( e.w ), samples ( e.samples ) {};
    01135
    01137                         void set_statistics ( const vec &w0, const epdf* pdf0 );
    01139                         void set_statistics ( const epdf* pdf0 , int n ) {set_statistics ( ones ( n ) /n,pdf0 );};
    01141                         void set_samples ( const epdf* pdf0 );
    01143                         void set_parameters ( int n0, bool copy=true ) {n=n0; w.set_size ( n0,copy );samples.set_size ( n0,copy );};
    01145                         vec& _w()  {return w;};
    01147                         const vec& _w() const {return w;};
    01149                         Array<vec>& _samples() {return samples;};
    01151                         const Array<vec>& _samples() const {return samples;};
    01153                         ivec resample ( RESAMPLING_METHOD method=SYSTEMATIC );
    01155                         vec sample() const {it_error ( "Not implemented" );return 0;}
    01157                         double evallog ( const vec &val ) const {it_error ( "Not implemented" );return 0.0;}
    01158                         vec mean() const
    01159                         {
    01160                                 vec pom=zeros ( dim );
    01161                                 for ( int i=0;i<n;i++ ) {pom+=samples ( i ) *w ( i );}
    01162                                 return pom;
    01163                         }
    01164                         vec variance() const
    01165                         {
    01166                                 vec pom=zeros ( dim );
    01167                                 for ( int i=0;i<n;i++ ) {pom+=pow ( samples ( i ),2 ) *w ( i );}
    01168                                 return pom-pow ( mean(),2 );
    01169                         }
    01171                         void qbounds ( vec &lb, vec &ub, double perc=0.95 ) const
    01172                         {
    01173                                 // lb in inf so than it will be pushed below;
    01174                                 lb.set_size ( dim );
    01175                                 ub.set_size ( dim );
    01176                                 lb = std::numeric_limits<double>::infinity();
    01177                                 ub = -std::numeric_limits<double>::infinity();
    01178                                 int j;
    01179                                 for ( int i=0;i<n;i++ )
    01180                                 {
    01181                                         for ( j=0;j<dim; j++ )
    01182                                         {
    01183                                                 if ( samples ( i ) ( j ) <lb ( j ) ) {lb ( j ) =samples ( i ) ( j );}
    01184                                                 if ( samples ( i ) ( j ) >ub ( j ) ) {ub ( j ) =samples ( i ) ( j );}
    01185                                         }
    01186                                 }
    01187                         }
    01188         };
    01189
    01190
    01192
    01193         template<class sq_T>
    01194         void enorm<sq_T>::set_parameters ( const vec &mu0, const sq_T &R0 )
    01195         {
    01196 //Fixme test dimensions of mu0 and R0;
    01197                 mu = mu0;
    01198                 R = R0;
    01199                 validate();
    01200         };
    01201
    01202         template<class sq_T>
    01203         void enorm<sq_T>::from_setting(const Setting &set){
    01204                 epdf::from_setting(set); //reads rv
    01205
    01206                 UI::get(mu,set,"mu", UI::compulsory);
    01207                 mat Rtmp;// necessary for conversion
    01208                 UI::get(Rtmp,set,"R", UI::compulsory);
    01209                 R=Rtmp; // conversion
    01210                 validate();
    01211         }
    01212
    01213         template<class sq_T>
    01214         void enorm<sq_T>::dupdate ( mat &v, double nu )
    01215         {
    01216                 //
    01217         };
    01218
    01219 // template<class sq_T>
    01220 // void enorm<sq_T>::tupdate ( double phi, mat &vbar, double nubar ) {
    01221 //      //
    01222 // };
    01223
    01224         template<class sq_T>
    01225         vec enorm<sq_T>::sample() const
    01226         {
    01227                 vec x ( dim );
    01228 #pragma omp critical
    01229                 NorRNG.sample_vector ( dim,x );
    01230                 vec smp = R.sqrt_mult ( x );
    01231
    01232                 smp += mu;
    01233                 return smp;
    01234         };
    01235
    01236 // template<class sq_T>
    01237 // double enorm<sq_T>::eval ( const vec &val ) const {
    01238 //      double pdfl,e;
    01239 //      pdfl = evallog ( val );
    01240 //      e = exp ( pdfl );
    01241 //      return e;
    01242 // };
    01243
    01244         template<class sq_T>
    01245         double enorm<sq_T>::evallog_nn ( const vec &val ) const
    01246         {
    01247                 // 1.83787706640935 = log(2pi)
    01248                 double tmp=-0.5* ( R.invqform ( mu-val ) );// - lognc();
    01249                 return  tmp;
    01250         };
    01251
    01252         template<class sq_T>
    01253         inline double enorm<sq_T>::lognc () const
    01254         {
    01255                 // 1.83787706640935 = log(2pi)
    01256                 double tmp=0.5* ( R.cols() * 1.83787706640935 +R.logdet() );
    01257                 return tmp;
    01258         };
    01259
    01260         template<class sq_T>
    01261         void mlnorm<sq_T>::set_parameters ( const mat &A0, const vec &mu0, const sq_T &R0 )
    01262         {
    01263                 it_assert_debug ( A0.rows() ==mu0.length(),"" );
    01264                 it_assert_debug ( A0.rows() ==R0.rows(),"" );
    01265
    01266                 iepdf->set_parameters(zeros(A0.rows()), R0);
    01267                 A = A0;
    01268                 mu_const = mu0;
    01269                 dimc = A0.cols();
    01270         }
    01271
    01272 // template<class sq_T>
    01273 // vec mlnorm<sq_T>::samplecond (const  vec &cond, double &lik ) {
    01274 //      this->condition ( cond );
    01275 //      vec smp = epdf.sample();
    01276 //      lik = epdf.eval ( smp );
    01277 //      return smp;
    01278 // }
    01279
    01280 // template<class sq_T>
    01281 // mat mlnorm<sq_T>::samplecond (const vec &cond, vec &lik, int n ) {
    01282 //      int i;
    01283 //      int dim = rv.count();
    01284 //      mat Smp ( dim,n );
    01285 //      vec smp ( dim );
    01286 //      this->condition ( cond );
    01287 //
    01288 //      for ( i=0; i<n; i++ ) {
    01289 //              smp = epdf.sample();
    01290 //              lik ( i ) = epdf.eval ( smp );
    01291 //              Smp.set_col ( i ,smp );
    01292 //      }
    01293 //
    01294 //      return Smp;
    01295 // }
    01296
    01297         template<class sq_T>
    01298         void mlnorm<sq_T>::condition ( const vec &cond )
    01299         {
    01300                 _mu = A*cond + mu_const;
    01301 //R is already assigned;
    01302         }
    01303
    01304         template<class sq_T>
    01305         enorm<sq_T>* enorm<sq_T>::marginal ( const RV &rvn ) const
    01306         {
    01307                 it_assert_debug ( isnamed(), "rv description is not assigned" );
    01308                 ivec irvn = rvn.dataind ( rv );
    01309
    01310                 sq_T Rn ( R,irvn ); //select rows and columns of R
    01311
    01312                 enorm<sq_T>* tmp = new enorm<sq_T>;
    01313                 tmp->set_rv ( rvn );
    01314                 tmp->set_parameters ( mu ( irvn ), Rn );
    01315                 return tmp;
    01316         }
    01317
    01318         template<class sq_T>
    01319         mpdf* enorm<sq_T>::condition ( const RV &rvn ) const
    01320         {
    01321
    01322                 it_assert_debug ( isnamed(),"rvs are not assigned" );
    01323
    01324                 RV rvc = rv.subt ( rvn );
    01325                 it_assert_debug ( ( rvc._dsize() +rvn._dsize() ==rv._dsize() ),"wrong rvn" );
    01326                 //Permutation vector of the new R
    01327                 ivec irvn = rvn.dataind ( rv );
    01328                 ivec irvc = rvc.dataind ( rv );
    01329                 ivec perm=concat ( irvn , irvc );
    01330                 sq_T Rn ( R,perm );
    01331
    01332                 //fixme - could this be done in general for all sq_T?
    01333                 mat S=Rn.to_mat();
    01334                 //fixme
    01335                 int n=rvn._dsize()-1;
    01336                 int end=R.rows()-1;
    01337                 mat S11 = S.get ( 0,n, 0, n );
    01338                 mat S12 = S.get ( 0, n , rvn._dsize(), end );
    01339                 mat S22 = S.get ( rvn._dsize(), end, rvn._dsize(), end );
    01340
    01341                 vec mu1 = mu ( irvn );
    01342                 vec mu2 = mu ( irvc );
    01343                 mat A=S12*inv ( S22 );
    01344                 sq_T R_n ( S11 - A *S12.T() );
    01345
    01346                 mlnorm<sq_T>* tmp=new mlnorm<sq_T> ( );
    01347                 tmp->set_rv ( rvn ); tmp->set_rvc ( rvc );
    01348                 tmp->set_parameters ( A,mu1-A*mu2,R_n );
    01349                 return tmp;
    01350         }
    01351
    01353
    01354         template<class sq_T>
    01355         std::ostream &operator<< ( std::ostream &os,  mlnorm<sq_T> &ml )
    01356         {
    01357                 os << "A:"<< ml.A<<endl;
    01358                 os << "mu:"<< ml.mu_const<<endl;
    01359                 os << "R:" << ml.iepdf->_R().to_mat() <<endl;
    01360                 return os;
    01361         };
    01362
    01363 }
    01364 #endif //EF_H