Changeset 170

Show
Ignore:
Timestamp:
09/24/08 13:07:50 (16 years ago)
Author:
smidl
Message:

Mixtures of EF and related changes to libEF and BM

Files:
23 modified

Legend:

Unmodified
Added
Removed

  • CMakeLists.txt

    r169 r170  
    2727        SET(CMAKE_BUILD_TYPE Debug) 
    2828        SET(CMAKE_CXX_FLAGS_DEBUG "-g -Wall -Wno-unknown-pragmas") 
    29         SET(ITPP_DIR "/home/smidl/work/git/trunk/itpp/.libs") 
    3029ENDIF(WIN32) 
    3130 

  • bdm/estim/arx.cpp

    r162 r170  
    33using namespace itpp; 
    44 
    5 void ARX::bayes ( const vec &dt ) { 
     5void ARX::bayes ( const vec &dt, const double w ) { 
    66        double lnc; 
    77 
    8         if (frg<1.0) { 
    9                 V*=frg; 
    10                 nu*=frg; 
    11                 if (evalll) { 
     8        if ( frg<1.0 ) { 
     9                est.pow ( frg ); 
     10                if ( evalll ) { 
    1211                        last_lognc = est.lognc(); 
    1312                } 
    1413        } 
    15         V.opupdt ( dt,1.0); 
    16         nu+=1.0; 
     14        V.opupdt ( dt,w ); 
     15        nu+=w; 
    1716 
    1817        if ( evalll ) { 
    … …  
    2019                ll = lnc - last_lognc; 
    2120                last_lognc = lnc; 
    22                 tll +=ll; 
    2321        } 
    2422} 
    2523 
     24double ARX::logpred ( const vec &dt ) const { 
     25        egiw pred ( est ); 
     26        ldmat &V=pred._V(); 
     27        double &nu=pred._nu(); 
     28 
     29        double lll; 
     30 
     31        if ( frg<1.0 ) { 
     32                pred.pow ( frg ); 
     33                lll = pred.lognc(); 
     34        } 
     35        else  
     36                if(evalll){lll=last_lognc;}else{lll=pred.lognc();} 
     37 
     38        V.opupdt ( dt,1.0 ); 
     39        nu+=1.0; 
     40 
     41        return pred.lognc()-lll; 
     42} 
     43 
     44BM* ARX::_copy_ ( bool changerv) { 
     45        ARX* Tmp=new ARX ( *this ); 
     46        if ( changerv ) {Tmp->rv.newids(); Tmp->est._renewrv(Tmp->rv);} 
     47        return Tmp; 
     48} 
     49 
     50void ARX::set_statistics ( const BMEF* B0 ) { 
     51        const ARX* A0=dynamic_cast<const ARX*> ( B0 ); 
     52 
     53        it_assert_debug ( V.rows() ==A0->V.rows(),"ARX::set_statistics Statistics  differ" ); 
     54        set_parameters ( A0->V,A0->nu ); 
     55} 
    2656/*! \brief Return the best structure 
    2757@param Eg a copy of GiW density that is being examined 
    … …  
    4878 
    4979 
    50 cout << "bb:(" << indeces <<") ll=" << Egll <<endl; 
     80        cout << "bb:(" << indeces <<") ll=" << Egll <<endl; 
    5181 
    5282        //try to remove only one rv 
    … …  
    6292                tmpll = Eg.lognc()-Eg0.lognc(); // likelihood is difference of norm. coefs. 
    6393 
    64 cout << "i=(" << i <<") ll=" << tmpll <<endl; 
    65                  
     94                cout << "i=(" << i <<") ll=" << tmpll <<endl; 
     95 
    6696                // 
    6797                if ( tmpll > Egll ) { //increase of the likelihood 

  • bdm/estim/arx.h

    r168 r170  
    3434Extension for time-variant parameters \f$\theta_t,r_t\f$ may be achived using exponential forgetting (Kulhavy and Zarrop, 1993). In such a case, the forgetting factor \c frg \f$\in <0,1>\f$ should be given in the constructor. Time-invariant parameters are estimated for \c frg = 1. 
    3535*/ 
    36 class ARX: public BM { 
     36class ARX: public BMEF { 
    3737protected: 
    3838        //! Posterior estimate of \f$\theta,r\f$ in the form of Normal-inverse Wishart density 
    … …  
    4242        //! cached value of est.nu 
    4343        double &nu; 
    44         //! forgetting factor  
    45         double frg; 
    46         //! cached value of lognc() in the previous step (used in evaluation of \c ll ) 
    47         double last_lognc; 
    48         //! total likelihood 
    49         double tll; 
    5044public: 
    5145        //! Full constructor 
    52         ARX (const RV &rv, const mat &V0, const double &nu0, const double frg0=1.0) : BM(rv),est(rv,V0,nu0), V(est._V()), nu(est._nu()), frg(frg0){last_lognc=est.lognc();tll=0.0;}; 
     46        ARX ( const RV &rv, const mat &V0, const double &nu0, const double frg0=1.0 ) : BMEF ( rv,frg0 ),est ( rv,V0,nu0 ), V ( est._V() ), nu ( est._nu() ) 
     47        {last_lognc=est.lognc();}; 
     48 
     49        //!Copy constructor 
     50        ARX ( const ARX &A0 ) : BMEF ( A0),est ( rv,A0.V,A0.nu ), V ( est._V() ), nu ( est._nu() ) {}; 
     51 
     52        //!Auxiliary function 
     53        BM* _copy_(bool changerv=false); 
     54         
     55//      //! Set parameters given by moments, \c mu (mean of theta), \c R (mean of R) and \c C (variance of theta) 
     56//      void set_parameters ( const vec &mu, const mat &R, const mat &C, double dfm){}; 
    5357        //! Set sufficient statistics 
    54         void set_parameters(const mat &V0, const double &nu0){est._V()=V0;est._nu()=nu0;last_lognc=est.lognc();tll=last_lognc;} 
     58        void set_parameters ( const ldmat &V0, const double &nu0 ) 
     59        {est._V() =V0;est._nu() =nu0;last_lognc=est.lognc();} 
     60        void set_statistics ( const BMEF* BM0 ); 
    5561        //! Returns sufficient statistics 
    56         void get_parameters(mat &V0, double &nu0){V0=est._V().to_mat(); nu0=est._nu();} 
     62        void get_parameters ( mat &V0, double &nu0 ) {V0=est._V().to_mat(); nu0=est._nu();} 
    5763        //! Here \f$dt = [y_t psi_t] \f$. 
    58         void bayes ( const vec &dt ); 
    59         epdf& _epdf() {return est;} 
     64        void bayes ( const vec &dt, const double w ); 
     65        void bayes ( const vec &dt ) {bayes ( dt,1.0 );}; 
     66        const epdf& _epdf() const {return est;} 
     67        double logpred ( const vec &dt ) const; 
     68        void flatten (BMEF* B ) { 
     69                ARX* A=dynamic_cast<ARX*>(B); 
     70                // nu should be equal to B.nu 
     71                est.pow ( A->nu/nu); 
     72                if(evalll){last_lognc=est.lognc();} 
     73        } 
     74 
    6075        //! Brute force structure estimation.\return indeces of accepted regressors. 
    61         ivec structure_est(egiw Eg0); 
    62         //!access function 
    63         double _tll(){return tll;} 
     76        ivec structure_est ( egiw Eg0 ); 
    6477}; 
    6578 

  • bdm/estim/libKF.h

    r132 r170  
    121121        void bayes ( const vec &dt ); 
    122122        //!access function 
    123         epdf& _epdf() {return est;} 
     123        const epdf& _epdf() const {return est;} 
    124124        //!access function 
    125125        mat& __K() {return _K;} 
    … …  
    186186        void set_est (vec mu0, mat P0){mu=mu0;P=P0;}; 
    187187        //!dummy! 
    188         epdf& _epdf(){return E;}; 
     188        const epdf& _epdf()const{return E;}; 
    189189}; 
    190190 

  • bdm/estim/libPF.h

    r141 r170  
    7171                eEmp &E; 
    7272                vec &_w; 
    73                 Array<epdf*> Coms; 
     73                Array<const epdf*> Coms; 
    7474        public: 
    7575                mpfepdf ( eEmp &E0, const RV &rvc ) : 
    … …  
    8080                }; 
    8181 
    82                 void set_elements ( int &i, double wi, epdf* ep ) 
     82                void set_elements ( int &i, double wi, const epdf* ep ) 
    8383                {_w ( i ) =wi; Coms ( i ) =ep;}; 
    8484 
    … …  
    112112                for ( int i=0;i<n;i++ ) { 
    113113                        Bms[i] = new BM_T ( BMcond0 ); //copy constructor 
    114                         epdf& pom=Bms[i]->_epdf(); 
     114                        const epdf& pom=Bms[i]->_epdf(); 
    115115                        jest.set_elements ( i,1.0/n,&pom ); 
    116116                } 
    … …  
    121121 
    122122        void bayes ( const vec &dt ); 
    123         epdf& _epdf() {return jest;} 
     123        const epdf& _epdf() const {return jest;} 
    124124        //! Set postrior of \c rvc to samples from epdf0. Statistics of Bms are not re-computed! Use only for initialization! 
    125125        void set_est ( const epdf& epdf0 ) { 
    … …  
    197197                                delete Bms[i]; 
    198198                                Bms[i] = new BM_T ( *Bms[ind ( i ) ] ); //copy constructor 
    199                                 epdf& pom=Bms[i]->_epdf(); 
     199                                const epdf& pom=Bms[i]->_epdf(); 
    200200                                jest.set_elements ( i,1.0/n,&pom ); 
    201201                        } 

  • bdm/math/libDC.cpp

    r168 r170  
    157157        double x = 0.0, sum; 
    158158        int i,j; 
    159         vec s(v.length()); 
    160         vec S=L*v; 
    161159 
    162160        for ( i=0; i<D.length(); i++ ) { //rows of L 
    … …  
    164162                for ( j=0; j<=i; j++ ){sum+=L( i,j )*v( j );} 
    165163                x +=D( i )*sum*sum; 
    166                 s(i)=sum; 
    167164        }; 
    168165        return x; 

  • bdm/stat/emix.h

    r168 r170  
    139139protected: 
    140140        //! Components (epdfs) 
    141         Array<epdf*> epdfs; 
     141        Array<const epdf*> epdfs; 
    142142        //! Array of indeces 
    143143        Array<ivec> rvinds; 
    144144public: 
    145         eprod ( const Array<epdf*> epdfs0 ) : epdf ( RV() ),epdfs ( epdfs0 ),rvinds ( epdfs.length() ) { 
     145        eprod ( const Array<const epdf*> epdfs0 ) : epdf ( RV() ),epdfs ( epdfs0 ),rvinds ( epdfs.length() ) { 
    146146                bool independent=true; 
    147147                for ( int i=0;i<epdfs.length();i++ ) { 
    … …  
    175175                return tmp; 
    176176        } 
     177        //!access function 
     178        const epdf* operator () (int i) const {it_assert_debug(i<epdfs.length(),"wrong index");return epdfs(i);} 
    177179}; 
    178180 

  • bdm/stat/libBM.cpp

    r168 r170  
    180180        return pom; 
    181181} 
     182 
     183void RV::newids(){ids=linspace ( RVcounter+1, RVcounter+len ),RVcounter+=len;} 
     184 
     185void BM::bayesB(const mat &Data){ 
     186        for(int t=0;t<Data.cols();t++){bayes(Data.get_col(t));} 
     187} 

  • bdm/stat/libBM.h

    r168 r170  
    103103        //!access function 
    104104        std::string name ( int at ) {return names ( at );}; 
     105        //!Assign unused ids to this rv  
     106        void newids(); 
    105107}; 
    106108 
    … …  
    144146        epdf ( const RV &rv0 ) :rv ( rv0 ) {}; 
    145147 
    146         //! Returns the required moment of the epdf 
     148//      //! Returns the required moment of the epdf 
    147149//      virtual vec moment ( const int order = 1 ); 
     150         
    148151        //! Returns a sample, \f$x\f$ from density \f$epdf(rv)\f$ 
    149152        virtual vec sample () const =0; 
    … …  
    156159        virtual double evalpdflog ( const vec &val ) const =0; 
    157160 
     161        //! Compute log-probability of multiple values argument \c val 
     162        virtual vec evalpdflog ( const mat &Val ) const { 
     163                vec x ( Val.cols() ); 
     164                for ( int i=0;i<Val.cols();i++ ) {x ( i ) =evalpdflog( Val.get_col(i) ) ;} 
     165                return x; 
     166        } 
     167 
    158168        //! return expected value 
    159169        virtual vec mean() const =0; 
    … …  
    162172        virtual ~epdf() {}; 
    163173        //! access function, possibly dangerous! 
    164         RV& _rv() {return rv;} 
     174        const RV& _rv() const {return rv;} 
     175        //! modifier function - useful when copying epdfs 
     176        void _renewrv(const RV &in_rv){rv=in_rv;} 
    165177}; 
    166178 
    … …  
    270282 
    271283        //!Default constructor 
    272         BM ( const RV &rv0 ) :rv ( rv0 ), ll ( 0 ),evalll ( true ) {//Fixme: test rv 
     284        BM ( const RV &rv0, double ll0=0,bool evalll0=true ) :rv ( rv0 ), ll ( ll0 ),evalll ( evalll0) {//Fixme: test rv 
    273285        }; 
     286        //!Copy constructor 
     287        BM (const BM &B) : rv(B.rv), ll(B.ll), evalll(B.evalll) {} 
    274288 
    275289        /*! \brief Incremental Bayes rule 
    … …  
    278292        virtual void bayes ( const vec &dt ) = 0; 
    279293        //! Batch Bayes rule (columns of Dt are observations) 
    280         void bayes ( mat Dt ); 
     294        virtual void bayesB (const mat &Dt ); 
    281295        //! Returns a pointer to the epdf representing posterior density on parameters. Use with care! 
    282         virtual epdf& _epdf() =0; 
    283  
     296        virtual const epdf& _epdf() const =0; 
     297 
     298        //! Evaluates predictive log-likelihood of the given data record 
     299        //! I.e. marginal likelihood of the data with the posterior integrated out. 
     300        virtual double logpred(const vec &dt)const{it_error("Not implemented");return 0.0;} 
     301         
    284302        //! Destructor for future use; 
    285303        virtual ~BM() {}; 
    … …  
    290308        //!access function  
    291309        void set_evalll(bool evl0){evalll=evl0;} 
     310         
     311        //! Copy function required in vectors, Arrays of BM etc. Have to be DELETED manually! 
     312        //! Prototype: BM* _copy_(){BM Tmp*=new Tmp(this*);  return Tmp; } 
     313        virtual BM* _copy_(bool changerv=false){it_error("function _copy_ not implemented for this BM"); return NULL;}; 
    292314}; 
    293315 

  • bdm/stat/libEF.cpp

    r168 r170  
    1212using std::cout; 
    1313 
     14void BMEF::bayes ( const vec &dt ) {this->bayes ( dt,1.0 );}; 
     15 
    1416vec egiw::sample() const { 
    1517        it_warning ( "Function not implemented" ); 
    … …  
    1719} 
    1820 
    19 double egiw::evalpdflog ( const vec &val ) const { 
     21double egiw::evalpdflog_nn ( const vec &val ) const { 
     22        int vend = val.length()-1; 
    2023 
    2124        if ( xdim==1 ) { //same as the following, just quicker. 
    22                 double r = val ( val.length()-1 ); //last entry! 
     25                double r = val ( vend ); //last entry! 
    2326                vec Psi ( nPsi+xdim ); 
    2427                Psi ( 0 ) = -1.0; 
    25                 Psi.set_subvector ( 1,val ); // fill the rest 
    26  
    27                 return -0.5* ( nu*log ( r ) + V.qform ( Psi ) /r ) + lognc(); 
     28                Psi.set_subvector ( 1,val ( 0,vend-1 ) ); // fill the rest 
     29 
     30                double Vq=V.qform ( Psi ); 
     31                return -0.5* ( nu*log ( r ) + Vq /r ); 
    2832        } 
    2933        else { 
    30                 int vend = val.length()-1; 
    3134                mat Th= reshape ( val ( 0,nPsi*xdim-1 ),nPsi,xdim ); 
    3235                fsqmat R ( reshape ( val ( nPsi*xdim,vend ),xdim,xdim ) ); 
    … …  
    3437                fsqmat iR ( xdim ); 
    3538                R.inv ( iR ); 
    36                  
    37                 return -0.5* ( nu*R.logdet() + trace ( iR.to_mat() *Tmp.T()*V.to_mat() *Tmp ) ) + lognc(); 
     39 
     40                return -0.5* ( nu*R.logdet() + trace ( iR.to_mat() *Tmp.T() *V.to_mat() *Tmp ) ); 
    3841        } 
    3942} 
    … …  
    5558                     - 0.5*xdim* ( m*log2 + 0.5* ( xdim-1 ) *log2pi )  - lg; 
    5659 
    57         return nkG+nkW; 
     60        return -nkG-nkW; 
    5861} 
    5962 
    6063vec egiw::mean() const { 
     64 
     65        if ( xdim==1 ) { 
     66                const mat &L= V._L(); 
     67                const vec &D= V._D(); 
     68                int end = L.rows()-1; 
     69 
     70                vec m ( rv.count() ); 
     71                mat iLsub = ltuinv ( L ( xdim,end,xdim,end ) ); 
     72 
     73                vec L0 = L.get_col ( 0 ); 
     74 
     75                m.set_subvector ( 0,iLsub*L0 ( 1,end ) ); 
     76                m ( end ) = D ( 0 ) / ( nu -nPsi -2*xdim -2 ); 
     77                return m; 
     78        } else { 
     79                mat M; 
     80                mat R; 
     81                mean_mat(M,R); 
     82                return cvectorize (concat_vertical(M,R)); 
     83        } 
     84 
     85} 
     86void egiw::mean_mat(mat &M, mat&R) const { 
    6187        const mat &L= V._L(); 
    6288        const vec &D= V._D(); 
    6389        int end = L.rows()-1; 
    64  
    65         vec m(rv.count()); 
    66         mat iLsub = ltuinv ( L ( xdim,end,xdim,end ) ); 
     90                 
     91        ldmat ldR ( L ( 0,xdim-1,0,xdim-1 ), D ( 0,xdim-1 ) / ( nu -nPsi -2*xdim -2 ) ); //exp val of R 
     92        mat iLsub=ltuinv ( L ( xdim,end,xdim,end ) ); 
    6793         
    68         if ( xdim==1 ) { 
    69                 vec L0 = L.get_col ( 0 ); 
    70  
    71                 m.set_subvector ( 0,iLsub*L0 ( 1,end ) ); 
    72                 m ( end ) = D ( 0 ) / ( nu -nPsi -2*xdim -2); 
    73         } 
    74         else { 
    75                 ldmat ldR(L(0,xdim,0,xdim), D(0,xdim)/( nu -nPsi -2*xdim -2)); //exp val of R 
    76                 mat iLsub=ltuinv ( L ( xdim,end,xdim,end ) ); 
    77                 mat M = iLsub*L(xdim,end,0,xdim); 
    78                  
    79                 m = cvectorize(concat_vertical(M,ldR.to_mat())); 
    80         } 
    81         return m; 
    82 } 
     94        // set mean value 
     95        M= iLsub*L ( xdim,end,0,xdim-1 ); 
     96        R= ldR.to_mat()  ; 
     97} 
     98 
     99void egiw::pow(double p){V*=p;nu*=p;} 
    83100 
    84101vec egamma::sample() const { 

  • bdm/stat/libEF.h

    r168 r170  
    4141        //! default constructor 
    4242        eEF ( const RV &rv ) :epdf ( rv ) {}; 
    43         //! logarithm of the normalizing constant, \f$\mathcal{I}\f$  
    44         virtual double lognc()const =0; 
     43        //! logarithm of the normalizing constant, \f$\mathcal{I}\f$ 
     44        virtual double lognc() const =0; 
    4545        //!TODO decide if it is really needed 
    46         virtual void tupdate ( double phi, mat &vbar, double nubar ) {}; 
    47         //!TODO decide if it is really needed 
    48         virtual void dupdate ( mat &v,double nu=1.0 ) {}; 
     46        virtual void dupdate ( mat &v ) {it_error ( "Not implemneted" );}; 
     47        //!Evaluate normalized log-probability 
     48        virtual double evalpdflog_nn ( const vec &val ) const{it_error ( "Not implemneted" );return 0.0;}; 
     49        //!Evaluate normalized log-probability 
     50        virtual double evalpdflog ( const vec &val ) const {return evalpdflog_nn ( val )-lognc();} 
     51        //!Evaluate normalized log-probability for many samples 
     52        virtual vec evalpdflog ( const mat &Val ) const { 
     53                vec x ( Val.cols() ); 
     54                for ( int i=0;i<Val.cols();i++ ) {x ( i ) =evalpdflog_nn ( Val.get_col ( i ) ) ;} 
     55                return x-lognc(); 
     56        } 
     57        //!Power of the density, used e.g. to flatten the density 
     58        virtual void pow ( double p ) {it_error ( "Not implemented" );}; 
    4959}; 
    5060 
    … …  
    6070        //! Default constructor 
    6171        mEF ( const RV &rv0, const RV &rvc0 ) :mpdf ( rv0,rvc0 ) {}; 
     72}; 
     73 
     74//! Estimator for Exponential family 
     75class BMEF : public BM { 
     76protected: 
     77        //! forgetting factor 
     78        double frg; 
     79        //! cached value of lognc() in the previous step (used in evaluation of \c ll ) 
     80        double last_lognc; 
     81public: 
     82        //! Default constructor 
     83        BMEF ( const RV &rv, double frg0=1.0 ) :BM ( rv ), frg ( frg0 ) {} 
     84        //! Copy constructor 
     85        BMEF ( const BMEF &B ) :BM ( B ), frg ( B.frg ), last_lognc ( B.last_lognc ) {} 
     86        //!get statistics from another model 
     87        virtual void set_statistics ( const BMEF* BM0 ) {it_error ( "Not implemented" );}; 
     88        //! Weighted update of sufficient statistics (Bayes rule) 
     89        virtual void bayes ( const vec &data, const double w ) {}; 
     90        //original Bayes 
     91        void bayes ( const vec &dt ); 
     92        //!Flatten the posterior 
     93        virtual void flatten ( BMEF * B) {it_error ( "Not implemented" );} 
    6294}; 
    6395 
    … …  
    99131        //! returns a pointer to the internal mean value. Use with Care! 
    100132        vec& _mu() {return mu;} 
    101          
     133 
    102134        //! access function 
    103         void set_mu(const vec mu0) { mu=mu0;} 
     135        void set_mu ( const vec mu0 ) { mu=mu0;} 
    104136 
    105137        //! returns pointers to the internal variance and its inverse. Use with Care! 
    … …  
    128160public: 
    129161        //!Default constructor, assuming 
    130         egiw(RV rv, mat V0, double nu0): eEF(rv), V(V0), nu(nu0) { 
    131                 xdim = rv.count()/V.rows(); 
    132                 it_assert_debug(rv.count()==xdim*V.rows(),"Incompatible V0."); 
     162        egiw ( RV rv, mat V0, double nu0 ) : eEF ( rv ), V ( V0 ), nu ( nu0 ) { 
     163                xdim = rv.count() /V.rows(); 
     164                it_assert_debug ( rv.count() ==xdim*V.rows(),"Incompatible V0." ); 
     165                nPsi = V.rows()-xdim; 
     166        } 
     167        //!Full constructor for V in ldmat form 
     168        egiw ( RV rv, ldmat V0, double nu0 ) : eEF ( rv ), V ( V0 ), nu ( nu0 ) { 
     169                xdim = rv.count() /V.rows(); 
     170                it_assert_debug ( rv.count() ==xdim*V.rows(),"Incompatible V0." ); 
    133171                nPsi = V.rows()-xdim; 
    134172        } 
    … …  
    136174        vec sample() const; 
    137175        vec mean() const; 
     176        void mean_mat ( mat &M, mat&R ) const; 
    138177        //! In this instance, val= [theta, r]. For multivariate instances, it is stored columnwise val = [theta_1 theta_2 ... r_1 r_2 ] 
    139         double evalpdflog ( const vec &val ) const; 
     178        double evalpdflog_nn ( const vec &val ) const; 
    140179        double lognc () const; 
    141180 
    … …  
    145184        //! returns a pointer to the internal statistics. Use with Care! 
    146185        double& _nu() {return nu;} 
    147  
     186        void pow ( double p ); 
     187}; 
     188 
     189/*! \brief Dirichlet posterior density 
     190Continuous Dirichlet density of \f$n\f$-dimensional variable \f$x\f$ 
     191\f[ 
     192f(x|\beta) = \frac{\Gamma[\gamma]}{\prod_{i=1}^{n}\Gamma(\beta_i)} \prod_{i=1}^{n]x_i^(\beta_i-1) 
     193\f] 
     194where \f$\gamma=\sum_i beta_i\f$. 
     195*/ 
     196class eDirich: public eEF { 
     197protected: 
     198        //!sufficient statistics 
     199        vec beta; 
     200public: 
     201        //!Default constructor 
     202        eDirich ( const RV &rv, const vec &beta0 ) : eEF ( rv ),beta ( beta0 ) {it_assert_debug ( rv.count() ==beta.length(),"Incompatible statistics" ); }; 
     203        //! Copy constructor 
     204        eDirich ( const eDirich &D0 ) : eEF ( D0.rv ),beta ( D0.beta ) {}; 
     205        vec sample() const {it_error ( "Not implemented" );return vec_1 ( 0.0 );}; 
     206        vec mean() const {return beta/sum ( beta );}; 
     207        //! In this instance, val= [theta, r]. For multivariate instances, it is stored columnwise val = [theta_1 theta_2 ... r_1 r_2 ] 
     208        double evalpdflog_nn ( const vec &val ) const {return ( beta-1 ) *log ( val );}; 
     209        double lognc () const { 
     210                double gam=sum ( beta ); 
     211                double lgb=0.0; 
     212                for ( int i=0;i<beta.length();i++ ) {lgb+=lgamma ( beta ( i ) );} 
     213                return lgb-lgamma ( gam ); 
     214        }; 
     215        //!access function 
     216        vec& _beta() {return beta;} 
     217}; 
     218 
     219//! Estimator for Multinomial density 
     220class multiBM : public BMEF { 
     221protected: 
     222        //! Conjugate prior and posterior 
     223        eDirich est; 
     224        vec &beta; 
     225public: 
     226        //!Default constructor 
     227        multiBM ( const RV &rv, const vec beta0 ) : BMEF ( rv ),est ( rv,beta0 ),beta ( est._beta() ) {last_lognc=est.lognc();} 
     228        //!Copy constructor 
     229        multiBM ( const multiBM &B ) : BMEF ( B ),est ( rv,B.beta ),beta ( est._beta() ) {} 
     230 
     231        void set_statistics ( const BM* mB0 ) {const multiBM* mB=dynamic_cast<const multiBM*> ( mB0 ); beta=mB->beta;} 
     232        void bayes ( const vec &dt ) { 
     233                if ( frg<1.0 ) {beta*=frg;last_lognc=est.lognc();} 
     234                beta+=dt; 
     235                if ( evalll ) {ll=est.lognc()-last_lognc;} 
     236        } 
     237        double logpred ( const vec &dt ) const { 
     238                eDirich pred ( est ); 
     239                vec &beta = pred._beta(); 
     240                 
     241                double lll; 
     242                if ( frg<1.0 ) 
     243                        {beta*=frg;lll=pred.lognc();} 
     244                else 
     245                        if ( evalll ) {lll=last_lognc;} 
     246                        else{lll=pred.lognc();} 
     247 
     248                beta+=dt; 
     249                return pred.lognc()-lll; 
     250        } 
     251        void flatten (BMEF* B ) { 
     252                eDirich* E=dynamic_cast<eDirich*>(B); 
     253                // sum(beta) should be equal to sum(B.beta) 
     254                const vec &Eb=E->_beta(); 
     255                est.pow ( sum(beta)/sum(Eb) ); 
     256                if(evalll){last_lognc=est.lognc();} 
     257        } 
     258        const epdf& _epdf() const {return est;}; 
     259        //!access funct 
    148260}; 
    149261 
    … …  
    151263 \brief Gamma posterior density 
    152264 
    153  Multvariate Gamma density as product of independent univariate densities. 
     265 Multivariate Gamma density as product of independent univariate densities. 
    154266 \f[ 
    155267 f(x|\alpha,\beta) = \prod f(x_i|\alpha_i,\beta_i) 
    … …  
    213325        double evalpdflog ( const vec &val ) const  {return lnk;} 
    214326        vec sample() const { 
    215                 vec smp ( rv.count() );  
    216                 #pragma omp critical 
     327                vec smp ( rv.count() ); 
     328#pragma omp critical 
    217329                UniRNG.sample_vector ( rv.count(),smp ); 
    218                 return low+elem_mult(distance,smp); 
     330                return low+elem_mult ( distance,smp ); 
    219331        } 
    220332        //! set values of \c low and \c high 
    … …  
    408520double enorm<sq_T>::evalpdflog ( const vec &val ) const { 
    409521        // 1.83787706640935 = log(2pi) 
    410         return  -0.5* (  +R.invqform ( mu-val ) ) - lognc(); 
     522        return  -0.5* ( +R.invqform ( mu-val ) ) - lognc(); 
    411523}; 
    412524 
    … …  
    414526inline double enorm<sq_T>::lognc () const { 
    415527        // 1.83787706640935 = log(2pi) 
    416         return -0.5* ( R.cols() * 1.83787706640935 +R.logdet()); 
    417 }; 
    418  
    419 template<class sq_T> 
    420 mlnorm<sq_T>::mlnorm ( RV &rv0,RV &rvc0 ) :mEF ( rv0,rvc0 ),epdf ( rv0 ),A ( rv0.count(),rv0.count() ),_mu(epdf._mu()) { ep =&epdf; 
     528        return -0.5* ( R.cols() * 1.83787706640935 +R.logdet() ); 
     529}; 
     530 
     531template<class sq_T> 
     532mlnorm<sq_T>::mlnorm ( RV &rv0,RV &rvc0 ) :mEF ( rv0,rvc0 ),epdf ( rv0 ),A ( rv0.count(),rv0.count() ),_mu ( epdf._mu() ) { 
     533        ep =&epdf; 
    421534} 
    422535 

  • pmsm/pmsm_mix.cpp

    r162 r170  
    8383        M.set_est ( evolQR._epdf() ); 
    8484 
    85         epdf& Efix_ep = Efix._epdf(); 
    86         epdf& M_ep = M._epdf(); 
     85        const epdf& Efix_ep = Efix._epdf(); 
     86        const epdf& M_ep = M._epdf(); 
    8787         
    8888        //LOG 

  • pmsm/pmsm_sim.cpp

    r162 r170  
    7272        // 
    7373 
    74         epdf& KFEep = KFE._epdf(); 
    75         epdf& Mep = M._epdf(); 
     74        const epdf& KFEep = KFE._epdf(); 
     75        const epdf& Mep = M._epdf(); 
    7676 
    7777        mat Xt=zeros ( Ndat ,9 ); //true state from simulator 

  • pmsm/pmsm_sim2.cpp

    r162 r170  
    107107        evolQ.set_parameters ( 100000.0, Qdiag, 0.9999 ); 
    108108        // 
    109         epdf& KFEep = KFE._epdf(); 
    110         epdf& Mep = M._epdf(); 
     109        const epdf& KFEep = KFE._epdf(); 
     110        const epdf& Mep = M._epdf(); 
    111111 
    112112        int X_log = L.add(rx,"X"); 

  • pmsm/pmsm_unkQ.cpp

    r72 r170  
    6262        } 
    6363 
    64         epdf& KFEep = KFE._epdf(); 
     64        const epdf& KFEep = KFE._epdf(); 
    6565 
    6666        //simulator values 

  • pmsm/pmsm_unkQpf.cpp

    r162 r170  
    6262        MPF<EKF_unQ> M ( rx,rQ,evolQ,evolQ,100,KFE ); 
    6363 
    64         epdf& KFEep = KFE._epdf(); 
    65         epdf& Mep = M._epdf(); 
     64        const epdf& KFEep = KFE._epdf(); 
     65        const epdf& Mep = M._epdf(); 
    6666        // initialize 
    6767        evolQ.set_parameters ( 1.0 ); //sigma = 1/10 mu 
    6868        evolQ.condition ( "0.5 0.5" ); 
    69         epdf& pfinit=evolQ._epdf(); 
     69        const epdf& pfinit=evolQ._epdf(); 
    7070        M.set_est ( pfinit ); 
    7171        evolQ.set_parameters ( 1000.0 ); //sigma = 1/10 mu 

  • pmsm/sim_var.cpp

    r162 r170  
    7272        Edi.set_parameters ( &fxu,&hxu,Q,R); 
    7373 
    74         epdf& Efix_ep = Efix._epdf(); 
    75         epdf& Eop_ep = Eop._epdf(); 
    76         epdf& Edi_ep = Edi._epdf(); 
     74        const epdf& Efix_ep = Efix._epdf(); 
     75        const epdf& Eop_ep = Eop._epdf(); 
     76        const epdf& Edi_ep = Edi._epdf(); 
    7777         
    7878        //LOG 

  • pmsm/sim_var_arx.cpp

    r162 r170  
    4545        ARX Ar_a ( rgr,V0,nu0 ,0.95 ); 
    4646        ARX Ar_b ( rgr,V0,nu0 ,0.95 ); 
    47         epdf& pA= Ar_a._epdf(); 
    48         epdf& pB= Ar_b._epdf(); 
     47        const epdf& pA= Ar_a._epdf(); 
     48        const epdf& pB= Ar_b._epdf(); 
    4949 
    5050        RV rta ( "{th_a }",vec_1 ( rglen ) ); 

  • tests/arx_test.cpp

    r162 r170  
    1111        vec th("0.8 -0.3 0.4 0.01"); 
    1212        int ord=th.length(); 
    13         double sqr=1; 
     13        double sqr=0.1; 
    1414         
    1515        //Test constructor 
    16         mat V0 = 0.001*eye(ord+1); V0(0.0)*= 10; // 
     16        mat V0 = 0.00001*eye(ord+1); V0(0.0)*= 10; // 
    1717        double nu0 = ord+1; 
    1818         
    19         RV thr("{theta_and_r }"); 
     19        RV thr("{theta_and_r }",vec_1(ord+1)); 
    2020        ARX Ar(thr, V0, nu0); 
    21         epdf& Ar_ep = Ar._epdf(); 
    22          
     21        const epdf& Ar_ep = Ar._epdf(); 
     22                                 
    2323        //Test estimation 
    24         int ndat = 1000; 
     24        int ndat = 10000; 
    2525        int t,j; 
    2626        vec Yt(ndat); 

  • tests/test0.cpp

    r162 r170  
    44#include "../bdm/math/chmat.h" 
    55 
     6#include <vector> 
    67using namespace itpp; 
    78 

  • tests/testKF.cpp

    r162 r170  
    5757        KF.set_parameters(A,B,C,D,chmat(R),chmat(Q)); 
    5858        KF.set_est(mu0,chmat(P0) ); //prediction! 
    59         epdf& KFep = KF._epdf(); 
     59        const epdf& KFep = KF._epdf(); 
    6060        mat Xt(dimx,Ndat);       
    6161        Xt.set_col( 0,KFep.mean() ); 
    … …  
    6666        KFf.set_parameters(A,B,C,D,ldmat(R),ldmat(Q)); 
    6767        KFf.set_est(mu0,ldmat(P0) ); 
    68         epdf& KFfep = KFf._epdf(); 
     68        const   epdf& KFfep = KFf._epdf(); 
    6969        mat Xtf(dimx,Ndat);      
    7070        Xtf.set_col( 0,KFfep.mean() ); 
    … …  
    8282        KFE.set_parameters(&fxu,&hxu,Q,R); 
    8383        KFE.set_est(mu0,chmat(P0)); 
    84         epdf& KFEep = KFE._epdf(); 
     84        const epdf& KFEep = KFE._epdf(); 
    8585        mat XtE(dimx,Ndat);      
    8686        XtE.set_col( 0,KFEep.mean() ); 

  • tests/testKF_QR.cpp

    r162 r170  
    6565        epdf& pfinit=evolQR._epdf(); 
    6666        KF_QR.set_est(pfinit); 
    67         epdf& mpost=KF_QR._epdf(); 
    68         epdf& mposttr=KFtr._epdf(); 
     67        const epdf& mpost=KF_QR._epdf(); 
     68        const epdf& mposttr=KFtr._epdf(); 
    6969 
    7070        XQRt.set_col( 0,mpost.mean()); 

  • tests/testSmp.cpp

    r162 r170  
    9999        cout << "======= EProd ======== " << endl; 
    100100        // we have to change eG.rv to y 
    101         eG._rv()= y;  
     101        egamma eGy(x); 
     102        eGy.set_parameters(a,b); 
    102103        //create array  
    103104        Array<mpdf*> A(2); 
    104105        mepdf meN(eN); 
    105         mepdf meG(eG); 
     106        mepdf meG(eGy); 
    106107        A(0) = &meN; 
    107108        A(1) = &meG; 
    … …  
    114115        vec v0=vec(0); 
    115116        Smp = eP.samplecond(v0,lik,N); 
    116         disp(concat(eN.mean(),eG.mean()), epV,Smp); 
     117        disp(concat(eN.mean(),eGy.mean()), epV,Smp); 
    117118         
    118119        //Exit program: