- Main Page
- Related Pages
- Modules
- 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 "../base/bdmbase.h"
    00018 #include "../math/chmat.h"
    00019
    00020 namespace bdm
    00021 {
    00022
    00023
    00025         extern Uniform_RNG UniRNG;
    00027         extern Normal_RNG NorRNG;
    00029         extern Gamma_RNG GamRNG;
    00030
    00037         class eEF : public epdf
    00038         {
    00039                 public:
    00040 //      eEF() :epdf() {};
    00042                         eEF ( ) :epdf ( ) {};
    00044                         virtual double lognc() const =0;
    00046                         virtual void dupdate ( mat &v ) {it_error ( "Not implemented" );};
    00048                         virtual double evallog_nn ( const vec &val ) const{it_error ( "Not implemented" );return 0.0;};
    00050                         virtual double evallog ( const vec &val ) const {
    00051                                 double tmp;
    00052                                 tmp= evallog_nn ( val )-lognc();
    00053                                 it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
    00054                                 return tmp;}
    00056                         virtual vec evallog ( const mat &Val ) const
    00057                         {
    00058                                 vec x ( Val.cols() );
    00059                                 for ( int i=0;i<Val.cols();i++ ) {x ( i ) =evallog_nn ( Val.get_col ( i ) ) ;}
    00060                                 return x-lognc();
    00061                         }
    00063                         virtual void pow ( double p ) {it_error ( "Not implemented" );};
    00064         };
    00065
    00072         class mEF : public mpdf
    00073         {
    00074
    00075                 public:
    00077                         mEF ( ) :mpdf ( ) {};
    00078         };
    00079
    00081         class BMEF : public BM
    00082         {
    00083                 protected:
    00085                         double frg;
    00087                         double last_lognc;
    00088                 public:
    00090                         BMEF ( double frg0=1.0 ) :BM ( ), frg ( frg0 ) {}
    00092                         BMEF ( const BMEF &B ) :BM ( B ), frg ( B.frg ), last_lognc ( B.last_lognc ) {}
    00094                         virtual void set_statistics ( const BMEF* BM0 ) {it_error ( "Not implemented" );};
    00096                         virtual void bayes ( const vec &data, const double w ) {};
    00097                         //original Bayes
    00098                         void bayes ( const vec &dt );
    00100                         virtual void flatten ( const BMEF * B ) {it_error ( "Not implemented" );}
    00102 //      virtual void flatten ( double nu0 ) {it_error ( "Not implemented" );}
    00103
    00104                         BMEF* _copy_ () const {it_error ( "function _copy_ not implemented for this BM" ); return NULL;};
    00105         };
    00106
    00107         template<class sq_T>
    00108         class mlnorm;
    00109
    00115         template<class sq_T>
    00116         class enorm : public eEF
    00117         {
    00118                 protected:
    00120                         vec mu;
    00122                         sq_T R;
    00123                 public:
    00126
    00127                         enorm ( ) :eEF ( ), mu ( ),R ( ) {};
    00128                         enorm ( const vec &mu,const sq_T &R ) {set_parameters ( mu,R );}
    00129                         void set_parameters ( const vec &mu,const sq_T &R );
    00130                         void from_setting(const Setting &root);
    00131                         void validate() {
    00132                                 it_assert(mu.length()==R.rows(),"parameters mismatch");
    00133                                 dim = mu.length();
    00134                         }
    00136
    00139
    00141                         void dupdate ( mat &v,double nu=1.0 );
    00142
    00143                         vec sample() const;
    00144                         mat sample ( int N ) const;
    00145                         double evallog_nn ( const vec &val ) const;
    00146                         double lognc () const;
    00147                         vec mean() const {return mu;}
    00148                         vec variance() const {return diag ( R.to_mat() );}
    00149 //      mlnorm<sq_T>* condition ( const RV &rvn ) const ; <=========== fails to cmpile. Why?
    00150                         mpdf* condition ( const RV &rvn ) const ;
    00151                         enorm<sq_T>* marginal ( const RV &rv ) const;
    00152 //                      epdf* marginal ( const RV &rv ) const;
    00154
    00157
    00158                         vec& _mu() {return mu;}
    00159                         void set_mu ( const vec mu0 ) { mu=mu0;}
    00160                         sq_T& _R() {return R;}
    00161                         const sq_T& _R() const {return R;}
    00163
    00164         };
    00165         UIREGISTER(enorm<chmat>);
    00166         UIREGISTER(enorm<ldmat>);
    00167         UIREGISTER(enorm<fsqmat>);
    00168
    00169
    00176         class egiw : public eEF
    00177         {
    00178                 protected:
    00180                         ldmat V;
    00182                         double nu;
    00184                         int dimx;
    00186                         int nPsi;
    00187                 public:
    00190                         egiw() :eEF() {};
    00191                         egiw ( int dimx0, ldmat V0, double nu0=-1.0 ) :eEF() {set_parameters ( dimx0,V0, nu0 );};
    00192
    00193                         void set_parameters ( int dimx0, ldmat V0, double nu0=-1.0 )
    00194                         {
    00195                                 dimx=dimx0;
    00196                                 nPsi = V0.rows()-dimx;
    00197                                 dim = dimx* ( dimx+nPsi ); // size(R) + size(Theta)
    00198
    00199                                 V=V0;
    00200                                 if ( nu0<0 )
    00201                                 {
    00202                                         nu = 0.1 +nPsi +2*dimx +2; // +2 assures finite expected value of R
    00203                                         // terms before that are sufficient for finite normalization
    00204                                 }
    00205                                 else
    00206                                 {
    00207                                         nu=nu0;
    00208                                 }
    00209                         }
    00211
    00212                         vec sample() const;
    00213                         vec mean() const;
    00214                         vec variance() const;
    00215
    00217                         vec est_theta() const;
    00218
    00220                         ldmat est_theta_cov() const;
    00221
    00222                         void mean_mat ( mat &M, mat&R ) const;
    00224                         double evallog_nn ( const vec &val ) const;
    00225                         double lognc () const;
    00226                         void pow ( double p ) {V*=p;nu*=p;};
    00227
    00230
    00231                         ldmat& _V() {return V;}
    00232                         const ldmat& _V() const {return V;}
    00233                         double& _nu()  {return nu;}
    00234                         const double& _nu() const {return nu;}
    00235                         void from_setting(const Setting &set){
    00236                                 set.lookupValue("nu",nu);
    00237                                 set.lookupValue("dimx",dimx);
    00238                                 mat V;
    00239                                 UI::get(V,set,"V");
    00240                                 set_parameters(dimx, V, nu);
    00241                                 RV* rv=UI::build<RV>(set,"rv");
    00242                                 set_rv(*rv);
    00243                                 delete rv;
    00244                         }
    00246         };
    00247         UIREGISTER(egiw);
    00248
    00257         class eDirich: public eEF
    00258         {
    00259                 protected:
    00261                         vec beta;
    00262                 public:
    00265
    00266                         eDirich () : eEF ( ) {};
    00267                         eDirich ( const eDirich &D0 ) : eEF () {set_parameters ( D0.beta );};
    00268                         eDirich ( const vec &beta0 ) {set_parameters ( beta0 );};
    00269                         void set_parameters ( const vec &beta0 )
    00270                         {
    00271                                 beta= beta0;
    00272                                 dim = beta.length();
    00273                         }
    00275
    00276                         vec sample() const {it_error ( "Not implemented" );return vec_1 ( 0.0 );};
    00277                         vec mean() const {return beta/sum(beta);};
    00278                         vec variance() const {double gamma =sum(beta); return elem_mult ( beta, ( beta+1 ) ) / ( gamma* ( gamma+1 ) );}
    00280                         double evallog_nn ( const vec &val ) const
    00281                         {
    00282                                 double tmp; tmp= ( beta-1 ) *log ( val );               it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
    00283                                 return tmp;
    00284                         };
    00285                         double lognc () const
    00286                         {
    00287                                 double tmp;
    00288                                 double gam=sum ( beta );
    00289                                 double lgb=0.0;
    00290                                 for ( int i=0;i<beta.length();i++ ) {lgb+=lgamma ( beta ( i ) );}
    00291                                 tmp= lgb-lgamma ( gam );
    00292                                 it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
    00293                                 return tmp;
    00294                         };
    00296                         vec& _beta()  {return beta;}
    00298         };
    00299
    00301         class multiBM : public BMEF
    00302         {
    00303                 protected:
    00305                         eDirich est;
    00307                         vec &beta;
    00308                 public:
    00310                         multiBM ( ) : BMEF ( ),est ( ),beta ( est._beta() )
    00311                         {
    00312                                 if ( beta.length() >0 ) {last_lognc=est.lognc();}
    00313                                 else{last_lognc=0.0;}
    00314                         }
    00316                         multiBM ( const multiBM &B ) : BMEF ( B ),est ( B.est ),beta ( est._beta() ) {}
    00318                         void set_statistics ( const BM* mB0 ) {const multiBM* mB=dynamic_cast<const multiBM*> ( mB0 ); beta=mB->beta;}
    00319                         void bayes ( const vec &dt )
    00320                         {
    00321                                 if ( frg<1.0 ) {beta*=frg;last_lognc=est.lognc();}
    00322                                 beta+=dt;
    00323                                 if ( evalll ) {ll=est.lognc()-last_lognc;}
    00324                         }
    00325                         double logpred ( const vec &dt ) const
    00326                         {
    00327                                 eDirich pred ( est );
    00328                                 vec &beta = pred._beta();
    00329
    00330                                 double lll;
    00331                                 if ( frg<1.0 )
    00332                                         {beta*=frg;lll=pred.lognc();}
    00333                                 else
    00334                                         if ( evalll ) {lll=last_lognc;}
    00335                                         else{lll=pred.lognc();}
    00336
    00337                                 beta+=dt;
    00338                                 return pred.lognc()-lll;
    00339                         }
    00340                         void flatten ( const BMEF* B )
    00341                         {
    00342                                 const multiBM* E=dynamic_cast<const multiBM*> ( B );
    00343                                 // sum(beta) should be equal to sum(B.beta)
    00344                                 const vec &Eb=E->beta;//const_cast<multiBM*> ( E )->_beta();
    00345                                 beta*= ( sum ( Eb ) /sum ( beta ) );
    00346                                 if ( evalll ) {last_lognc=est.lognc();}
    00347                         }
    00348                         const epdf& posterior() const {return est;};
    00349                         const eDirich* _e() const {return &est;};
    00350                         void set_parameters ( const vec &beta0 )
    00351                         {
    00352                                 est.set_parameters ( beta0 );
    00353                                 if ( evalll ) {last_lognc=est.lognc();}
    00354                         }
    00355         };
    00356
    00366         class egamma : public eEF
    00367         {
    00368                 protected:
    00370                         vec alpha;
    00372                         vec beta;
    00373                 public :
    00376                         egamma ( ) :eEF ( ), alpha ( 0 ), beta ( 0 ) {};
    00377                         egamma ( const vec &a, const vec &b ) {set_parameters ( a, b );};
    00378                         void set_parameters ( const vec &a, const vec &b ) {alpha=a,beta=b;dim = alpha.length();};
    00380
    00381                         vec sample() const;
    00383 //      mat sample ( int N ) const;
    00384                         double evallog ( const vec &val ) const;
    00385                         double lognc () const;
    00387                         vec& _alpha() {return alpha;}
    00388                         vec& _beta() {return beta;}
    00389                         vec mean() const {return elem_div ( alpha,beta );}
    00390                         vec variance() const {return elem_div ( alpha,elem_mult ( beta,beta ) ); }
    00391
    00400                         void from_setting(const Setting &set){
    00401                                 epdf::from_setting(set); // reads rv
    00402                                 UI::get(alpha,set,"alpha");
    00403                                 UI::get(beta,set,"beta");
    00404                                 validate();
    00405                         }
    00406                         void validate(){
    00407                                 it_assert(alpha.length() ==beta.length(), "parameters do not match");
    00408                                 dim =alpha.length();
    00409                         }
    00410         };
    00411 UIREGISTER(egamma);
    00428         class eigamma : public egamma
    00429         {
    00430                 protected:
    00431                 public :
    00436
    00437                         vec sample() const {return 1.0/egamma::sample();};
    00439                         vec mean() const {return elem_div ( beta,alpha-1 );}
    00440                         vec variance() const {vec mea=mean(); return elem_div ( elem_mult ( mea,mea ),alpha-2 );}
    00441         };
    00442         /*
    00444         class emix : public epdf {
    00445         protected:
    00446                 int n;
    00447                 vec &w;
    00448                 Array<epdf*> Coms;
    00449         public:
    00451                 emix ( const RV &rv, vec &w0): epdf(rv), n(w0.length()), w(w0), Coms(n) {};
    00452                 void set_parameters( int &i, double wi, epdf* ep){w(i)=wi;Coms(i)=ep;}
    00453                 vec mean(){vec pom; for(int i=0;i<n;i++){pom+=Coms(i)->mean()*w(i);} return pom;};
    00454                 vec sample() {it_error ( "Not implemented" );return 0;}
    00455         };
    00456         */
    00457
    00459
    00460         class euni: public epdf
    00461         {
    00462                 protected:
    00464                         vec low;
    00466                         vec high;
    00468                         vec distance;
    00470                         double nk;
    00472                         double lnk;
    00473                 public:
    00476                         euni ( ) :epdf ( ) {}
    00477                         euni ( const vec &low0, const vec &high0 ) {set_parameters ( low0,high0 );}
    00478                         void set_parameters ( const vec &low0, const vec &high0 )
    00479                         {
    00480                                 distance = high0-low0;
    00481                                 it_assert_debug ( min ( distance ) >0.0,"bad support" );
    00482                                 low = low0;
    00483                                 high = high0;
    00484                                 nk = prod ( 1.0/distance );
    00485                                 lnk = log ( nk );
    00486                                 dim = low.length();
    00487                         }
    00489
    00490                         double eval ( const vec &val ) const  {return nk;}
    00491                         double evallog ( const vec &val ) const  {return lnk;}
    00492                         vec sample() const
    00493                         {
    00494                                 vec smp ( dim );
    00495 #pragma omp critical
    00496                                 UniRNG.sample_vector ( dim ,smp );
    00497                                 return low+elem_mult ( distance,smp );
    00498                         }
    00500                         vec mean() const {return ( high-low ) /2.0;}
    00501                         vec variance() const {return ( pow ( high,2 ) +pow ( low,2 ) +elem_mult ( high,low ) ) /3.0;}
    00510                         void from_setting(const Setting &set){
    00511                                 epdf::from_setting(set); // reads rv and rvc
    00512                                 UI::get(high,set,"high");
    00513                                 UI::get(low,set,"low");
    00514                         }
    00515         };
    00516
    00517
    00523         template<class sq_T>
    00524         class mlnorm : public mEF
    00525         {
    00526                 protected:
    00528                         enorm<sq_T> epdf;
    00529                         mat A;
    00530                         vec mu_const;
    00531                         vec& _mu; //cached epdf.mu;
    00532                 public:
    00535                         mlnorm ( ) :mEF (),epdf ( ),A ( ),_mu ( epdf._mu() ) {ep =&epdf; };
    00536                         mlnorm ( const  mat &A, const vec &mu0, const sq_T &R ) :epdf ( ),_mu ( epdf._mu() )
    00537                         {
    00538                                 ep =&epdf; set_parameters ( A,mu0,R );
    00539                         };
    00541                         void set_parameters ( const  mat &A, const vec &mu0, const sq_T &R );
    00544                         void condition ( const vec &cond );
    00545
    00547                         vec& _mu_const() {return mu_const;}
    00549                         mat& _A() {return A;}
    00551                         mat _R() {return epdf._R().to_mat();}
    00552
    00553                         template<class sq_M>
    00554                         friend std::ostream &operator<< ( std::ostream &os,  mlnorm<sq_M> &ml );
    00555         };
    00556
    00558         template<class sq_T>
    00559         class mgnorm : public mEF
    00560         {
    00561                 protected:
    00563                         enorm<sq_T> epdf;
    00564                         vec &mu;
    00565                         fnc* g;
    00566                 public:
    00568                         mgnorm() :mu ( epdf._mu() ) {ep=&epdf;}
    00570                         void set_parameters ( fnc* g0, const sq_T &R0 ) {g=g0; epdf.set_parameters ( zeros ( g->dimension() ), R0 );}
    00571                         void condition ( const vec &cond ) {mu=g->eval ( cond );};
    00572
    00573
    00601                         void from_setting( const Setting &set )
    00602                         {
    00603                                 fnc* g = UI::build<fnc>( set, "g" );
    00604
    00605                                 mat R;
    00606                                 if ( set.exists( "dR" ) )
    00607                                 {
    00608                                         vec dR;
    00609                                         UI::get( dR, set, "dR" );
    00610                                         R=diag(dR);
    00611                                 }
    00612                                 else
    00613                                         UI::get( R, set, "R");
    00614
    00615                                 set_parameters(g,R);
    00616                         }
    00617
    00618                         /*void mgnorm::to_setting( Setting &set ) const
    00619                         {
    00620                                 Transport::to_setting( set );
    00621
    00622                                 Setting &kilometers_setting = set.add("kilometers", Setting::TypeInt );
    00623                                 kilometers_setting = kilometers;
    00624
    00625                                 UI::save( passengers, set, "passengers" );
    00626                         }*/
    00627
    00628         };
    00629
    00630         UIREGISTER(mgnorm<chmat>);
    00631
    00632
    00640         class mlstudent : public mlnorm<ldmat>
    00641         {
    00642                 protected:
    00643                         ldmat Lambda;
    00644                         ldmat &_R;
    00645                         ldmat Re;
    00646                 public:
    00647                         mlstudent ( ) :mlnorm<ldmat> (),
    00648                                         Lambda (),      _R ( epdf._R() ) {}
    00649                         void set_parameters ( const mat &A0, const vec &mu0, const ldmat &R0, const ldmat& Lambda0 )
    00650                         {
    00651                                 it_assert_debug ( A0.rows() ==mu0.length(),"" );
    00652                                 it_assert_debug ( R0.rows() ==A0.rows(),"" );
    00653
    00654                                 epdf.set_parameters ( mu0,Lambda ); //
    00655                                 A = A0;
    00656                                 mu_const = mu0;
    00657                                 Re=R0;
    00658                                 Lambda = Lambda0;
    00659                         }
    00660                         void condition ( const vec &cond )
    00661                         {
    00662                                 _mu = A*cond + mu_const;
    00663                                 double zeta;
    00664                                 //ugly hack!
    00665                                 if ( ( cond.length() +1 ) ==Lambda.rows() )
    00666                                 {
    00667                                         zeta = Lambda.invqform ( concat ( cond, vec_1 ( 1.0 ) ) );
    00668                                 }
    00669                                 else
    00670                                 {
    00671                                         zeta = Lambda.invqform ( cond );
    00672                                 }
    00673                                 _R = Re;
    00674                                 _R*= ( 1+zeta );// / ( nu ); << nu is in Re!!!!!!
    00675                         };
    00676
    00677         };
    00687         class mgamma : public mEF
    00688         {
    00689                 protected:
    00691                         egamma epdf;
    00693                         double k;
    00695                         vec &_beta;
    00696
    00697                 public:
    00699                         mgamma ( ) : mEF ( ), epdf (), _beta ( epdf._beta() ) {ep=&epdf;};
    00701                         void set_parameters ( double k, const vec &beta0 );
    00702                         void condition ( const vec &val ) {_beta=k/val;};
    00712                         void from_setting(const Setting &set){
    00713                                 mpdf::from_setting(set); // reads rv and rvc
    00714                                 vec betatmp; // ugly but necessary
    00715                                 UI::get(betatmp,set,"beta");
    00716                                 set.lookupValue("k",k);
    00717                                 set_parameters(k,betatmp);
    00718                         }
    00719         };
    00720         UIREGISTER(mgamma);
    00721
    00731         class migamma : public mEF
    00732         {
    00733                 protected:
    00735                         eigamma epdf;
    00737                         double k;
    00739                         vec &_alpha;
    00741                         vec &_beta;
    00742
    00743                 public:
    00746                         migamma ( ) : mEF (), epdf ( ), _alpha ( epdf._alpha() ), _beta ( epdf._beta() ) {ep=&epdf;};
    00747                         migamma ( const migamma &m ) : mEF (), epdf ( m.epdf ), _alpha ( epdf._alpha() ), _beta ( epdf._beta() ) {ep=&epdf;};
    00749
    00751                         void set_parameters ( int len, double k0 )
    00752                         {
    00753                                 k=k0;
    00754                                 epdf.set_parameters ( ( 1.0/ ( k*k ) +2.0 ) *ones ( len ) /*alpha*/, ones ( len ) /*beta*/ );
    00755                                 dimc = dimension();
    00756                         };
    00757                         void condition ( const vec &val )
    00758                         {
    00759                                 _beta=elem_mult ( val, ( _alpha-1.0 ) );
    00760                         };
    00761         };
    00762
    00763
    00775         class mgamma_fix : public mgamma
    00776         {
    00777                 protected:
    00779                         double l;
    00781                         vec refl;
    00782                 public:
    00784                         mgamma_fix ( ) : mgamma ( ),refl () {};
    00786                         void set_parameters ( double k0 , vec ref0, double l0 )
    00787                         {
    00788                                 mgamma::set_parameters ( k0, ref0 );
    00789                                 refl=pow ( ref0,1.0-l0 );l=l0;
    00790                                 dimc=dimension();
    00791                         };
    00792
    00793                         void condition ( const vec &val ) {vec mean=elem_mult ( refl,pow ( val,l ) ); _beta=k/mean;};
    00794         };
    00795
    00796
    00809         class migamma_ref : public migamma
    00810         {
    00811                 protected:
    00813                         double l;
    00815                         vec refl;
    00816                 public:
    00818                         migamma_ref ( ) : migamma (),refl ( ) {};
    00820                         void set_parameters ( double k0 , vec ref0, double l0 )
    00821                         {
    00822                                 migamma::set_parameters ( ref0.length(), k0 );
    00823                                 refl=pow ( ref0,1.0-l0 );
    00824                                 l=l0;
    00825                                 dimc = dimension();
    00826                         };
    00827
    00828                         void condition ( const vec &val )
    00829                         {
    00830                                 vec mean=elem_mult ( refl,pow ( val,l ) );
    00831                                 migamma::condition ( mean );
    00832                         };
    00833
    00854                         void from_setting( const Setting &set );
    00855
    00856                         // TODO dodelat void to_setting( Setting &set ) const;
    00857         };
    00858
    00859
    00860         UIREGISTER(migamma_ref);
    00861
    00871         class elognorm: public enorm<ldmat>
    00872         {
    00873                 public:
    00874                         vec sample() const {return exp ( enorm<ldmat>::sample() );};
    00875                         vec mean() const {vec var=enorm<ldmat>::variance();return exp ( mu - 0.5*var );};
    00876
    00877         };
    00878
    00890         class mlognorm : public mpdf
    00891         {
    00892                 protected:
    00893                         elognorm eno;
    00895                         double sig2;
    00897                         vec &mu;
    00898                 public:
    00900                         mlognorm ( ) : eno (), mu ( eno._mu() ) {ep=&eno;};
    00902                         void set_parameters ( int size, double k )
    00903                         {
    00904                                 sig2 = 0.5*log ( k*k+1 );
    00905                                 eno.set_parameters ( zeros ( size ),2*sig2*eye ( size ) );
    00906
    00907                                 dimc = size;
    00908                         };
    00909
    00910                         void condition ( const vec &val )
    00911                         {
    00912                                 mu=log ( val )-sig2;//elem_mult ( refl,pow ( val,l ) );
    00913                         };
    00914
    00933                         void from_setting( const Setting &set );
    00934
    00935                         // TODO dodelat void to_setting( Setting &set ) const;
    00936
    00937         };
    00938
    00939         UIREGISTER(mlognorm);
    00940
    00944         class eWishartCh : public epdf
    00945         {
    00946                 protected:
    00948                         chmat Y;
    00950                         int p;
    00952                         double delta;
    00953                 public:
    00954                         void set_parameters ( const mat &Y0, const double delta0 ) {Y=chmat ( Y0 );delta=delta0; p=Y.rows(); dim = p*p; }
    00955                         mat sample_mat() const
    00956                         {
    00957                                 mat X=zeros ( p,p );
    00958
    00959                                 //sample diagonal
    00960                                 for ( int i=0;i<p;i++ )
    00961                                 {
    00962                                         GamRNG.setup ( 0.5* ( delta-i ) , 0.5 ); // no +1 !! index if from 0
    00963 #pragma omp critical
    00964                                         X ( i,i ) =sqrt ( GamRNG() );
    00965                                 }
    00966                                 //do the rest
    00967                                 for ( int i=0;i<p;i++ )
    00968                                 {
    00969                                         for ( int j=i+1;j<p;j++ )
    00970                                         {
    00971 #pragma omp critical
    00972                                                 X ( i,j ) =NorRNG.sample();
    00973                                         }
    00974                                 }
    00975                                 return X*Y._Ch();// return upper triangular part of the decomposition
    00976                         }
    00977                         vec sample () const
    00978                         {
    00979                                 return vec ( sample_mat()._data(),p*p );
    00980                         }
    00982                         void setY ( const mat &Ch0 ) {copy_vector ( dim,Ch0._data(), Y._Ch()._data() );}
    00984                         void _setY ( const vec &ch0 ) {copy_vector ( dim, ch0._data(), Y._Ch()._data() ); }
    00986                         const chmat& getY()const {return Y;}
    00987         };
    00988
    00989         class eiWishartCh: public epdf
    00990         {
    00991                 protected:
    00992                         eWishartCh W;
    00993                         int p;
    00994                         double delta;
    00995                 public:
    00996                         void set_parameters ( const mat &Y0, const double delta0) {
    00997                                 delta = delta0;
    00998                                 W.set_parameters ( inv ( Y0 ),delta0 );
    00999                                 dim = W.dimension(); p=Y0.rows();
    01000                         }
    01001                         vec sample() const {mat iCh; iCh=inv ( W.sample_mat() ); return vec ( iCh._data(),dim );}
    01002                         void _setY ( const vec &y0 )
    01003                         {
    01004                                 mat Ch ( p,p );
    01005                                 mat iCh ( p,p );
    01006                                 copy_vector ( dim, y0._data(), Ch._data() );
    01007
    01008                                 iCh=inv ( Ch );
    01009                                 W.setY ( iCh );
    01010                         }
    01011                         virtual double evallog ( const vec &val ) const {
    01012                                 chmat X(p);
    01013                                 const chmat& Y=W.getY();
    01014
    01015                                 copy_vector(p*p,val._data(),X._Ch()._data());
    01016                                 chmat iX(p);X.inv(iX);
    01017                                 // compute
    01018 //                              \frac{ |\Psi|^{m/2}|X|^{-(m+p+1)/2}e^{-tr(\Psi X^{-1})/2} }{ 2^{mp/2}\Gamma_p(m/2)},
    01019                                 mat M=Y.to_mat()*iX.to_mat();
    01020
    01021                                 double log1 = 0.5*p*(2*Y.logdet())-0.5*(delta+p+1)*(2*X.logdet())-0.5*trace(M);
    01022                                 //Fixme! Multivariate gamma omitted!! it is ok for sampling, but not otherwise!!
    01023
    01024 /*                              if (0) {
    01025                                         mat XX=X.to_mat();
    01026                                         mat YY=Y.to_mat();
    01027
    01028                                         double log2 = 0.5*p*log(det(YY))-0.5*(delta+p+1)*log(det(XX))-0.5*trace(YY*inv(XX));
    01029                                         cout << log1 << "," << log2 << endl;
    01030                                 }*/
    01031                                 return log1;
    01032                         };
    01033
    01034         };
    01035
    01036         class rwiWishartCh : public mpdf
    01037         {
    01038                 protected:
    01039                         eiWishartCh eiW;
    01041                         double sqd;
    01042                         //reference point for diagonal
    01043                         vec refl;
    01044                         double l;
    01045                         int p;
    01046                 public:
    01047                         void set_parameters ( int p0, double k, vec ref0, double l0  )
    01048                         {
    01049                                 p=p0;
    01050                                 double delta = 2/(k*k)+p+3;
    01051                                 sqd=sqrt ( delta-p-1 );
    01052                                 l=l0;
    01053                                 refl=pow(ref0,1-l);
    01054
    01055                                 eiW.set_parameters ( eye ( p ),delta );
    01056                                 ep=&eiW;
    01057                                 dimc=eiW.dimension();
    01058                         }
    01059                         void condition ( const vec &c ) {
    01060                                 vec z=c;
    01061                                 int ri=0;
    01062                                 for(int i=0;i<p*p;i+=(p+1)){//trace diagonal element
    01063                                         z(i) = pow(z(i),l)*refl(ri);
    01064                                         ri++;
    01065                                 }
    01066
    01067                                 eiW._setY ( sqd*z );
    01068                         }
    01069         };
    01070
    01072         enum RESAMPLING_METHOD { MULTINOMIAL = 0, STRATIFIED = 1, SYSTEMATIC = 3 };
    01078         class eEmp: public epdf
    01079         {
    01080                 protected :
    01082                         int n;
    01084                         vec w;
    01086                         Array<vec> samples;
    01087                 public:
    01090                         eEmp ( ) :epdf ( ),w ( ),samples ( ) {};
    01092                         eEmp ( const eEmp &e ) : epdf ( e ), w ( e.w ), samples ( e.samples ) {};
    01094
    01096                         void set_statistics ( const vec &w0, const epdf* pdf0 );
    01098                         void set_statistics ( const epdf* pdf0 , int n ) {set_statistics ( ones ( n ) /n,pdf0 );};
    01100                         void set_samples ( const epdf* pdf0 );
    01102                         void set_parameters ( int n0, bool copy=true ) {n=n0; w.set_size ( n0,copy );samples.set_size ( n0,copy );};
    01104                         vec& _w()  {return w;};
    01106                         const vec& _w() const {return w;};
    01108                         Array<vec>& _samples() {return samples;};
    01110                         const Array<vec>& _samples() const {return samples;};
    01112                         ivec resample ( RESAMPLING_METHOD method=SYSTEMATIC );
    01114                         vec sample() const {it_error ( "Not implemented" );return 0;}
    01116                         double evallog ( const vec &val ) const {it_error ( "Not implemented" );return 0.0;}
    01117                         vec mean() const
    01118                         {
    01119                                 vec pom=zeros ( dim );
    01120                                 for ( int i=0;i<n;i++ ) {pom+=samples ( i ) *w ( i );}
    01121                                 return pom;
    01122                         }
    01123                         vec variance() const
    01124                         {
    01125                                 vec pom=zeros ( dim );
    01126                                 for ( int i=0;i<n;i++ ) {pom+=pow ( samples ( i ),2 ) *w ( i );}
    01127                                 return pom-pow ( mean(),2 );
    01128                         }
    01130                         void qbounds ( vec &lb, vec &ub, double perc=0.95 ) const
    01131                         {
    01132                                 // lb in inf so than it will be pushed below;
    01133                                 lb.set_size ( dim );
    01134                                 ub.set_size ( dim );
    01135                                 lb = std::numeric_limits<double>::infinity();
    01136                                 ub = -std::numeric_limits<double>::infinity();
    01137                                 int j;
    01138                                 for ( int i=0;i<n;i++ )
    01139                                 {
    01140                                         for ( j=0;j<dim; j++ )
    01141                                         {
    01142                                                 if ( samples ( i ) ( j ) <lb ( j ) ) {lb ( j ) =samples ( i ) ( j );}
    01143                                                 if ( samples ( i ) ( j ) >ub ( j ) ) {ub ( j ) =samples ( i ) ( j );}
    01144                                         }
    01145                                 }
    01146                         }
    01147         };
    01148
    01149
    01151
    01152         template<class sq_T>
    01153         void enorm<sq_T>::set_parameters ( const vec &mu0, const sq_T &R0 )
    01154         {
    01155 //Fixme test dimensions of mu0 and R0;
    01156                 mu = mu0;
    01157                 R = R0;
    01158                 validate();
    01159         };
    01160
    01161         template<class sq_T>
    01162         void enorm<sq_T>::from_setting(const Setting &set){
    01163                 epdf::from_setting(set); //reads rv
    01164
    01165                 UI::get(mu,set,"mu");
    01166                 mat Rtmp;// necessary for conversion
    01167                 UI::get(Rtmp,set,"R");
    01168                 R=Rtmp; // conversion
    01169                 validate();
    01170         }
    01171
    01172         template<class sq_T>
    01173         void enorm<sq_T>::dupdate ( mat &v, double nu )
    01174         {
    01175                 //
    01176         };
    01177
    01178 // template<class sq_T>
    01179 // void enorm<sq_T>::tupdate ( double phi, mat &vbar, double nubar ) {
    01180 //      //
    01181 // };
    01182
    01183         template<class sq_T>
    01184         vec enorm<sq_T>::sample() const
    01185         {
    01186                 vec x ( dim );
    01187 #pragma omp critical
    01188                 NorRNG.sample_vector ( dim,x );
    01189                 vec smp = R.sqrt_mult ( x );
    01190
    01191                 smp += mu;
    01192                 return smp;
    01193         };
    01194
    01195         template<class sq_T>
    01196         mat enorm<sq_T>::sample ( int N ) const
    01197         {
    01198                 mat X ( dim,N );
    01199                 vec x ( dim );
    01200                 vec pom;
    01201                 int i;
    01202
    01203                 for ( i=0;i<N;i++ )
    01204                 {
    01205 #pragma omp critical
    01206                         NorRNG.sample_vector ( dim,x );
    01207                         pom = R.sqrt_mult ( x );
    01208                         pom +=mu;
    01209                         X.set_col ( i, pom );
    01210                 }
    01211
    01212                 return X;
    01213         };
    01214
    01215 // template<class sq_T>
    01216 // double enorm<sq_T>::eval ( const vec &val ) const {
    01217 //      double pdfl,e;
    01218 //      pdfl = evallog ( val );
    01219 //      e = exp ( pdfl );
    01220 //      return e;
    01221 // };
    01222
    01223         template<class sq_T>
    01224         double enorm<sq_T>::evallog_nn ( const vec &val ) const
    01225         {
    01226                 // 1.83787706640935 = log(2pi)
    01227                 double tmp=-0.5* ( R.invqform ( mu-val ) );// - lognc();
    01228                 return  tmp;
    01229         };
    01230
    01231         template<class sq_T>
    01232         inline double enorm<sq_T>::lognc () const
    01233         {
    01234                 // 1.83787706640935 = log(2pi)
    01235                 double tmp=0.5* ( R.cols() * 1.83787706640935 +R.logdet() );
    01236                 return tmp;
    01237         };
    01238
    01239         template<class sq_T>
    01240         void mlnorm<sq_T>::set_parameters ( const mat &A0, const vec &mu0, const sq_T &R0 )
    01241         {
    01242                 it_assert_debug ( A0.rows() ==mu0.length(),"" );
    01243                 it_assert_debug ( A0.rows() ==R0.rows(),"" );
    01244
    01245                 epdf.set_parameters ( zeros ( A0.rows() ),R0 );
    01246                 A = A0;
    01247                 mu_const = mu0;
    01248                 dimc=A0.cols();
    01249         }
    01250
    01251 // template<class sq_T>
    01252 // vec mlnorm<sq_T>::samplecond (const  vec &cond, double &lik ) {
    01253 //      this->condition ( cond );
    01254 //      vec smp = epdf.sample();
    01255 //      lik = epdf.eval ( smp );
    01256 //      return smp;
    01257 // }
    01258
    01259 // template<class sq_T>
    01260 // mat mlnorm<sq_T>::samplecond (const vec &cond, vec &lik, int n ) {
    01261 //      int i;
    01262 //      int dim = rv.count();
    01263 //      mat Smp ( dim,n );
    01264 //      vec smp ( dim );
    01265 //      this->condition ( cond );
    01266 //
    01267 //      for ( i=0; i<n; i++ ) {
    01268 //              smp = epdf.sample();
    01269 //              lik ( i ) = epdf.eval ( smp );
    01270 //              Smp.set_col ( i ,smp );
    01271 //      }
    01272 //
    01273 //      return Smp;
    01274 // }
    01275
    01276         template<class sq_T>
    01277         void mlnorm<sq_T>::condition ( const vec &cond )
    01278         {
    01279                 _mu = A*cond + mu_const;
    01280 //R is already assigned;
    01281         }
    01282
    01283         template<class sq_T>
    01284         enorm<sq_T>* enorm<sq_T>::marginal ( const RV &rvn ) const
    01285         {
    01286                 it_assert_debug ( isnamed(), "rv description is not assigned" );
    01287                 ivec irvn = rvn.dataind ( rv );
    01288
    01289                 sq_T Rn ( R,irvn ); //select rows and columns of R
    01290
    01291                 enorm<sq_T>* tmp = new enorm<sq_T>;
    01292                 tmp->set_rv ( rvn );
    01293                 tmp->set_parameters ( mu ( irvn ), Rn );
    01294                 return tmp;
    01295         }
    01296
    01297         template<class sq_T>
    01298         mpdf* enorm<sq_T>::condition ( const RV &rvn ) const
    01299         {
    01300
    01301                 it_assert_debug ( isnamed(),"rvs are not assigned" );
    01302
    01303                 RV rvc = rv.subt ( rvn );
    01304                 it_assert_debug ( ( rvc._dsize() +rvn._dsize() ==rv._dsize() ),"wrong rvn" );
    01305                 //Permutation vector of the new R
    01306                 ivec irvn = rvn.dataind ( rv );
    01307                 ivec irvc = rvc.dataind ( rv );
    01308                 ivec perm=concat ( irvn , irvc );
    01309                 sq_T Rn ( R,perm );
    01310
    01311                 //fixme - could this be done in general for all sq_T?
    01312                 mat S=Rn.to_mat();
    01313                 //fixme
    01314                 int n=rvn._dsize()-1;
    01315                 int end=R.rows()-1;
    01316                 mat S11 = S.get ( 0,n, 0, n );
    01317                 mat S12 = S.get ( 0, n , rvn._dsize(), end );
    01318                 mat S22 = S.get ( rvn._dsize(), end, rvn._dsize(), end );
    01319
    01320                 vec mu1 = mu ( irvn );
    01321                 vec mu2 = mu ( irvc );
    01322                 mat A=S12*inv ( S22 );
    01323                 sq_T R_n ( S11 - A *S12.T() );
    01324
    01325                 mlnorm<sq_T>* tmp=new mlnorm<sq_T> ( );
    01326                 tmp->set_rv ( rvn ); tmp->set_rvc ( rvc );
    01327                 tmp->set_parameters ( A,mu1-A*mu2,R_n );
    01328                 return tmp;
    01329         }
    01330
    01332
    01333         template<class sq_T>
    01334         std::ostream &operator<< ( std::ostream &os,  mlnorm<sq_T> &ml )
    01335         {
    01336                 os << "A:"<< ml.A<<endl;
    01337                 os << "mu:"<< ml.mu_const<<endl;
    01338                 os << "R:" << ml.epdf._R().to_mat() <<endl;
    01339                 return os;
    01340         };
    01341
    01342 }
    01343 #endif //EF_H