Changeset 811 for library/bdm/estim

Timestamp:

02/21/10 20:58:51 (14 years ago)

Author:

smidl

Message:

extensions and stuff for MPF

Location:

library/bdm/estim

Files:

• particles.cpp (modified) (1 diff)
• particles.h (modified) (3 diffs)

library/bdm/estim/particles.cpp

@@ -175 +175 @@
 
 
+void MPF_ARXg::bayes ( const vec &yt, const vec &cond ) {
+        // follows PF::bayes in most places!!!
+        int i;
+        int n = pf->__w().length();
+        vec &lls = pf->_lls();
+        Array<vec> &samples = pf->__samples();
+        
+        // generate samples - time step
+        for (int i =0;i<n; i++){
+                mat M; chmat R;
+                BMsp(i)->posterior().sample_mat(M,R);
+                vec diff=R._Ch().T()*randn(samples(i).length());
+                samples(i) = g->eval(samples(i)) + diff;
+                ////////////
+//              samples(i) = cond;
+                /////////
+                BMsp(i)->bayes(diff);
+        }
+        // weight them - data step
+        enorm<ldmat> ooo;
+        ooo.set_parameters(zeros(2),0.1*eye(2));
+        ooo.validate();
+        
+        #pragma parallel for
+        for ( i = 0; i < n; i++ ) {
+                vec tmp=yt - h->eval(samples(i));
+                BMso ( i ) -> bayes ( tmp );
+                lls ( i ) += BMso ( i )->_ll();
+                /////
+                ooo._mu()=h->eval(samples(i));
+                lls(i) = ooo.evallog_nn(yt);
+        }
+        
+        pf->bayes_weights();
+        
+        ivec ind;
+        if ( pf->do_resampling() ) {
+                pf->resample ( ind );
+                
+                #pragma omp parallel for
+                for ( i = 0; i < n; i++ ) {
+                        if ( ind ( i ) != i ) {//replace the current Bm by a new one
+                                delete BMso ( i );
+                                BMso ( i ) = (ARX*) BMso ( ind ( i ) )->_copy(); //copy constructor
+                                delete BMsp ( i );
+                                BMsp ( i ) = (ARX*) BMsp ( ind ( i ) )->_copy(); //copy constructor
+                        }
+                };
+        }
+};
+
+void MPF_ARXg::from_setting(const libconfig::Setting& set) {
+        BM::from_setting(set);
+        
+        pf = new PF;
+        // prior must be set before BM
+        pf->prior_from_set ( set );
+        pf->resmethod_from_set ( set );
+        
+        string opt;
+        if ( UI::get ( opt, set, "options", UI::optional ) ) {
+                set_options ( opt );
+        }
+        
+        // read functions g and h
+        g=UI::build<fnc>(set,"g");
+        h=UI::build<fnc>(set,"h");
+        
+        set_dim( g->dimension());
+        dimy = h->dimension();
+        
+        shared_ptr<ARX> arxo=UI::build<ARX>(set,"arxo");
+        shared_ptr<ARX> arxp=UI::build<ARX>(set,"arxp");
+        int n = pf->__samples().length();
+        BMso.set_length(n);
+        BMsp.set_length(n);
+        for(int i=0; i<n; i++){
+                BMso(i)=arxo->_copy();
+                BMsp(i)=arxp->_copy();
+        }
+        
+        yrv = arxo->_yrv();
+        //rvc = arxo->_rvc();
+        
+        validate();
+        
+}
+
+
 }
 //MPF::MPF:{}

library/bdm/estim/particles.h

@@ -15 +15 @@
 
 
-#include "../stat/exp_family.h"
+#include "../estim/arx_ext.h"
 
 namespace bdm {
@@ -183 +183 @@
                 }
                 n = est._w().length();
-                //validate();
+                lls=zeros(n);
         }
 
@@ -363 +363 @@
 UIREGISTER ( MPF );
 
+/*! ARXg for estimation of state-space variances
+*/
+class MPF_ARXg :public BM{
+        protected:
+        shared_ptr<PF> pf;
+        //! pointer to Array of BMs
+        Array<ARX*> BMso;
+        //! pointer to Array of BMs
+        Array<ARX*> BMsp;
+        //!parameter evolution
+        shared_ptr<fnc> g;
+        //!observation function
+        shared_ptr<fnc> h;
+        
+        public:
+                void bayes(const vec &yt, const vec &cond );
+                void from_setting(const Setting &set) ;
+                void validate() {
+                        bdm_assert(g->dimensionc()==g->dimension(),"not supported yet");
+                        bdm_assert(h->dimensionc()==g->dimension(),"not supported yet");                        
+                }
+
+                double logpred(const vec &cond) const NOT_IMPLEMENTED(0.0);
+                epdf* epredictor() const NOT_IMPLEMENTED(NULL);
+                pdf* predictor() const NOT_IMPLEMENTED(NULL);
+                const epdf& posterior() const {return pf->posterior();};
+                
+                void log_register( logger &L, const string &prefix ){
+                        BM::log_register(L,prefix);
+                        logrec->ids.set_size ( 3 );
+                        logrec->ids(1)= L.add_vector(RV("Q",dimension()*dimension()), prefix+L.prefix_sep()+"Q");
+                        logrec->ids(2)= L.add_vector(RV("R",dimensiony()*dimensiony()), prefix+L.prefix_sep()+"R");
+                        
+                };
+                void log_write() const {
+                        BM::log_write();
+                        mat mQ=zeros(dimension(),dimension());
+                        mat pom=zeros(dimension(),dimension());
+                        mat mR=zeros(dimensiony(),dimensiony());
+                        mat pom2=zeros(dimensiony(),dimensiony());
+                        mat dum;
+                        const vec w=pf->posterior()._w();
+                        for (int i=0; i<w.length(); i++){
+                                BMsp(i)->posterior().mean_mat(dum,pom);
+                                mQ += w(i) * pom;
+                                BMso(i)->posterior().mean_mat(dum,pom2);
+                                mR += w(i) * pom2;
+                                
+                        }
+                        logrec->L.log_vector ( logrec->ids ( 1 ), cvectorize(mQ) );
+                        logrec->L.log_vector ( logrec->ids ( 2 ), cvectorize(mR) );
+                        
+                }
+};
+UIREGISTER(MPF_ARXg);
+
+
 }
 #endif // KF_H