root/library/bdm/estim/particles.cpp @ 1226

Revision 1196, 6.1 kB (checked in by smidl, 14 years ago)

monitor Neff + test enorm

  • Property svn:eol-style set to native
Line 
1#include "particles.h"
2
3namespace bdm {
4
5using std::endl;
6
7
8void PF::bayes_weights() {
9//
10    double mlls = max ( lls );
11    // compute weights
12    for ( int  i = 0; i < n; i++ ) {
13        w ( i ) *= exp ( lls ( i ) - mlls ); // multiply w by likelihood
14    }
15
16    //renormalize
17    double sw = sum ( w );
18    if ( !std::isfinite ( sw ) ) {
19        for ( int i = 0; i < n; i++ ) {
20            if ( !std::isfinite ( w ( i ) ) ) {
21                w ( i ) = 0;
22            }
23        }
24        sw = sum ( w );
25        if ( !std::isfinite ( sw ) || sw == 0.0 ) {
26            bdm_error ( "Particle filter is lost; no particle is good enough." );
27        }
28    }
29    w /= sw;
30}
31
32void PF::bayes ( const vec &yt, const vec &cond ) {
33
34    /*  if (auxiliary){
35                ...
36        }*/
37    // weight them - data step
38    for (int i = 0; i < n; i++ ) {
39        particles(i)->bayes(yt,cond);
40        lls ( i ) = particles(i)->_ll(); //+= because lls may have something from gensmp!
41    }
42
43    bayes_weights();
44
45        if (log_level[logweights]) {
46                log_level.store( logweights, w);
47        }
48        if (log_level[logNeff]) {
49                log_level.store( logNeff, 1.0 / ( w * w ));
50        }
51       
52    if ( do_resampling() ) {
53        resample ( );
54    }
55
56}
57
58// void PF::set_est ( const epdf &epdf0 ) {
59//      int i;
60//
61//      for ( i=0;i<n;i++ ) {
62//              _samples ( i ) = epdf0.sample();
63//      }
64// }
65
66// vec MPF::mpfepdf::mean() const {
67//      const vec &w = pf->posterior()._w();
68//      vec pom = zeros ( BMs ( 0 )->posterior ().dimension() );
69//      //compute mean of BMs
70//      for ( int i = 0; i < w.length(); i++ ) {
71//              pom += BMs ( i )->posterior().mean() * w ( i );
72//      }
73//      return concat ( pf->posterior().mean(), pom );
74// }
75//
76// vec MPF::mpfepdf::variance() const {
77//      const vec &w = pf->posterior()._w();
78//
79//      vec pom = zeros ( BMs ( 0 )->posterior ().dimension() );
80//      vec pom2 = zeros ( BMs ( 0 )->posterior ().dimension() );
81//      vec mea;
82//
83//      for ( int i = 0; i < w.length(); i++ ) {
84//              // save current mean
85//              mea = BMs ( i )->posterior().mean();
86//              pom += mea * w ( i );
87//              //compute variance
88//              pom2 += ( BMs ( i )->posterior().variance() + pow ( mea, 2 ) ) * w ( i );
89//      }
90//      return concat ( pf->posterior().variance(), pom2 - pow ( pom, 2 ) );
91// }
92//
93// void MPF::mpfepdf::qbounds ( vec &lb, vec &ub, double perc ) const {
94//      //bounds on particles
95//      vec lbp;
96//      vec ubp;
97//      pf->posterior().qbounds ( lbp, ubp );
98//
99//      //bounds on Components
100//      int dimC = BMs ( 0 )->posterior().dimension();
101//      int j;
102//      // temporary
103//      vec lbc ( dimC );
104//      vec ubc ( dimC );
105//      // minima and maxima
106//      vec Lbc ( dimC );
107//      vec Ubc ( dimC );
108//      Lbc = std::numeric_limits<double>::infinity();
109//      Ubc = -std::numeric_limits<double>::infinity();
110//
111//      for ( int i = 0; i < BMs.length(); i++ ) {
112//              // check Coms
113//              BMs ( i )->posterior().qbounds ( lbc, ubc );
114//              //save either minima or maxima
115//              for ( j = 0; j < dimC; j++ ) {
116//                      if ( lbc ( j ) < Lbc ( j ) ) {
117//                              Lbc ( j ) = lbc ( j );
118//                      }
119//                      if ( ubc ( j ) > Ubc ( j ) ) {
120//                              Ubc ( j ) = ubc ( j );
121//                      }
122//              }
123//      }
124//      lb = concat ( lbp, Lbc );
125//      ub = concat ( ubp, Ubc );
126// }
127//
128//
129//
130// void MPF::bayes ( const vec &yt, const vec &cond ) {
131//      // follows PF::bayes in most places!!!
132//      int i;
133//      int n = pf->__w().length();
134//      vec &lls = pf->_lls();
135//      Array<vec> &samples = pf->__samples();
136//
137//      // generate samples - time step
138//      pf->bayes_gensmp ( vec ( 0 ) );
139//      // weight them - data step
140// #pragma parallel for
141//      for ( i = 0; i < n; i++ ) {
142//              vec bm_cond ( BMs ( i )->dimensionc() );
143//              this2bm.fill_cond ( yt, cond, bm_cond );
144//              pf2bm.filldown ( samples ( i ), bm_cond );
145//              BMs ( i ) -> bayes ( this2bm.pushdown ( yt ), bm_cond );
146//              lls ( i ) += BMs ( i )->_ll();
147//      }
148//
149//      pf->bayes_weights();
150//
151//      ivec ind;
152//      if ( pf->do_resampling() ) {
153//              pf->resample ( ind );
154//
155// #pragma omp parallel for
156//              for ( i = 0; i < n; i++ ) {
157//                      if ( ind ( i ) != i ) {//replace the current Bm by a new one
158//                              delete BMs ( i );
159//                              BMs ( i ) = (BM*) BMs ( ind ( i ) )->_copy(); //copy constructor
160//                      }
161//              };
162//      }
163// };
164//
165//
166// void MPF_ARXg::bayes ( const vec &yt, const vec &cond ) {
167//      // follows PF::bayes in most places!!!
168//      int i;
169//      int n = pf->__w().length();
170//      vec &lls = pf->_lls();
171//      Array<vec> &samples = pf->__samples();
172//
173//      // generate samples - time step
174//      for (int i =0;i<n; i++){
175//              mat M; chmat R;
176//              BMsp(i)->posterior().sample_mat(M,R);
177//              vec diff=R._Ch().T()*randn(samples(i).length());
178//              samples(i) = g->eval(samples(i)) + diff;
179//              ////////////
180// //           samples(i) = cond;
181//              /////////
182//              BMsp(i)->bayes(diff);
183//      }
184//      // weight them - data step
185//      enorm<ldmat> ooo;
186//      ooo.set_parameters(zeros(2),0.1*eye(2));
187//      ooo.validate();
188//
189//      #pragma parallel for
190//      for ( i = 0; i < n; i++ ) {
191//              vec tmp=yt - h->eval(samples(i));
192//              BMso ( i ) -> bayes ( tmp );
193//              lls ( i ) += BMso ( i )->_ll();
194//              /////
195//              ooo._mu()=h->eval(samples(i));
196//              lls(i) = ooo.evallog_nn(yt);
197//      }
198//
199//      pf->bayes_weights();
200//
201//      ivec ind;
202//      if ( pf->do_resampling() ) {
203//              pf->resample ( ind );
204//
205//              #pragma omp parallel for
206//              for ( i = 0; i < n; i++ ) {
207//                      if ( ind ( i ) != i ) {//replace the current Bm by a new one
208//                              delete BMso ( i );
209//                              BMso ( i ) = (ARX*) BMso ( ind ( i ) )->_copy(); //copy constructor
210//                              delete BMsp ( i );
211//                              BMsp ( i ) = (ARX*) BMsp ( ind ( i ) )->_copy(); //copy constructor
212//                      }
213//              };
214//      }
215// };
216//
217// void MPF_ARXg::from_setting(const libconfig::Setting& set) {
218//      BM::from_setting(set);
219//
220//      pf = new PF;
221//      // prior must be set before BM
222//      pf->prior_from_set ( set );
223//      pf->resmethod_from_set ( set );
224//
225//      // read functions g and h
226//      g=UI::build<fnc>(set,"g");
227//      h=UI::build<fnc>(set,"h");
228//
229//      set_dim( g->dimension());
230//      dimy = h->dimension();
231//
232//      shared_ptr<ARX> arxo=UI::build<ARX>(set,"arxo");
233//      shared_ptr<ARX> arxp=UI::build<ARX>(set,"arxp");
234//      int n = pf->__samples().length();
235//      BMso.set_length(n);
236//      BMsp.set_length(n);
237//      for(int i=0; i<n; i++){
238//              BMso(i)=arxo->_copy();
239//              BMsp(i)=arxp->_copy();
240//      }
241//
242//      yrv = arxo->_yrv();
243//      //rvc = arxo->_rvc();
244//
245//      validate();
246//
247// }
248
249
250}
251//MPF::MPF:{}
Note: See TracBrowser for help on using the browser.