Changeset 283 for bdm

Show
Ignore:
Timestamp:
02/24/09 14:14:01 (15 years ago)
Author:
smidl
Message:

get rid of BMcond + adaptation in doprava/

Location:
bdm
Files:
13 modified

Legend:

Unmodified
Added
Removed

  • bdm/estim/arx.cpp

    r270 r283  
    4444} 
    4545 
    46 ARX* ARX::_copy_ ( ) { 
     46ARX* ARX::_copy_ ( ) const { 
    4747        ARX* Tmp=new ARX ( *this ); 
    4848        return Tmp; 

  • bdm/estim/arx.h

    r271 r283  
    5555        ARX ( const double frg0=1.0 ) : BMEF ( frg0 ),est (), V ( est._V() ), nu ( est._nu() ) {}; 
    5656        ARX ( const ARX &A0 ) : BMEF (),est ( A0.est ), V ( est._V() ), nu ( est._nu() ) {}; 
    57         ARX* _copy_(); 
     57        ARX* _copy_() const; 
    5858        void set_parameters ( double frg0 ) {frg=frg0;} 
    5959        void set_statistics ( int dimx0, const ldmat V0, double nu0=-1.0 ) {est.set_parameters ( dimx0,V0,nu0 );last_lognc=est.lognc();dimx=dimx0;} 

  • bdm/estim/ekf_templ.h

    r279 r283  
    1919 
    2020//!Extended Kalman filter with unknown \c Q and \c R 
    21 class EKFful_unQR : public EKFfull , public BMcond { 
     21class EKFful_unQR : public EKFfull { 
    2222public: 
    23         //! Default constructor 
    24         EKFful_unQR ( ) :EKFfull ( ),BMcond ( ) {}; 
    2523        void condition ( const vec &QR0 ) { 
    2624                Q=diag(QR0(0,dimx-1)); 
    … …  
    3028 
    3129//!Extended Kalman filter in Choleski form with unknown \c Q 
    32 class EKFCh_unQ : public EKFCh , public BMcond { 
     30class EKFCh_unQ : public EKFCh { 
    3331public: 
    34         //! Default constructor 
    35         EKFCh_unQ (  ) :EKFCh ( ),BMcond ( ) {}; 
    3632        void condition ( const vec &Q0 ) { 
    3733                Q.setD ( Q0,0 ); 
    … …  
    4238 
    4339//!Extended Kalman filter with unknown parameters in \c IM 
    44 class EKFCh_cond : public EKFCh , public BMcond { 
     40class EKFCh_cond : public EKFCh { 
    4541        public: 
    46         //! Default constructor 
    47                 EKFCh_cond ( ) :EKFCh ( ),BMcond ( ) {}; 
    4842                void condition ( const vec &val ) { 
    4943                        pfxu->condition ( val ); 

  • bdm/estim/libKF.cpp

    r279 r283  
    6969 
    7070        dimx = pfxu0->_dimx(); 
    71         dimy = phxu0->_dimy(); 
     71        dimy = phxu0->dimension(); 
    7272        dimu = pfxu0->_dimu(); 
    7373 
    … …  
    169169        phxu = phxu0; 
    170170 
    171         dimx = pfxu0->_dimy(); 
    172         dimy = phxu0->_dimy(); 
     171        dimx = pfxu0->dimension(); 
     172        dimy = phxu0->dimension(); 
    173173        dimu = pfxu0->_dimu(); 
    174174        // set size of mu - just for constant terms of A and C  
    … …  
    249249} 
    250250 
    251 void KFcondQR::condition ( const vec &QR ) { 
    252         it_assert_debug ( QR.length() == ( dimx+dimy ),"KFcondRQ: conditioning by incompatible vector" ); 
    253  
    254         Q.setD ( QR ( 0, dimx-1 ) ); 
    255         R.setD ( QR ( dimx, -1 ) ); 
    256 }; 
    257  
    258 void KFcondR::condition ( const vec &R0 ) { 
    259         it_assert_debug ( R0.length() == ( dimy ),"KFcondR: conditioning by incompatible vector" ); 
    260  
    261         R.setD ( R0 ); 
    262 }; 
    263  
    264  
    265 } 
     251} 

  • bdm/estim/libKF.h

    r279 r283  
    1919#include "../math/chmat.h" 
    2020 
    21 namespace bdm{ 
     21namespace bdm { 
    2222 
    2323/*! 
    … …  
    4949        friend std::ostream &operator<< ( std::ostream &os, const KalmanFull &kf ); 
    5050        //! For EKFfull; 
    51         KalmanFull(){}; 
     51        KalmanFull() {}; 
    5252}; 
    5353 
    … …  
    5656* \brief Kalman filter with covariance matrices in square root form. 
    5757 
    58 Parameter evolution model:\f[ x_t = A x_{t-1} + B u_t + Q^{1/2} e_t \f]  
     58Parameter evolution model:\f[ x_t = A x_{t-1} + B u_t + Q^{1/2} e_t \f] 
    5959Observation model: \f[ y_t = C x_{t-1} + C u_t + Q^{1/2} w_t. \f] 
    60 Where $e_t$ and $w_t$ are independent vectors Normal(0,1)-distributed disturbances.  
     60Where $e_t$ and $w_t$ are independent vectors Normal(0,1)-distributed disturbances. 
    6161*/ 
    6262template<class sq_T> 
    … …  
    7777        mat A; 
    7878        //! Matrix B 
    79         mat B;  
     79        mat B; 
    8080        //! Matrix C 
    8181        mat C; 
    … …  
    112112        //! Set estimate values, used e.g. in initialization. 
    113113        void set_est ( const vec &mu0, const sq_T &P0 ) { 
    114                 sq_T pom(dimy); 
     114                sq_T pom ( dimy ); 
    115115                est.set_parameters ( mu0,P0 ); 
    116                 P0.mult_sym(C,pom); 
     116                P0.mult_sym ( C,pom ); 
    117117                fy.set_parameters ( C*mu0, pom ); 
    118118        }; 
    … …  
    133133Trivial example: 
    134134\include kalman_simple.cpp 
    135   
    136 */ 
    137 class KalmanCh : public Kalman<chmat>{ 
     135 
     136*/ 
     137class KalmanCh : public Kalman<chmat> { 
    138138protected: 
    139139//! pre array (triangular matrix) 
    140 mat preA; 
     140        mat preA; 
    141141//! post array (triangular matrix) 
    142 mat postA; 
    143  
    144 public: 
    145         //! Default constructor 
    146         KalmanCh ():Kalman<chmat>(),preA(),postA(){}; 
     142        mat postA; 
     143 
     144public: 
     145        //! copy constructor 
     146        BM* _copy_() const { 
     147                KalmanCh* K=new KalmanCh; 
     148                K->set_parameters ( A,B,C,D,Q,R ); 
     149                K->set_statistics ( _mu,_P ); 
     150                return K; 
     151        } 
    147152        //! Set parameters with check of relevance 
    148153        void set_parameters ( const mat &A0,const mat &B0,const mat &C0,const mat &D0,const chmat &Q0,const chmat &R0 ); 
    … …  
    150155                est.set_parameters ( mu0,P0 ); 
    151156        }; 
    152          
    153          
     157 
     158 
    154159        /*!\brief  Here dt = [yt;ut] of appropriate dimensions 
    155          
     160 
    156161        The following equality hold::\f[ 
    157 \left[\begin{array}{cc} 
    158 R^{0.5}\\ 
    159 P_{t|t-1}^{0.5}C' & P_{t|t-1}^{0.5}CA'\\ 
    160  & Q^{0.5}\end{array}\right]<\mathrm{orth.oper.}>=\left[\begin{array}{cc} 
    161 R_{y}^{0.5} & KA'\\ 
    162  & P_{t+1|t}^{0.5}\\ 
    163 \\\end{array}\right]\f] 
    164  
    165 Thus this object evaluates only predictors! Not filtering densities. 
     162        \left[\begin{array}{cc} 
     163        R^{0.5}\\ 
     164        P_{t|t-1}^{0.5}C' & P_{t|t-1}^{0.5}CA'\\ 
     165        & Q^{0.5}\end{array}\right]<\mathrm{orth.oper.}>=\left[\begin{array}{cc} 
     166        R_{y}^{0.5} & KA'\\ 
     167        & P_{t+1|t}^{0.5}\\ 
     168        \\\end{array}\right]\f] 
     169 
     170        Thus this object evaluates only predictors! Not filtering densities. 
    166171        */ 
    167172        void bayes ( const vec &dt ); 
    … …  
    179184        //! Observation Model h(x,u) 
    180185        diffbifn* phxu; 
    181          
    182         enorm<fsqmat> E;  
     186 
     187        enorm<fsqmat> E; 
    183188public: 
    184189        //! Default constructor 
    … …  
    189194        void bayes ( const vec &dt ); 
    190195        //! set estimates 
    191         void set_est (vec mu0, mat P0){mu=mu0;P=P0;}; 
     196        void set_est ( vec mu0, mat P0 ) {mu=mu0;P=P0;}; 
    192197        //!dummy! 
    193         const epdf& posterior()const{return E;}; 
    194         const enorm<fsqmat>* _e()const{return &E;}; 
    195         const mat _R(){return P;} 
     198        const epdf& posterior() const{return E;}; 
     199        const enorm<fsqmat>* _e() const{return &E;}; 
     200        const mat _R() {return P;} 
    196201}; 
    197202 
    … …  
    210215        //! Default constructor 
    211216        EKF ( RV rvx, RV rvy, RV rvu ); 
     217        //! copy constructor 
     218        EKF<sq_T>* _copy_() const { return new EKF<sq_T>(this); } 
    212219        //! Set nonlinear functions for mean values and covariance matrices. 
    213220        void set_parameters ( diffbifn* pfxu, diffbifn* phxu, const sq_T Q0, const sq_T R0 ); 
    … …  
    223230 
    224231class EKFCh : public KalmanCh { 
    225         protected: 
     232protected: 
    226233        //! Internal Model f(x,u) 
    227234        diffbifn* pfxu; 
    … …  
    229236        diffbifn* phxu; 
    230237public: 
     238        //! copy constructor duplicated - calls different set_parameters 
     239        BM* _copy_() const { 
     240                EKFCh* E=new EKFCh; 
     241                E->set_parameters ( pfxu,phxu,Q,R ); 
     242                E->set_statistics ( _mu,_P ); 
     243                return E; 
     244        } 
    231245        //! Set nonlinear functions for mean values and covariance matrices. 
    232246        void set_parameters ( diffbifn* pfxu, diffbifn* phxu, const chmat Q0, const chmat R0 ); 
    … …  
    239253*/ 
    240254 
    241 class KFcondQR : public Kalman<ldmat>, public BMcond { 
     255class KFcondQR : public Kalman<ldmat> { 
    242256//protected: 
    243257public: 
     258        void condition ( const vec &QR ) { 
     259                it_assert_debug ( QR.length() == ( dimx+dimy ),"KFcondRQ: conditioning by incompatible vector" ); 
     260 
     261                Q.setD ( QR ( 0, dimx-1 ) ); 
     262                R.setD ( QR ( dimx, -1 ) ); 
     263        }; 
     264}; 
     265 
     266/*! 
     267\brief Kalman Filter with conditional diagonal matrices R and Q. 
     268*/ 
     269 
     270class KFcondR : public Kalman<ldmat> { 
     271//protected: 
     272public: 
    244273        //!Default constructor 
    245         KFcondQR ( ) : Kalman<ldmat> ( ),BMcond ( ) {}; 
    246  
    247         void condition ( const vec &RQ ); 
    248 }; 
    249  
    250 /*! 
    251 \brief Kalman Filter with conditional diagonal matrices R and Q. 
    252 */ 
    253  
    254 class KFcondR : public Kalman<ldmat>, public BMcond { 
    255 //protected: 
    256 public: 
    257         //!Default constructor 
    258         KFcondR ( ) : Kalman<ldmat> ( ),BMcond ( ) {}; 
    259  
    260         void condition ( const vec &R ); 
     274        KFcondR ( ) : Kalman<ldmat> ( ) {}; 
     275 
     276        void condition ( const vec &R0 ) { 
     277                it_assert_debug ( R0.length() == ( dimy ),"KFcondR: conditioning by incompatible vector" ); 
     278 
     279                R.setD ( R0 ); 
     280        }; 
     281 
    261282}; 
    262283 
    … …  
    267288                dimx ( K0.dimx ), dimy ( K0.dimy ),dimu ( K0.dimu ), 
    268289                A ( K0.A ), B ( K0.B ), C ( K0.C ), D ( K0.D ), 
    269                 Q(K0.Q), R(K0.R), 
    270                 est ( K0.est ), fy ( K0.fy ), _yp(fy._mu()),_Ry(fy._R()), _mu(est._mu()), _P(est._R()) { 
     290                Q ( K0.Q ), R ( K0.R ), 
     291                est ( K0.est ), fy ( K0.fy ), _yp ( fy._mu() ),_Ry ( fy._R() ), _mu ( est._mu() ), _P ( est._R() ) { 
    271292 
    272293// copy values in pointers 
    … …  
    279300 
    280301template<class sq_T> 
    281 Kalman<sq_T>::Kalman ( ) : BM (), est ( ), fy (),  _yp(fy._mu()), _Ry(fy._R()), _mu(est._mu()), _P(est._R()) { 
     302Kalman<sq_T>::Kalman ( ) : BM (), est ( ), fy (),  _yp ( fy._mu() ), _Ry ( fy._R() ), _mu ( est._mu() ), _P ( est._R() ) { 
    282303}; 
    283304 
    … …  
    287308        dimy = C0.rows(); 
    288309        dimu = B0.cols(); 
    289          
     310 
    290311        it_assert_debug ( A0.cols() ==dimx, "Kalman: A is not square" ); 
    291312        it_assert_debug ( B0.rows() ==dimx, "Kalman: B is not compatible" ); 
    … …  
    307328        it_assert_debug ( dt.length() == ( dimy+dimu ),"KalmanFull::bayes wrong size of dt" ); 
    308329 
    309         sq_T iRy(dimy); 
     330        sq_T iRy ( dimy ); 
    310331        vec u = dt.get ( dimy,dimy+dimu-1 ); 
    311332        vec y = dt.get ( 0,dimy-1 ); 
    … …  
    328349        sq_T pom ( ( int ) Pfull.rows() ); 
    329350        iRy.mult_sym_t ( C*Pfull,pom ); 
    330         (_P ) -= pom; // P = P -PC'iRy*CP; 
    331         (_yp ) = C* _mu  +D*u; //y prediction 
    332         (_mu ) += _K* ( y- _yp ); 
     351        ( _P ) -= pom; // P = P -PC'iRy*CP; 
     352        ( _yp ) = C* _mu  +D*u; //y prediction 
     353        ( _mu ) += _K* ( y- _yp ); 
    333354 
    334355 
    … …  
    340361 
    341362}; 
    342   
     363 
    343364 
    344365 
    … …  
    369390        it_assert_debug ( dt.length() == ( dimy+dimu ),"KalmanFull::bayes wrong size of dt" ); 
    370391 
    371         sq_T iRy(dimy,dimy); 
     392        sq_T iRy ( dimy,dimy ); 
    372393        vec u = dt.get ( dimy,dimy+dimu-1 ); 
    373394        vec y = dt.get ( 0,dimy-1 ); 
    … …  
    393414        sq_T pom ( ( int ) Pfull.rows() ); 
    394415        iRy.mult_sym_t ( C*Pfull,pom ); 
    395         (_P ) -= pom; // P = P -PC'iRy*CP; 
     416        ( _P ) -= pom; // P = P -PC'iRy*CP; 
    396417        _yp = phxu->eval ( _mu,u ); //y prediction 
    397418        ( _mu ) += _K* ( y-_yp ); 

  • bdm/estim/libPF.cpp

    r278 r283  
    11#include "libPF.h" 
    22 
    3 namespace bdm{ 
     3namespace bdm { 
    44 
    55using std::endl; 
    … …  
    1313        for ( i=0;i<n;i++ ) { 
    1414                //generate new samples from paramater evolution model; 
    15                 _samples ( i ) = par->samplecond ( _samples ( i )); 
    16                 lls(i )= par->_e()->evallog(_samples(i)); 
     15                _samples ( i ) = par->samplecond ( _samples ( i ) ); 
     16                lls ( i ) = par->_e()->evallog ( _samples ( i ) ); 
    1717                lls ( i ) *= obs->evallogcond ( dt,_samples ( i ) ); 
    1818 
    … …  
    3030        _w ( i ) /=sum; //? 
    3131 
    32         ind = est.resample(); 
     32        ind = est.resample(resmethod); 
    3333 
    3434} 
    3535 
    36 void PF::set_est ( const epdf &epdf0 ) { 
    37         int i; 
     36// void PF::set_est ( const epdf &epdf0 ) { 
     37//      int i; 
     38//  
     39//      for ( i=0;i<n;i++ ) { 
     40//              _samples ( i ) = epdf0.sample(); 
     41//      } 
     42// } 
    3843 
    39         for ( i=0;i<n;i++ ) { 
    40                 _samples ( i ) = epdf0.sample(); 
    41         } 
    42 } 
    4344 
    4445} 

  • bdm/estim/libPF.h

    r281 r283  
    4040        mpdf *obs; 
    4141 
     42        //! which resampling method will be used 
     43        RESAMPLING_METHOD resmethod; 
     44 
    4245        //! \name Options 
    4346        //!@{ 
    4447 
    45         //! Log resampling times 
    46         bool opt_L_res; 
    4748        //! Log all samples 
    4849        bool opt_L_smp; 
     50        //! Log all samples 
     51        bool opt_L_wei; 
    4952        //!@} 
    5053 
    … …  
    5255        //! \name Constructors 
    5356        //!@{ 
    54         PF ( ) :est(), _w ( est._w() ),_samples ( est._samples() ) {}; 
    55         PF ( mpdf *par0, mpdf *obs0, epdf *epdf0, int n0 ) : 
    56                         est ( ),_w ( est._w() ),_samples ( est._samples() ) 
    57         { set_parameters ( par0,obs0,n0 ); set_statistics ( ones ( n0 ),epdf0 ); }; 
    58         void set_parameters ( mpdf *par0, mpdf *obs0, int n0 ) 
    59         { par = par0; obs=obs0; n=n0; est.set_n ( n );}; 
    60         void set_statistics ( const vec w0, epdf *epdf0 ) {est.set_parameters ( w0,epdf0 );}; 
     57        PF ( ) :est(), _w ( est._w() ),_samples ( est._samples() ), opt_L_smp ( false ), opt_L_wei ( false ) {LIDs.set_size ( 5 );}; 
     58        /*      PF ( mpdf *par0, mpdf *obs0, epdf *epdf0, int n0 ) : 
     59                                est ( ),_w ( est._w() ),_samples ( est._samples() ),opt_L_smp(false), opt_L_wei(false) 
     60                { set_parameters ( par0,obs0,n0 ); set_statistics ( ones ( n0 ),epdf0 ); };*/ 
     61        void set_parameters ( mpdf *par0, mpdf *obs0, int n0, RESAMPLING_METHOD rm=SYSTEMATIC ) 
     62        { par = par0; obs=obs0; n=n0; resmethod= rm;}; 
     63        void set_statistics ( const vec w0, epdf *epdf0 ) {est.set_statistics ( w0,epdf0 );}; 
    6164        //!@} 
    6265        //! Set posterior density by sampling from epdf0 
    63         void set_est ( const epdf &epdf0 ); 
    64         void set_options(const string &opt){ 
    65                 opt_L_res= ( s.find ( "logres" ) !=string::npos ); 
    66                 opt_L_smp= ( s.find ( "logsmp" ) !=string::npos );               
     66//      void set_est ( const epdf &epdf0 ); 
     67        void set_options ( const string &opt ) { 
     68                BM::set_options(opt); 
     69                opt_L_wei= ( opt.find ( "logweights" ) !=string::npos ); 
     70                opt_L_smp= ( opt.find ( "logsamples" ) !=string::npos ); 
    6771        } 
    6872        void bayes ( const vec &dt ); 
    … …  
    7882 
    7983template<class BM_T> 
    80  
    8184class MPF : public PF { 
    82         BM_T* Bms[10000]; 
     85        Array<BM_T*> BMs; 
    8386 
    8487        //! internal class for MPDF providing composition of eEmp with external components 
    … …  
    9497                                Coms ( _w.length() ) { 
    9598                }; 
     99                //! read statistics from MPF 
     100                void read_statistics ( Array<BM_T*> &A ) { 
     101                        dim = E.dimension() +A ( 0 )->posterior().dimension(); 
     102                        for ( int i=0; i<_w.length() ;i++ ) {Coms ( i ) = A ( i )->_e();} 
     103                } 
     104                //! needed in resampling 
    96105                void set_elements ( int &i, double wi, const epdf* ep ) 
    97106                {_w ( i ) =wi; Coms ( i ) =ep;}; 
    98107 
    99                 void set_n ( int n ) {E.set_n ( n ); Coms.set_length ( n );} 
     108                void set_parameters ( int n ) { 
     109                        E.set_parameters ( n, false ); 
     110                        Coms.set_length ( n ); 
     111                } 
    100112                vec mean() const { 
    101113                        // ugly 
    … …  
    114126                        return concat ( E.variance(),pom2-pow ( pom,2 ) ); 
    115127                } 
     128                void qbounds ( vec &lb, vec &ub, double perc=0.95 ) const { 
     129                        //bounds on particles 
     130                        vec lbp; 
     131                        vec ubp; 
     132                        E.qbounds ( lbp,ubp ); 
     133 
     134                        //bounds on Components 
     135                        int dimC=Coms ( 0 )->dimension(); 
     136                        int j; 
     137                        // temporary 
     138                        vec lbc(dimC); 
     139                        vec ubc(dimC); 
     140                        // minima and maxima 
     141                        vec Lbc(dimC); 
     142                        vec Ubc(dimC); 
     143                        Lbc = std::numeric_limits<double>::infinity(); 
     144                        Ubc = -std::numeric_limits<double>::infinity(); 
     145 
     146                        for ( int i=0;i<_w.length();i++ ) { 
     147                                // check Coms 
     148                                Coms ( i )->qbounds ( lbc,ubc ); 
     149                                for ( j=0;j<dimC; j++ ) { 
     150                                        if ( lbc ( j ) <Lbc ( j ) ) {Lbc ( j ) =lbc ( j );} 
     151                                        if ( ubc ( j ) >Ubc ( j ) ) {Ubc ( j ) =ubc ( j );} 
     152                                } 
     153                        } 
     154                        lb=concat(lbp,Lbc); 
     155                        ub=concat(ubp,Ubc); 
     156                } 
    116157 
    117158                vec sample() const {it_error ( "Not implemented" );return 0;} 
    … …  
    125166        //! Log means of BMs 
    126167        bool opt_L_mea; 
    127          
     168 
    128169public: 
    129170        //! Default constructor. 
    130         MPF ( mpdf *par0, mpdf *obs0, int n, const BM_T &BMcond0 ) : PF (), jest ( est ) { 
    131                 PF::set_parameters ( par0,obs0,n ); 
    132                 jest.set_n ( n ); 
    133                 // 
    134                 //TODO test if rv and BMcond.rv are compatible. 
    135 //              rv.add ( rvlin ); 
    136                 // 
    137  
    138                 if ( n>10000 ) {it_error ( "increase 10000 here!" );} 
    139  
    140                 for ( int i=0;i<n;i++ ) { 
    141                         Bms[i] = new BM_T ( BMcond0 ); //copy constructor 
    142                         const epdf& pom=Bms[i]->posterior(); 
    143                         jest.set_elements ( i,1.0/n,&pom ); 
    144                 } 
     171        MPF () : PF (), jest ( est ) {}; 
     172        void set_parameters ( mpdf *par0, mpdf *obs0, int n0, RESAMPLING_METHOD rm=SYSTEMATIC ) { 
     173                PF::set_parameters ( par0, obs0, n0, rm ); 
     174                jest.set_parameters ( n0 );//duplication of rm 
     175                BMs.set_length ( n0 ); 
     176        } 
     177        void set_statistics ( epdf *epdf0, const BM_T* BMcond0 ) { 
     178 
     179                PF::set_statistics ( ones ( n ) /n, epdf0 ); 
     180                // copy 
     181                for ( int i=0;i<n;i++ ) { BMs ( i ) = new BM_T ( *BMcond0 ); BMs ( i )->condition ( _samples ( i ) );} 
     182 
     183                jest.read_statistics ( BMs ); 
     184                //options 
    145185        }; 
    146  
    147         ~MPF() { 
    148         } 
    149186 
    150187        void bayes ( const vec &dt ); 
    151188        const epdf& posterior() const {return jest;} 
    152189        const epdf* _e() const {return &jest;} //Fixme: is it useful? 
    153         //! Set postrior of \c rvc to samples from epdf0. Statistics of Bms are not re-computed! Use only for initialization! 
    154         void set_est ( const epdf& epdf0 ) { 
    155                 PF::set_est ( epdf0 );  // sample params in condition 
    156                 // copy conditions to BMs 
    157  
    158                 for ( int i=0;i<n;i++ ) {Bms[i]->condition ( _samples ( i ) );} 
    159         } 
    160         void set_options(const string &opt){ 
    161                 PF:set_options(opt); 
    162                 opt_L_mea = ( s.find ( "logmeans" ) !=string::npos ); 
    163         } 
    164          
     190        //! Set postrior of \c rvc to samples from epdf0. Statistics of BMs are not re-computed! Use only for initialization! 
     191        /*      void set_est ( const epdf& epdf0 ) { 
     192                        PF::set_est ( epdf0 );  // sample params in condition 
     193                        // copy conditions to BMs 
     194 
     195                        for ( int i=0;i<n;i++ ) {BMs(i)->condition ( _samples ( i ) );} 
     196                }*/ 
     197        void set_options ( const string &opt ) { 
     198                PF::set_options ( opt ); 
     199                opt_L_mea = ( opt.find ( "logmeans" ) !=string::npos ); 
     200        } 
     201 
    165202        //!Access function 
    166         BM* _BM ( int i ) {return Bms[i];} 
     203        BM* _BM ( int i ) {return BMs ( i );} 
    167204}; 
    168205 
    … …  
    180217                _samples ( i ) = par->samplecond ( _samples ( i ) ); 
    181218                llsP ( i ) = par->_e()->evallog ( _samples ( i ) ); 
    182                 Bms[i]->condition ( _samples ( i ) ); 
    183                 Bms[i]->bayes ( dt ); 
    184                 lls ( i ) = Bms[i]->_ll(); // lls above is also in proposal her must be lls(i) =, not +=!! 
     219                BMs ( i )->condition ( _samples ( i ) ); 
     220                BMs ( i )->bayes ( dt ); 
     221                lls ( i ) = BMs ( i )->_ll(); // lls above is also in proposal her must be lls(i) =, not +=!! 
    185222                if ( lls ( i ) >mlls ) mlls=lls ( i ); //find maximum likelihood (for numerical stability) 
    186223        } 
    … …  
    204241        double eff = 1.0/ ( _w*_w ); 
    205242        if ( eff < ( 0.3*n ) ) { 
    206                 ind = est.resample(); 
     243                ind = est.resample ( resmethod ); 
    207244                // Resample Bms! 
    208245 
    … …  
    215252                                // poor-man's solution: replicate constructor here 
    216253                                // copied from MPF::MPF 
    217                                 delete Bms[i]; 
    218                                 Bms[i] = new BM_T ( *Bms[ind ( i ) ] ); //copy constructor 
    219                                 const epdf& pom=Bms[i]->posterior(); 
     254                                delete BMs ( i ); 
     255                                BMs ( i ) = new BM_T ( *BMs ( ind ( i ) ) ); //copy constructor 
     256                                const epdf& pom=BMs ( i )->posterior(); 
    220257                                jest.set_elements ( i,1.0/n,&pom ); 
    221258                        } 

  • bdm/estim/merger.cpp

    r278 r283  
    3333        it_assert_debug ( rv.equal ( g0->_rv() ),"Incompatible g0" ); 
    3434        //Empirical density - samples 
    35         eSmp.set_parameters ( ones ( Ns ), g0 ); 
     35        eSmp.set_statistics ( ones ( Ns ), g0 ); 
    3636        Array<vec> &Smp = eSmp._samples(); //aux 
    3737        vec &w = eSmp._w(); //aux 

  • bdm/estim/merger.h

    r278 r283  
    7171        } 
    7272//! Set internal parameters used in approximation 
    73         void set_parameters ( double beta0, int Ns0, int Nc0 ) {beta=beta0;Ns=Ns0;Nc=Nc0;eSmp.set_n(Ns0,false);} 
     73        void set_parameters ( double beta0, int Ns0, int Nc0 ) {beta=beta0;Ns=Ns0;Nc=Nc0;eSmp.set_parameters(Ns0,false);} 
    7474//!Initialize the proposal density. This function must be called before merge()! 
    7575        void init() { ////////////// NOT FINISHED 

  • bdm/stat/libBM.h

    r281 r283  
    126126        int length() const {return len;} ; 
    127127        int id ( int at ) const{return ids ( at );}; 
    128         int size ( int at ) const {return RV_SIZES ( ids(at) );}; 
     128        int size ( int at ) const {return RV_SIZES ( ids ( at ) );}; 
    129129        int time ( int at ) const{return times ( at );}; 
    130         std::string name ( int at ) const {return RV_NAMES ( ids(at) );}; 
     130        std::string name ( int at ) const {return RV_NAMES ( ids ( at ) );}; 
    131131        void set_time ( int at, int time0 ) {times ( at ) =time0;}; 
    132132        //!@} 
    … …  
    204204 
    205205        //! access function 
    206         int _dimy() const{return dimy;} 
     206        int dimension() const{return dimy;} 
    207207}; 
    208208 
    … …  
    259259        //! return expected variance (not covariance!) 
    260260        virtual vec variance() const {it_error ( "not implemneted" );return vec ( 0 );}; 
     261        //! Lower and upper bounds of \c percentage % quantile, returns mean-2*sigma as default 
     262        virtual void qbounds ( vec &lb, vec &ub, double percentage=0.95 ) const { 
     263                vec mea=mean(); vec std=sqrt(variance());  
     264                lb = mea-2*std; ub=mea+2*std; 
     265        }; 
    261266        //!@} 
    262267 
    … …  
    592597/*! \brief Bayesian Model of a system, i.e. all uncertainty is modeled by probabilities. 
    593598 
     599This object represents exact or approximate evaluation of the Bayes rule: 
     600\f[ 
     601f(\theta_t | d_1,\ldots,d_t) = \frac{f(y_t|\theta_t,\cdot) f(\theta_t|d_1,\ldots,d_{t-1})}{f(y_t|d_1,\ldots,d_{t-1})} 
     602\f] 
     603 
     604Access to the resulting posterior density is via function \c posterior(). 
     605 
     606As a "side-effect" it also evaluates log-likelihood of the data, which can be accessed via function _ll(). 
     607It can also evaluate predictors of future values of \f$y_t\f$, see functions epredictor() and predictor(). 
     608 
     609Alternatively, it can evaluate posterior density conditioned by a known constant, \f$ c_t \f$: 
     610\f[ 
     611f(\theta_t | c_t, d_1,\ldots,d_t) \propto  f(y_t,\theta_t|c_t,\cdot, d_1,\ldots,d_{t-1}) 
     612\f] 
     613 
     614The value of \f$ c_t \f$ is set by function condition(). 
     615 
    594616*/ 
    595617 
    … …  
    606628        //! @{ 
    607629 
    608         BM () :ll ( 0 ),evalll ( true) {}; 
     630        BM () :ll ( 0 ),evalll ( true ), LIDs ( 3 ), opt_L_bounds ( false ) {}; 
    609631        BM ( const BM &B ) :  drv ( B.drv ), ll ( B.ll ), evalll ( B.evalll ) {} 
    610632        //! Copy function required in vectors, Arrays of BM etc. Have to be DELETED manually! 
    611         //! Prototype: \code BM* _copy_(){return new BM(*this);} \endcode 
    612         virtual BM* _copy_ () {return NULL;}; 
     633        //! Prototype: \code BM* _copy_() const {return new BM(*this);} \endcode 
     634        virtual BM* _copy_ () const {return NULL;}; 
    613635        //!@} 
    614636 
    … …  
    634656        //!@} 
    635657 
     658        //! \name Extension to conditional BM 
     659        //! This extension is useful e.g. in Marginalized Particle Filter (\ref bdm::MPF). 
     660        //! Alternatively, it can be used for automated connection to DS when the condition is observed 
     661        //!@{ 
     662 
     663        //! Name of extension variable 
     664        RV rvc; 
     665        //! access function 
     666        const RV& _rvc() const {return rvc;} 
     667 
     668        //! Substitute \c val for \c rvc. 
     669        virtual void condition ( const vec &val ) {it_error ( "Not implemented!" );}; 
     670 
     671        //!@} 
     672 
     673 
    636674        //! \name Access to attributes 
    637675        //!@{ 
    … …  
    646684        //!@} 
    647685 
    648 }; 
    649  
    650 /*! 
    651 \brief Conditional Bayesian Filter 
    652  
    653 Evaluates conditional filtering density \f$f(rv|rvc,data)\f$ for a given \c rvc which is specified in each step by calling function \c condition. 
    654  
    655 This is an interface class used to assure that certain BM has operation \c condition . 
    656  
    657 */ 
    658  
    659 class BMcond :public bdmroot { 
    660 protected: 
    661         //!dimension of the conditioning variable 
    662         int dimc; 
    663         //! Identificator of the conditioning variable 
    664         RV rvc; 
    665 public: 
    666         //! Substitute \c val for \c rvc. 
    667         virtual void condition ( const vec &val ) =0; 
    668         //! Default constructor 
    669         BMcond ( ) :rvc ( ) {}; 
    670         //! Destructor for future use 
    671         virtual ~BMcond() {}; 
    672         //! access function 
    673         const RV& _rvc() const {return rvc;} 
     686        //! \name Logging of results 
     687        //!@{ 
     688 
     689        //! Set boolean options from a string 
     690        void set_options ( const string &opt ) { 
     691                opt_L_bounds= ( opt.find ( "logbounds" ) !=string::npos ); 
     692        } 
     693        //! IDs of storages in loggers 
     694        ivec LIDs; 
     695 
     696        //! Option for logging bounds 
     697        bool opt_L_bounds; 
     698        //! Add all logged variables to a logger 
     699        void log_add ( logger *L, const string &name="" ) { 
     700                // internal 
     701                RV r; 
     702                if ( posterior().isnamed() ) {r=posterior()._rv();} 
     703                else{r=RV ( "est", posterior().dimension() );}; 
     704 
     705                // Add mean value 
     706                LIDs ( 0 ) =L->add ( r,name ); 
     707                if ( opt_L_bounds ) { 
     708                        LIDs ( 1 ) =L->add ( r,name+"_lb" ); 
     709                        LIDs ( 2 ) =L->add ( r,name+"_ub" ); 
     710                } 
     711        } 
     712        void logit ( logger *L ) { 
     713                L->logit ( LIDs ( 0 ), posterior().mean() ); 
     714                if ( opt_L_bounds ) { 
     715                        vec ub,lb; 
     716                        posterior().qbounds(lb,ub); 
     717                        L->logit ( LIDs ( 1 ), lb ); 
     718                        L->logit ( LIDs ( 2 ), ub ); 
     719                } 
     720        } 
     721        //!@} 
    674722}; 
    675723 

  • bdm/stat/libDS.h

    r271 r283  
    2828*/ 
    2929class MemDS : public DS { 
     30        protected: 
    3031        //! internal matrix of data 
    3132        mat Data; 
    … …  
    4546        void step(); 
    4647        //!Default constructor 
     48        MemDS () {}; 
    4749        MemDS ( mat &Dat, ivec &rowid, ivec &delays ); 
     50}; 
     51 
     52/*! Read Data Matrix from an IT file 
     53 
     54*/ 
     55class FileDS: public MemDS { 
     56 
     57public: 
     58        FileDS ( const string &fname, const string &varname ) :MemDS() { 
     59                it_file it ( fname ); 
     60                it << Name ( varname );  
     61                it >> Data; 
     62                time =0; 
     63                //rowid and delays are ignored 
     64        } 
     65        void getdata ( vec &dt ) { 
     66                it_assert_debug ( dt.length() ==Data.rows(),"" ); 
     67                dt = Data.get_col(time); 
     68        }; 
     69        void getdata ( vec &dt, const ivec &indeces ){ 
     70                it_assert_debug ( dt.length() ==indeces.length(),"" ); 
     71                vec tmp(indeces.length()); 
     72                tmp = Data.get_col(time); 
     73                dt = tmp(indeces); 
     74        }; 
     75        //! returns number of data in the file; 
     76        int ndat(){return Data.cols();} 
    4877}; 
    4978 
    … …  
    96125        { model.set_parameters ( Th0, mu0, sqR0 );}; 
    97126        //! Set 
    98         void set_drv(RV &yrv, RV &urv, RV &rrv){ 
     127        void set_drv ( RV &yrv, RV &urv, RV &rrv ) { 
    99128                Rrv = rrv; 
    100129                Urv = urv; 
    101                 dt_size = yrv._dsize()+urv._dsize(); 
    102                  
    103                 RV drv = concat(yrv,urv); 
     130                dt_size = yrv._dsize() +urv._dsize(); 
     131 
     132                RV drv = concat ( yrv,urv ); 
    104133                Drv = drv; 
    105134                int td = rrv.mint(); 
    106                 H.set_size(drv._dsize()*(-td+1)); 
    107                 U.set_size(Urv._dsize()); 
    108                 for (int i=-1;i>=td;i--){ 
    109                         drv.t(-1); 
    110                         Drv.add(drv); //shift u1 
     135                H.set_size ( drv._dsize() * ( -td+1 ) ); 
     136                U.set_size ( Urv._dsize() ); 
     137                for ( int i=-1;i>=td;i-- ) { 
     138                        drv.t ( -1 ); 
     139                        Drv.add ( drv ); //shift u1 
    111140                } 
    112                 rgrlnk.set_connection(rrv,Drv); 
    113                  
     141                rgrlnk.set_connection ( rrv,Drv ); 
     142 
    114143                dtsize = Drv._dsize(); 
    115144                utsize = Urv._dsize(); 
    … …  
    121150        virtual void log_add ( logger &L ) { 
    122151                //DS::log_add ( L ); too long!! 
    123                 L_dt=L.add ( Drv(0,dt_size),"" ); 
     152                L_dt=L.add ( Drv ( 0,dt_size ),"" ); 
    124153                L_ut=L.add ( Urv,"" ); 
    125154 
    … …  
    131160        virtual void logit ( logger &L ) { 
    132161                //DS::logit ( L ); 
    133                 L.logit( L_dt, H.left(dt_size)); 
    134                 L.logit(L_ut, U); 
    135                  
     162                L.logit ( L_dt, H.left ( dt_size ) ); 
     163                L.logit ( L_ut, U ); 
     164 
    136165                mat &A =model._A(); 
    137166                mat R =model._R(); 

  • bdm/stat/libEF.cpp

    r281 r283  
    183183}; 
    184184 
    185 ivec eEmp::resample ( RESAMPLING_METHOD method ) { 
     185ivec eEmp::resample (RESAMPLING_METHOD method) { 
    186186        ivec ind=zeros_i ( n ); 
    187187        ivec N_babies = zeros_i ( n ); 
    … …  
    268268} 
    269269 
    270 void eEmp::set_parameters ( const vec &w0, const epdf* epdf0 ) { 
     270void eEmp::set_statistics ( const vec &w0, const epdf* epdf0 ) { 
    271271        //it_assert_debug(rv==epdf0->rv(),"Wrong epdf0"); 
    272272        dim = epdf0->dimension(); 

  • bdm/stat/libEF.h

    r281 r283  
    9494//      virtual void flatten ( double nu0 ) {it_error ( "Not implemented" );} 
    9595 
    96         BMEF* _copy_ ( bool changerv=false ) {it_error ( "function _copy_ not implemented for this BM" ); return NULL;}; 
     96        BMEF* _copy_ ( bool changerv=false ) const {it_error ( "function _copy_ not implemented for this BM" ); return NULL;}; 
    9797}; 
    9898 
    … …  
    348348 \f] 
    349349 
     350Vector \f$\beta\f$ has different meaning (in fact it is 1/beta as used in definition of iG) 
     351 
    350352 Inverse Gamma can be converted to Gamma using \f[ 
    351353 x\sim iG(a,b) => 1/x\sim G(a,1/b) 
    … …  
    354356 */ 
    355357 
    356 class eigamma : public eEF { 
    357 protected: 
    358         //!internal egamma 
    359         egamma eg; 
    360         //! Vector \f$\alpha\f$ 
    361         vec &alpha; 
    362         //! Vector \f$\beta\f$ (in fact it is 1/beta as used in definition of iG) 
    363         vec &beta; 
     358class eigamma : public egamma { 
     359protected: 
    364360public : 
    365361        //! \name Constructors 
    366         //!@{ 
    367         eigamma ( ) :eEF ( ), eg(),alpha ( eg._alpha() ), beta ( eg._beta() ) {}; 
    368         eigamma ( const vec &a, const vec &b ) :eEF ( ), eg(),alpha ( eg._alpha() ), beta ( eg._beta() ) {eg.set_parameters ( a,b );}; 
    369         void set_parameters ( const vec &a, const vec &b ) {eg.set_parameters ( a,b );}; 
    370         //!@} 
    371  
    372         vec sample() const {return 1.0/eg.sample();}; 
    373         //! TODO: is it used anywhere? 
    374 //      mat sample ( int N ) const; 
    375         double evallog ( const vec &val ) const {return eg.evallog ( val );}; 
    376         double lognc () const {return eg.lognc();}; 
     362        //! All constructors are inherited 
     363        //!@{ 
     364        //!@} 
     365 
     366        vec sample() const {return 1.0/egamma::sample();}; 
    377367        //! Returns poiter to alpha and beta. Potentially dangerous: use with care! 
    378368        vec mean() const {return elem_div ( beta,alpha-1 );} 
    379369        vec variance() const {vec mea=mean(); return elem_div ( elem_mult ( mea,mea ),alpha-2 );} 
    380         vec& _alpha() {return alpha;} 
    381         vec& _beta() {return beta;} 
    382370}; 
    383371/* 
    … …  
    490478        mgnorm() :mu ( epdf._mu() ) {ep=&epdf;} 
    491479        //!set mean function 
    492         void set_parameters ( fnc* g0, const sq_T &R0 ) {g=g0; epdf.set_parameters ( zeros ( g->_dimy() ), R0 );} 
     480        void set_parameters ( fnc* g0, const sq_T &R0 ) {g=g0; epdf.set_parameters ( zeros ( g->dimension() ), R0 );} 
    493481        void condition ( const vec &cond ) {mu=g->eval ( cond );}; 
    494482}; 
    … …  
    686674 
    687675        //! Set samples and weights 
    688         void set_parameters ( const vec &w0, const epdf* pdf0 ); 
     676        void set_statistics ( const vec &w0, const epdf* pdf0 ); 
     677        //! Set samples and weights 
     678        void set_statistics ( const epdf* pdf0 , int n ) {set_statistics ( ones ( n ) /n,pdf0 );}; 
    689679        //! Set sample 
    690680        void set_samples ( const epdf* pdf0 ); 
    691681        //! Set sample 
    692         void set_n ( int n0, bool copy=true ) {n=n0; w.set_size ( n0,copy );samples.set_size ( n0,copy );}; 
     682        void set_parameters ( int n0, bool copy=true ) {n=n0; w.set_size ( n0,copy );samples.set_size ( n0,copy );}; 
    693683        //! Potentially dangerous, use with care. 
    694684        vec& _w()  {return w;}; 
    … …  
    700690        const Array<vec>& _samples() const {return samples;}; 
    701691        //! Function performs resampling, i.e. removal of low-weight samples and duplication of high-weight samples such that the new samples represent the same density. 
    702         ivec resample ( RESAMPLING_METHOD method = SYSTEMATIC ); 
     692        ivec resample ( RESAMPLING_METHOD method=SYSTEMATIC ); 
    703693        //! inherited operation : NOT implemneted 
    704694        vec sample() const {it_error ( "Not implemented" );return 0;} 
    … …  
    714704                for ( int i=0;i<n;i++ ) {pom+=pow ( samples ( i ),2 ) *w ( i );} 
    715705                return pom-pow ( mean(),2 ); 
     706        } 
     707        //! For this class, qbounds are minimum and maximum value of the population! 
     708        void qbounds ( vec &lb, vec &ub, double perc=0.95 ) const { 
     709                // lb in inf so than it will be pushed below; 
     710                lb.set_size(dim); 
     711                ub.set_size(dim); 
     712                lb = std::numeric_limits<double>::infinity(); 
     713                ub = -std::numeric_limits<double>::infinity(); 
     714                int j; 
     715                for ( int i=0;i<n;i++ ) { 
     716                        for ( j=0;j<dim; j++ ) { 
     717                                if ( samples ( i ) ( j ) <lb ( j ) ) {lb ( j ) =samples ( i ) ( j );} 
     718                                if ( samples ( i ) ( j ) >ub ( j ) ) {ub ( j ) =samples ( i ) ( j );} 
     719                        } 
     720                } 
    716721        } 
    717722};