Changeset 974

Timestamp:

05/23/10 11:40:48 (14 years ago)

Author:

smidl

Message:

Noise Particle

Files:

• applications/bdmtoolbox/sandbox/mpf_arx_3.m (modified) (7 diffs)
• library/bdm/estim/particles.h (modified) (12 diffs)

applications/bdmtoolbox/sandbox/mpf_arx_3.m

@@ -8 +8 @@
 g.function = 'test_function';
 
-g.class = 'linfn';
-g.A     = eye(2);
-g.B     = [1;0];
+%     g.class = 'linfn';
+%     g.A     = eye(2);
+%     g.B     = [1;0];
 
 h.class = 'linfn';
@@ -17 +17 @@
 
 fx.class = 'mgnorm<chmat>';
-fx.R = [0.3 -0.2; -0.2 0.5];
+fx.R = [0.03 -0.02; -0.02 0.05];
 fx.g = g;
 fx.rv = x;
@@ -23 +23 @@
 
 fy.class = 'mgnorm<chmat>';
-fy.R = 0.1*eye(2);
+fy.R = 0.01*eye(2);
 fy.g = h;
 fy.rv = y;
@@ -39 +39 @@
 A.class = 'ARX';
 A.yrv = RV('vw',4);
+A.rv = RV('R',16);
 A.rgr = RV({});
 A.dimx=4;
 A.constant = 0;
-A.frg=0.99;
+A.frg=1.0;
+A.prior.class='egiw';
+A.prior.dimx=4;
+A.prior.V =1e-3*eye(4);
+
 
 M.class = 'NoiseParticle';
 M.g = g;
 M.h = h;
+M.yrv = y;
 M.rvx = x;
 M.rvxc = RVtimes(x,-1);
@@ -54 +60 @@
 PF.class='PF';
 PF.particle = M;
-PF.n = 100;
-PF.res_threshold = 1.0;
+PF.n = 1000;
+PF.res_threshold = 0.9;
 PF.prior.class = 'enorm<ldmat>';
 PF.prior.mu = [0.2;0.3];
@@ -61 +67 @@
 
 
-exper.ndat = 2000;
+exper.ndat = 200;
 O = estimator(DS,{PF},exper);
 %%%%%% ARX estimator conditioned on frg
@@ -69 +75 @@
 figure(1);
 hold off
-plot(M.Est0_Q_Q);
-hold on
-plot(ones(size(M.Est0_Q_Q,1),1)* DS.pdf.pdfs{2}.R(:)','--');
+plot(O.Est0_apost_mean_R);
 
 figure(2)
 hold off
-plot(M.Est0_apost_mean_);
+plot(O.Est0_apost_mean_x);
 hold on
-plot(M.DS_x,'--');
+plot(O.DS_dt_x,'--');
 

library/bdm/estim/particles.h

@@ -43 +43 @@
                         bm = m2.bm->_copy();
                         est_emp = m2.est_emp;
+                        est.validate();;
                         validate();
                 };
@@ -66 +67 @@
                 void validate() {
                         BM::validate();
+                        //est.validate(); --pdfs not known
                         bdm_assert(bm,"Internal BM is not given");
                 }
@@ -118 +120 @@
                 }
                 void validate(){
+                        MarginalizedParticleBase::validate();
+                        est_emp.set_rv(par->_rv());
                         if (est_emp.point.length()!=par->dimension())
                                 est_emp.set_point(zeros(par->dimension()));
@@ -324 +328 @@
                 shared_ptr<BM> bm0 = UI::build<BM>(set, "particle",UI::compulsory);
                 
-                shared_ptr<epdf> pri = UI::build<epdf> ( set, "prior", UI::compulsory );
                 n =0;
                 UI::get(n,set,"n",UI::optional);;
@@ -332 +335 @@
                         w = ones(n)/n;
                 }
