Changeset 723 for library/bdm/estim/kalman.h

Timestamp:

11/15/09 23:02:02 (15 years ago)

Author:

smidl

Message:

Big commit of LQG stuff

Files:

• library/bdm/estim/kalman.h (modified) (7 diffs)

library/bdm/estim/kalman.h

@@ -27 +27 @@
  * \brief Basic elements of linear state-space model
 
-Parameter evolution model:\f[ x_t = A x_{t-1} + B u_t + Q^{1/2} e_t \f]
-Observation model: \f[ y_t = C x_{t-1} + C u_t + Q^{1/2} w_t. \f]
-Where $e_t$ and $w_t$ are independent vectors Normal(0,1)-distributed disturbances.
+Parameter evolution model:\f[ x_{t+1} = A x_{t} + B u_t + Q^{1/2} e_t \f]
+Observation model: \f[ y_t = C x_{t} + C u_t + R^{1/2} w_t. \f]
+Where $e_t$ and $w_t$ are mutually independent vectors of Normal(0,1)-distributed disturbances.
  */
 template<class sq_T>
@@ -61 +61 @@
                         UI::get (C, set, "C", UI::compulsory);
                         UI::get (D, set, "D", UI::compulsory);
-                        mat Qtm, Rtm;
+                        mat Qtm, Rtm; // full matrices
                         if(!UI::get(Qtm, set, "Q", UI::optional)){
                                 vec dq;
@@ -72 +72 @@
                                 Rtm=diag(dr);
                         }
-                        R=Rtm; // automatic conversion
+                        R=Rtm; // automatic conversion to square-root form
                         Q=Qtm; 
                         
@@ -519 +519 @@
 
 */
-class StateFromARX: public StateSpace<fsqmat>{
+class StateFromARX: public StateSpace<chmat>{
         protected:
                 //! remember connection from theta ->A
                 datalink_part th2A;
-                //! xrv
-                RV xrv;
+                //! remember connection from theta ->B
+                datalink_part th2B;
                 //!function adds n diagonal elements from given starting point r,c
                 void diagonal_part(mat &A, int r, int c, int n){
                         for(int i=0; i<n; i++){A(r,c)=1.0; r++;c++;}
                 };
+                //! similar to ARX.have_constant
+                bool have_constant;
         public:
                 //! set up this object to match given mlnorm
-                void connect_mlnorm(const mlnorm<fsqmat> &ml){
+                //! Note that state-space and common mpdf use different meaning of \f$ _t \f$ in \f$ u_t \f$. 
+                //!While mlnorm typically assumes that \f$ u_t \rightarrow y_t \f$ in state space it is \f$ u_{t-1} \rightarrow y_t \f$
+                //! For consequences in notation of internal variable xt see arx2statespace_notes.lyx.
+                void connect_mlnorm(const mlnorm<chmat> &ml, RV &xrv, RV &urv){
+
                         //get ids of yrv                                
                         const RV &yrv = ml._rv();
@@ -537 +543 @@
                         const RV &rgr = ml._rvc();
                         RV rgr0 = rgr.remove_time();
-                        RV urv = rgr0.subt(yrv); 
+                        urv = rgr0.subt(yrv); 
                                                 
                         // create names for state variables
                         xrv=yrv; 
                         
-                        int y_multiplicity = rgr.mint(yrv);
-                        int y_block_size = yrv.length()*(y_multiplicity+1); // current yt + all delayed yts
-                        for (int m=0;m<y_multiplicity;m++){ 
+                        int y_multiplicity = -rgr.mint(yrv);
+                        int y_block_size = yrv.length()*(y_multiplicity); // current yt + all delayed yts
+                        for (int m=0;m<y_multiplicity - 1;m++){ // ========= -1 is important see arx2statespace_notes
                                 xrv.add(yrv.copy_t(-m-1)); //add delayed yt
                         }
+                        //! temporary RV for connection to ml.rvc, since notation of xrv and ml.rvc does not match
+                        RV xrv_ml = xrv.copy_t(-1);
+                        
                         // add regressors
                         ivec u_block_sizes(urv.length()); // no_blocks = yt + unique rgr
@@ -552 +561 @@
                                 RV R=urv.subselect(vec_1(r)); //non-delayed element of rgr
                                 int r_size=urv.size(r);
-                                int r_multiplicity=rgr.mint(R);
+                                int r_multiplicity=-rgr.mint(R);
                                 u_block_sizes(r)= r_size*r_multiplicity ;
                                 for(int m=0;m<r_multiplicity;m++){ 
                                         xrv.add(R.copy_t(-m-1)); //add delayed yt
+                                        xrv_ml.add(R.copy_t(-m-1)); //add delayed yt
                                 }
                         }
-                        
+                        // add constant
+                        if (any(ml._mu_const()!=0.0)){
+                                have_constant=true;
+                                xrv.add(RV("bdm_reserved_constant_one",1));
+                        } else {have_constant=false;}
+                                
+                
                         // get mapp
-                        th2A.set_connection(xrv, ml._rvc());
+                        th2A.set_connection(xrv_ml, ml._rvc());
+                        th2B.set_connection(urv, ml._rvc());
                         
                         //set matrix sizes
@@ -577 +594 @@
                         this->C=zeros(yrv._dsize(), xrv._dsize());
                         this->C.set_submatrix(0,0,eye(yrv._dsize()));
-                                this->D=zeros(yrv._dsize(),urv._dsize());
-                                this->R = zeros(yrv._dsize(),yrv._dsize());
-                                this->Q = zeros(xrv._dsize(),xrv._dsize());
-                                // Q is set by update
-                                                                
-                                update_from(ml);
-                                validate();
+                        this->D=zeros(yrv._dsize(),urv._dsize());
+                        this->R.setCh(zeros(yrv._dsize(),yrv._dsize()));
+                        this->Q.setCh(zeros(xrv._dsize(),xrv._dsize()));
+                        // Q is set by update
+                                                        
+                        update_from(ml);
+                        validate();
                 };
                 //! fast function to update parameters from ml - not checked for compatibility!!
-                void update_from(const mlnorm<fsqmat> &ml){
-                        
-                        vec theta = ml._A().get_row(0); // this 
-                        this->Q._M().set_submatrix(0,0,ml._R());
-                        
-                }
+                void update_from(const mlnorm<chmat> &ml){
+
+                        vec Arow=zeros(A.cols());
+                        vec Brow=zeros(B.cols());
+                        //  ROW- WISE EVALUATION =====
+                        for(int i=0;i<ml._rv()._dsize(); i++){ 
+                                
+                                vec theta = ml._A().get_row(i); 
+                                
+                                th2A.filldown(theta,Arow);
+                                if (have_constant){
+                                        // constant is always at the end no need for datalink
+                                        Arow(A.cols()-1) = ml._mu_const()(i);
+                                }
+                                this->A.set_row(i,Arow);
+                                
+                                th2B.filldown(theta,Brow);
+                                this->B.set_row(i,Brow);
+                        }
+                        this->Q._Ch().set_submatrix(0,0,ml.__R()._Ch());
+                        
+                };
+                //! access function
+                bool _have_constant() const {return have_constant;}
 };