Changeset 679 for library/bdm/estim/arx.cpp

Timestamp:

10/23/09 00:05:25 (15 years ago)

Author:

smidl

Message:

Major changes in BM -- OK is only test suite and tests/tutorial -- the rest is broken!!!

Files:

• library/bdm/estim/arx.cpp (modified) (16 diffs)

library/bdm/estim/arx.cpp

@@ -2 +2 @@
 namespace bdm {
 
-void ARX::bayes ( const vec &dt, const double w ) {
+        void ARX::bayes_weighted ( const vec &yt, const vec &cond, const double w ) {
         double lnc;
-                
+        //cache
+        ldmat &V=est._V(); 
+        double &nu=est._nu();
+
+        dyad.set_subvector(0,yt);
+        dyad.set_subvector(dimy,cond);
+        // possible "1" is there from the beginning
+        
         if ( frg < 1.0 ) {
                 est.pow ( frg ); // multiply V and nu
+
                 
                 //stabilize
                 ldmat V0=alter_est._V(); //$ copy
                 double &nu0=alter_est._nu();
+                
                 V0*=(1-frg);
                 V += V0; //stabilization
@@ -20 +29 @@
                 }
         }
-        if (have_constant) {
-                _dt.set_subvector(0,dt);
-                V.opupdt ( _dt, w );
-        } else {
-                V.opupdt ( dt, w );
-        }
+        V.opupdt ( dyad, w );
         nu += w;
 
@@ -36 +40 @@
 }
 
-double ARX::logpred ( const vec &dt ) const {
+double ARX::logpred ( const vec &yt ) const {
         egiw pred ( est );
         ldmat &V = pred._V();
@@ -42 +46 @@
 
         double lll;
+        vec dyad_p = dyad;
+        dyad_p.set_subvector(0,yt);
 
         if ( frg < 1.0 ) {
@@ -53 +59 @@
                 }
 
-        V.opupdt ( dt, 1.0 );
+        V.opupdt ( dyad_p, 1.0 );
         nu += 1.0;
         // log(sqrt(2*pi)) = 0.91893853320467
@@ -67 +73 @@
         const ARX* A0 = dynamic_cast<const ARX*> ( B0 );
 
-        bdm_assert_debug ( V.rows() == A0->V.rows(), "ARX::set_statistics Statistics  differ" );
-        set_statistics ( A0->dimx, A0->V, A0->nu );
+        bdm_assert_debug ( dimension() == A0->dimension(), "Statistics of different dimensions" );
+        set_statistics ( A0->dimensiony(), A0->posterior()._V(), A0->posterior()._nu() );
 }
 
 enorm<ldmat>* ARX::epredictor ( const vec &rgr ) const {
-        int dim = dimx;//est.dimension();
-        mat mu ( dim, V.rows() - dim );
-        mat R ( dim, dim );
+        mat mu ( dimy, posterior()._V().rows() - dimy );
+        mat R ( dimy, dimy );
 
         enorm<ldmat>* tmp;
         tmp = new enorm<ldmat> ( );
         //TODO: too hackish
-        if ( drv._dsize() > 0 ) {
+        if ( yrv._dsize() > 0 ) {
         }
 
@@ -91 +96 @@
 
 mlnorm<ldmat>* ARX::predictor ( ) const {
-        int dim = est.dimension();
-        
-        mat mu ( dim, V.rows() - dim );
-        mat R ( dim, dim );
+        mat mu ( dimy, posterior()._V().rows() - dimy );
+        mat R ( dimy, dimy );
         mlnorm<ldmat>* tmp;
         tmp = new mlnorm<ldmat> ( );
@@ -106 +109 @@
                 tmp->set_parameters ( mu.get_cols ( 0, mu.cols() - 2 ), mu.get_col ( mu.cols() - 1 ), ldmat ( R ) );
         } else {
-                tmp->set_parameters ( mu, zeros ( dim ), ldmat ( R ) );
+                tmp->set_parameters ( mu, zeros ( dimy ), ldmat ( R ) );
         }
         return tmp;
@@ -112 +115 @@
 
 mlstudent* ARX::predictor_student ( ) const {
-        int dim = est.dimension();
-
-        mat mu ( dim, V.rows() - dim );
-        mat R ( dim, dim );
+        const ldmat &V = posterior()._V();
+        
+        mat mu ( dimy, V.rows() - dimy );
+        mat R ( dimy, dimy );
         mlstudent* tmp;
         tmp = new mlstudent ( );
@@ -122 +125 @@
         mu = mu.T();
 
-        int xdim = dimx;
         int end = V._L().rows() - 1;
-        ldmat Lam ( V._L() ( xdim, end, xdim, end ), V._D() ( xdim, end ) );  //exp val of R
+        ldmat Lam ( V._L() ( dimy, end, dimy, end ), V._D() ( dimy, end ) );  //exp val of R
 
 
@@ -132 +134 @@
                         tmp->set_parameters ( mu.get_cols ( 0, mu.cols() - 2 ), mu.get_col ( mu.cols() - 1 ), ldmat ( R ), Lam );
                 } else {
-                        tmp->set_parameters ( mat ( dim, 0 ), mu.get_col ( mu.cols() - 1 ), ldmat ( R ), Lam );
+                        tmp->set_parameters ( mat ( dimy, dimc ), mu.get_col ( mu.cols() - 1 ), ldmat ( R ), Lam );
                 }
         } else {
                 // no constant term
-                tmp->set_parameters ( mu, zeros ( xdim ), ldmat ( R ), Lam );
+                tmp->set_parameters ( mu, zeros ( dimy ), ldmat ( R ), Lam );
         }
         return tmp;
@@ -214 +216 @@
 
 void ARX::from_setting ( const Setting &set ) {
-        shared_ptr<RV> yrv = UI::build<RV> ( set, "rv", UI::compulsory );
+        shared_ptr<RV> yrv_ = UI::build<RV> ( set, "rv", UI::compulsory );
         shared_ptr<RV> rrv = UI::build<RV> ( set, "rgr", UI::compulsory );
-        int ylen = yrv->_dsize();
+        dimy = yrv_->_dsize();
         // rgrlen - including constant!!!
-        int rgrlen = rrv->_dsize();
-        
-        dimx = ylen;
-        set_rv ( *yrv, *rrv );
+        dimc = rrv->_dsize();
+        
+        yrv = *yrv_;
+        rvc = *rrv;
         
         string opt;
@@ -233 +235 @@
                 have_constant=constant>0;
         }
-        if (have_constant) {rgrlen++;_dt=ones(rgrlen+ylen);}
+        int rgrlen = dimc+int(have_constant==true);
 
         //init
         shared_ptr<egiw> pri=UI::build<egiw>(set, "prior", UI::optional);
         if (pri) {
-                bdm_assert(pri->_dimx()==ylen,"prior is not compatible");
-                bdm_assert(pri->_V().rows()==ylen+rgrlen,"prior is not compatible");
+                bdm_assert(pri->_dimx()==dimy,"prior is not compatible");
+                bdm_assert(pri->_V().rows()==dimy+rgrlen,"prior is not compatible");
                 est.set_parameters( pri->_dimx(),pri->_V(), pri->_nu());
         }else{
-                est.set_parameters( ylen, zeros(ylen+rgrlen));
+                est.set_parameters( dimy, zeros(dimy+rgrlen));
                 set_prior_default(est);
         }
@@ -248 +250 @@
         shared_ptr<egiw> alt=UI::build<egiw>(set, "alternative", UI::optional);
         if (alt) {
-                bdm_assert(alt->_dimx()==ylen,"alternative is not compatible");
-                bdm_assert(alt->_V().rows()==ylen+rgrlen,"alternative is not compatible");
+                bdm_assert(alt->_dimx()==dimy,"alternative is not compatible");
+                bdm_assert(alt->_V().rows()==dimy+rgrlen,"alternative is not compatible");
                 alter_est.set_parameters( alt->_dimx(),alt->_V(), alt->_nu());
         } else {
@@ -264 +266 @@
         shared_ptr<RV> rv_par=UI::build<RV>(set, "rv_param",UI::optional );
         if (!rv_par){
-                est.set_rv ( RV ( "{theta r }", vec_2 ( ylen*rgrlen, ylen*ylen ) ) );
+                est.set_rv ( RV ( "{theta r }", vec_2 ( dimy*rgrlen, dimy*dimy ) ) );
         } else {
                 est.set_rv ( *rv_par );
@@ -270 +272 @@
         validate();
 }
-
-}
+}
+