-                set_prior(pri.get());
                 // set resampling method
                 resmethod_from_set ( set );
@@ -435 +437 @@
 
                 if ( particles(0)->_rv()._dsize() > 0 ) {
-                        bdm_assert (  particles(0)->_rv()._dsize() == est.dimension(), "Mismatch of RV and dimension of posterior" );
+                        bdm_assert (  particles(0)->_rv()._dsize() == est.dimension(), "Mismatch of RV " +particles(0)->_rv().to_string() + 
+                        " of size (" +num2str(particles(0)->_rv()._dsize())+"and dimension of posterior ("+num2str(est.dimension()) + ")" );
                 }
         }
@@ -445 +448 @@
                 for (int i = 0; i < n; i++ ) {
                         if ( ind ( i ) != i ) {
+                                delete particles(i);
                                 particles( i ) = particles( ind ( i ) )->_copy();
                         }
@@ -453 +457 @@
         Array<BM*>& _particles() {
                 return particles;
+        }
+        ~PF(){
+                for (int i=0; i<particles.length(); i++){delete particles(i);}
         }
 
@@ -476 +483 @@
                 //! function transforming xt,ut -> yt
                 shared_ptr<fnc> h; // pdf for non-linear part
-
-                RV rvv;
                 
                 RV rvx;
@@ -493 +498 @@
                 
         public:
-                BM* _copy() const{return new NoiseParticle();};
+                BM* _copy() const{return new NoiseParticle(*this);};
                 void bayes(const vec &dt, const vec &cond){
-                        epdf* pred_vw=bm->epredictor();
-                        shared_ptr<epdf> pred_v = pred_vw->marginal(rvv);
+                        shared_ptr<epdf> pred_vw=bm->epredictor();
+                        shared_ptr<epdf> pred_v = pred_vw->marginal(rvx);
                         
                         vec vt=pred_v->sample();
@@ -506 +511 @@
                         cond2g.filldown(cond,g_args);
                         vec xt = g->eval(g_args) + vt;
+                        est_emp.point=xt;
                         
                         // residue of observation
@@ -521 +527 @@
                         UI::get(g,set,"g",UI::compulsory);
                         UI::get(h,set,"h",UI::compulsory);
-                        RV rvx;
                         UI::get(rvx,set,"rvx",UI::compulsory);
                         est_emp.set_rv(rvx);
                         
-                        RV rvxc;
                         UI::get(rvxc,set,"rvxc",UI::compulsory);
-                        RV rvyc;
                         UI::get(rvyc,set,"rvyc",UI::compulsory);
                 }
                 void validate(){
                         MarginalizedParticleBase::validate();
+                        
+                        dimy = h->dimension();
+                        bm->set_yrv(concat(rvx,yrv));
+                        
+                        est_emp.set_rv(rvx);
+                        est_emp.set_dim(rvx._dsize());
+                        est.validate();
+                        // 
                         //check dimensions
-                        rvc = rvxc;
+                        rvc = rvxc.subt(rvx.copy_t(-1));
                         rvc.add( rvyc);
-                        rvc.subt(rvx);
-                
+                        rvc=rvc.subt(rvx);
+                        
+                        bdm_assert(g->dimension()==rvx._dsize(),"rvx is not described");
+                        bdm_assert(g->dimensionc()==rvxc._dsize(),"rvxc is not described");
+                        bdm_assert(h->dimension()==rvyc._dsize(),"rvyc is not described");
+                        
+                        bdm_assert(bm->dimensiony()==g->dimension()+h->dimension(), 
+                                           "Incompatible noise estimator of dimension " +
+                                                num2str(bm->dimensiony()) + " does not match dimension of g and h, " + 
+                                                num2str(g->dimension())+" and "+ num2str(h->dimension()) );
+                                                
+                        dimc = rvc._dsize();
+                                        
                         //establish datalinks
                         x2g.set_connection(rvxc, rvx.copy_t(-1));