[8] | 1 | /*! |
---|
| 2 | \file |
---|
| 3 | \brief Bayesian Filtering using stochastic sampling (Particle Filters) |
---|
| 4 | \author Vaclav Smidl. |
---|
| 5 | |
---|
| 6 | ----------------------------------- |
---|
| 7 | BDM++ - C++ library for Bayesian Decision Making under Uncertainty |
---|
| 8 | |
---|
| 9 | Using IT++ for numerical operations |
---|
| 10 | ----------------------------------- |
---|
| 11 | */ |
---|
| 12 | |
---|
[384] | 13 | #ifndef PARTICLES_H |
---|
| 14 | #define PARTICLES_H |
---|
[8] | 15 | |
---|
[262] | 16 | |
---|
[811] | 17 | #include "../estim/arx_ext.h" |
---|
[887] | 18 | #include "../stat/emix.h" |
---|
[8] | 19 | |
---|
[270] | 20 | namespace bdm { |
---|
[1064] | 21 | |
---|
[1085] | 22 | //! \brief Abstract class for Marginalized Particles |
---|
[1064] | 23 | class MarginalizedParticleBase : public BM { |
---|
| 24 | protected: |
---|
| 25 | //! discrte particle |
---|
| 26 | dirac est_emp; |
---|
| 27 | //! internal Bayes Model |
---|
| 28 | shared_ptr<BM> bm; |
---|
| 29 | |
---|
[1068] | 30 | //! \brief Internal class for custom posterior - product of empirical and exact part |
---|
[1064] | 31 | class eprod_2:public eprod_base { |
---|
| 32 | protected: |
---|
| 33 | MarginalizedParticleBase ∓ |
---|
| 34 | public: |
---|
[1192] | 35 | eprod_2(MarginalizedParticleBase &m):mp(m) {} |
---|
| 36 | const epdf* factor(int i) const { |
---|
[1064] | 37 | return (i==0) ? &mp.bm->posterior() : &mp.est_emp; |
---|
| 38 | } |
---|
| 39 | const int no_factors() const { |
---|
| 40 | return 2; |
---|
| 41 | } |
---|
| 42 | } est; |
---|
| 43 | |
---|
| 44 | public: |
---|
| 45 | MarginalizedParticleBase():est(*this) {}; |
---|
| 46 | MarginalizedParticleBase(const MarginalizedParticleBase &m2):BM(m2),est(*this) { |
---|
| 47 | bm = m2.bm->_copy(); |
---|
| 48 | est_emp = m2.est_emp; |
---|
| 49 | est.validate(); |
---|
| 50 | validate(); |
---|
| 51 | }; |
---|
| 52 | void bayes(const vec &dt, const vec &cond) NOT_IMPLEMENTED_VOID; |
---|
| 53 | |
---|
| 54 | const eprod_2& posterior() const { |
---|
| 55 | return est; |
---|
| 56 | } |
---|
| 57 | |
---|
| 58 | void set_prior(const epdf *pdf0) { |
---|
| 59 | const eprod *ep=dynamic_cast<const eprod*>(pdf0); |
---|
| 60 | if (ep) { // full prior |
---|
| 61 | bdm_assert(ep->no_factors()==2,"Incompatible prod"); |
---|
| 62 | bm->set_prior(ep->factor(0)); |
---|
| 63 | est_emp.set_point(ep->factor(1)->sample()); |
---|
| 64 | } else { |
---|
| 65 | // assume prior is only for emp; |
---|
| 66 | est_emp.set_point(pdf0->sample()); |
---|
| 67 | } |
---|
| 68 | } |
---|
[1077] | 69 | |
---|
| 70 | |
---|
| 71 | /*! Create object from the following structure |
---|
| 72 | |
---|
| 73 | \code |
---|
| 74 | class = "MarginalizedParticleBase"; |
---|
| 75 | bm = configuration of bdm::BM; % any offspring of BM, bdm::BM::from_setting |
---|
| 76 | --- inherited fields --- |
---|
| 77 | bdm::BM::from_setting |
---|
| 78 | \endcode |
---|
| 79 | */ |
---|
[1064] | 80 | void from_setting(const Setting &set) { |
---|
| 81 | BM::from_setting ( set ); |
---|
| 82 | bm = UI::build<BM> ( set, "bm", UI::compulsory ); |
---|
| 83 | } |
---|
| 84 | void validate() { |
---|
| 85 | BM::validate(); |
---|
| 86 | //est.validate(); --pdfs not known |
---|
| 87 | bdm_assert(bm,"Internal BM is not given"); |
---|
| 88 | } |
---|
[971] | 89 | }; |
---|
| 90 | |
---|
[1085] | 91 | //! \brief Particle with marginalized subspace, used in PF |
---|
[1064] | 92 | class MarginalizedParticle : public MarginalizedParticleBase { |
---|
| 93 | protected: |
---|
| 94 | //! pdf with for transitional par |
---|
| 95 | shared_ptr<pdf> par; // pdf for non-linear part |
---|
| 96 | //! link from this to bm |
---|
| 97 | shared_ptr<datalink_part> cond2bm; |
---|
| 98 | //! link from cond to par |
---|
| 99 | shared_ptr<datalink_part> cond2par; |
---|
| 100 | //! link from emp 2 par |
---|
| 101 | shared_ptr<datalink_part> emp2bm; |
---|
| 102 | //! link from emp 2 par |
---|
| 103 | shared_ptr<datalink_part> emp2par; |
---|
| 104 | |
---|
| 105 | public: |
---|
| 106 | BM* _copy() const { |
---|
| 107 | return new MarginalizedParticle(*this); |
---|
| 108 | }; |
---|
| 109 | void bayes(const vec &dt, const vec &cond) { |
---|
| 110 | vec par_cond(par->dimensionc()); |
---|
| 111 | cond2par->filldown(cond,par_cond); // copy ut |
---|
| 112 | emp2par->filldown(est_emp._point(),par_cond); // copy xt-1 |
---|
| 113 | |
---|
| 114 | //sample new particle |
---|
| 115 | est_emp.set_point(par->samplecond(par_cond)); |
---|
| 116 | //if (evalll) |
---|
| 117 | vec bm_cond(bm->dimensionc()); |
---|
| 118 | cond2bm->filldown(cond, bm_cond);// set e.g. ut |
---|
| 119 | emp2bm->filldown(est_emp._point(), bm_cond);// set e.g. ut |
---|
| 120 | bm->bayes(dt, bm_cond); |
---|
| 121 | ll=bm->_ll(); |
---|
| 122 | } |
---|
| 123 | |
---|
[1077] | 124 | /*! Create object from the following structure |
---|
| 125 | |
---|
[1064] | 126 | \code |
---|
| 127 | class = "MarginalizedParticle"; |
---|
[1077] | 128 | parameter_pdf = configuration of bdm::epdf; % any offspring of epdf, bdm::epdf::from_setting |
---|
| 129 | --- inherited fields --- |
---|
| 130 | bdm::MarginalizedParticleBase::from_setting |
---|
[1064] | 131 | \endcode |
---|
[1077] | 132 | */ |
---|
[1064] | 133 | void from_setting(const Setting &set) { |
---|
| 134 | MarginalizedParticleBase::from_setting ( set ); |
---|
| 135 | par = UI::build<pdf> ( set, "parameter_pdf", UI::compulsory ); |
---|
| 136 | } |
---|
| 137 | |
---|
[1167] | 138 | void to_setting(Setting &set)const { |
---|
[1064] | 139 | MarginalizedParticleBase::to_setting(set); |
---|
| 140 | UI::save(par,set,"parameter_pdf"); |
---|
[1167] | 141 | UI::save(bm,set,"bm"); |
---|
[1064] | 142 | } |
---|
| 143 | void validate() { |
---|
| 144 | MarginalizedParticleBase::validate(); |
---|
| 145 | est_emp.set_rv(par->_rv()); |
---|
| 146 | if (est_emp.point.length()!=par->dimension()) |
---|
| 147 | est_emp.set_point(zeros(par->dimension())); |
---|
| 148 | est.validate(); |
---|
| 149 | |
---|
| 150 | yrv = bm->_yrv(); |
---|
| 151 | dimy = bm->dimensiony(); |
---|
| 152 | set_rv( concat(bm->_rv(), par->_rv())); |
---|
| 153 | set_dim( par->dimension()+bm->dimension()); |
---|
| 154 | |
---|
| 155 | rvc = par->_rvc(); |
---|
| 156 | rvc.add(bm->_rvc()); |
---|
| 157 | rvc=rvc.subt(par->_rv()); |
---|
| 158 | rvc=rvc.subt(par->_rv().copy_t(-1)); |
---|
| 159 | rvc=rvc.subt(bm->_rv().copy_t(-1)); // |
---|
| 160 | |
---|
| 161 | cond2bm=new datalink_part; |
---|
| 162 | cond2par=new datalink_part; |
---|
| 163 | emp2bm =new datalink_part; |
---|
| 164 | emp2par =new datalink_part; |
---|
| 165 | cond2bm->set_connection(bm->_rvc(), rvc); |
---|
| 166 | cond2par->set_connection(par->_rvc(), rvc); |
---|
| 167 | emp2bm->set_connection(bm->_rvc(), par->_rv()); |
---|
| 168 | emp2par->set_connection(par->_rvc(), par->_rv().copy_t(-1)); |
---|
| 169 | |
---|
| 170 | dimc = rvc._dsize(); |
---|
| 171 | }; |
---|
[887] | 172 | }; |
---|
| 173 | UIREGISTER(MarginalizedParticle); |
---|
| 174 | |
---|
[1077] | 175 | //! Internal class which is used in PF |
---|
[1064] | 176 | class BootstrapParticle : public BM { |
---|
| 177 | dirac est; |
---|
| 178 | shared_ptr<pdf> par; |
---|
| 179 | shared_ptr<pdf> obs; |
---|
| 180 | shared_ptr<datalink_part> cond2par; |
---|
| 181 | shared_ptr<datalink_part> cond2obs; |
---|
| 182 | shared_ptr<datalink_part> xt2obs; |
---|
| 183 | shared_ptr<datalink_part> xtm2par; |
---|
| 184 | public: |
---|
| 185 | BM* _copy() const { |
---|
| 186 | return new BootstrapParticle(*this); |
---|
| 187 | }; |
---|
| 188 | void bayes(const vec &dt, const vec &cond) { |
---|
| 189 | vec par_cond(par->dimensionc()); |
---|
| 190 | cond2par->filldown(cond,par_cond); // copy ut |
---|
| 191 | xtm2par->filldown(est._point(),par_cond); // copy xt-1 |
---|
| 192 | |
---|
| 193 | //sample new particle |
---|
| 194 | est.set_point(par->samplecond(par_cond)); |
---|
| 195 | //if (evalll) |
---|
| 196 | vec obs_cond(obs->dimensionc()); |
---|
| 197 | cond2obs->filldown(cond, obs_cond);// set e.g. ut |
---|
| 198 | xt2obs->filldown(est._point(), obs_cond);// set e.g. ut |
---|
| 199 | ll=obs->evallogcond(dt,obs_cond); |
---|
| 200 | } |
---|
| 201 | const dirac& posterior() const { |
---|
| 202 | return est; |
---|
| 203 | } |
---|
| 204 | |
---|
| 205 | void set_prior(const epdf *pdf0) { |
---|
| 206 | est.set_point(pdf0->sample()); |
---|
| 207 | } |
---|
| 208 | |
---|
[1077] | 209 | /*! Create object from the following structure |
---|
[1064] | 210 | \code |
---|
| 211 | class = "BootstrapParticle"; |
---|
[1077] | 212 | parameter_pdf = configuration of bdm::epdf; % any offspring of epdf, bdm::epdf::from_setting |
---|
| 213 | observation_pdf = configuration of bdm::epdf; % any offspring of epdf, bdm::epdf::from_setting |
---|
| 214 | --- inherited fields --- |
---|
| 215 | bdm::BM::from_setting |
---|
[1064] | 216 | \endcode |
---|
| 217 | */ |
---|
| 218 | void from_setting(const Setting &set) { |
---|
| 219 | BM::from_setting ( set ); |
---|
| 220 | par = UI::build<pdf> ( set, "parameter_pdf", UI::compulsory ); |
---|
| 221 | obs = UI::build<pdf> ( set, "observation_pdf", UI::compulsory ); |
---|
| 222 | } |
---|
[1077] | 223 | |
---|
[1064] | 224 | void validate() { |
---|
| 225 | yrv = obs->_rv(); |
---|
| 226 | dimy = obs->dimension(); |
---|
| 227 | set_rv( par->_rv()); |
---|
| 228 | set_dim( par->dimension()); |
---|
| 229 | |
---|
| 230 | rvc = par->_rvc().subt(par->_rv().copy_t(-1)); |
---|
| 231 | rvc.add(obs->_rvc()); // |
---|
| 232 | |
---|
| 233 | cond2obs=new datalink_part; |
---|
| 234 | cond2par=new datalink_part; |
---|
| 235 | xt2obs =new datalink_part; |
---|
| 236 | xtm2par =new datalink_part; |
---|
| 237 | cond2obs->set_connection(obs->_rvc(), rvc); |
---|
| 238 | cond2par->set_connection(par->_rvc(), rvc); |
---|
| 239 | xt2obs->set_connection(obs->_rvc(), _rv()); |
---|
| 240 | xtm2par->set_connection(par->_rvc(), _rv().copy_t(-1)); |
---|
| 241 | |
---|
| 242 | dimc = rvc._dsize(); |
---|
| 243 | }; |
---|
[887] | 244 | }; |
---|
| 245 | UIREGISTER(BootstrapParticle); |
---|
| 246 | |
---|
| 247 | |
---|
[8] | 248 | /*! |
---|
[32] | 249 | * \brief Trivial particle filter with proposal density equal to parameter evolution model. |
---|
[8] | 250 | |
---|
[32] | 251 | Posterior density is represented by a weighted empirical density (\c eEmp ). |
---|
[8] | 252 | */ |
---|
[32] | 253 | |
---|
| 254 | class PF : public BM { |
---|
[1064] | 255 | //! \var log_level_enums logweights |
---|
| 256 | //! all weightes will be logged |
---|
[870] | 257 | |
---|
[1064] | 258 | //! \var log_level_enums logmeans |
---|
| 259 | //! means of particles will be logged |
---|
[1196] | 260 | LOG_LEVEL(PF,logweights,logmeans,logvars,logNeff); |
---|
[1064] | 261 | |
---|
| 262 | class pf_mix: public emix_base { |
---|
| 263 | Array<BM*> &bms; |
---|
| 264 | public: |
---|
| 265 | pf_mix(vec &w0, Array<BM*> &bms0):emix_base(w0),bms(bms0) {} |
---|
| 266 | const epdf* component(const int &i)const { |
---|
| 267 | return &(bms(i)->posterior()); |
---|
| 268 | } |
---|
| 269 | int no_coms() const { |
---|
| 270 | return bms.length(); |
---|
| 271 | } |
---|
| 272 | }; |
---|
[8] | 273 | protected: |
---|
[1064] | 274 | //!number of particles; |
---|
| 275 | int n; |
---|
| 276 | //!posterior density |
---|
| 277 | pf_mix est; |
---|
| 278 | //! weights; |
---|
| 279 | vec w; |
---|
| 280 | //! particles |
---|
| 281 | Array<BM*> particles; |
---|
| 282 | //! internal structure storing loglikelihood of predictions |
---|
| 283 | vec lls; |
---|
[737] | 284 | |
---|
[1064] | 285 | //! which resampling method will be used |
---|
| 286 | RESAMPLING_METHOD resmethod; |
---|
| 287 | //! resampling threshold; in this case its meaning is minimum ratio of active particles |
---|
| 288 | //! For example, for 0.5 resampling is performed when the numebr of active aprticles drops belo 50%. |
---|
| 289 | double res_threshold; |
---|
[737] | 290 | |
---|
[1064] | 291 | //! \name Options |
---|
| 292 | //!@{ |
---|
| 293 | //!@} |
---|
[281] | 294 | |
---|
[8] | 295 | public: |
---|
[1064] | 296 | //! \name Constructors |
---|
| 297 | //!@{ |
---|
| 298 | PF ( ) : est(w,particles) { }; |
---|
[737] | 299 | |
---|
[1064] | 300 | void set_parameters ( int n0, double res_th0 = 0.5, RESAMPLING_METHOD rm = SYSTEMATIC ) { |
---|
| 301 | n = n0; |
---|
| 302 | res_threshold = res_th0; |
---|
| 303 | resmethod = rm; |
---|
| 304 | }; |
---|
| 305 | void set_model ( const BM *particle0, const epdf *prior) { |
---|
| 306 | if (n>0) { |
---|
| 307 | particles.set_length(n); |
---|
| 308 | for (int i=0; i<n; i++) { |
---|
| 309 | particles(i) = particle0->_copy(); |
---|
| 310 | particles(i)->set_prior(prior); |
---|
| 311 | } |
---|
| 312 | } |
---|
| 313 | // set values for posterior |
---|
| 314 | est.set_rv ( particle0->posterior()._rv() ); |
---|
| 315 | }; |
---|
| 316 | void set_statistics ( const vec w0, const epdf &epdf0 ) { |
---|
| 317 | //est.set_statistics ( w0, epdf0 ); |
---|
| 318 | }; |
---|
[1066] | 319 | /* void set_statistics ( const eEmp &epdf0 ) { |
---|
| 320 | bdm_assert_debug ( epdf0._rv().equal ( par->_rv() ), "Incompatible input" ); |
---|
| 321 | est = epdf0; |
---|
| 322 | };*/ |
---|
[1064] | 323 | //!@} |
---|
[850] | 324 | |
---|
[1064] | 325 | //! bayes compute weights of the |
---|
| 326 | virtual void bayes_weights(); |
---|
| 327 | //! important part of particle filtering - decide if it is time to perform resampling |
---|
| 328 | virtual bool do_resampling() { |
---|
| 329 | double eff = 1.0 / ( w * w ); |
---|
| 330 | return eff < ( res_threshold*n ); |
---|
| 331 | } |
---|
| 332 | void bayes ( const vec &yt, const vec &cond ); |
---|
| 333 | //!access function |
---|
| 334 | vec& _lls() { |
---|
| 335 | return lls; |
---|
| 336 | } |
---|
| 337 | //!access function |
---|
| 338 | RESAMPLING_METHOD _resmethod() const { |
---|
| 339 | return resmethod; |
---|
| 340 | } |
---|
| 341 | //! return correctly typed posterior (covariant return) |
---|
| 342 | const pf_mix& posterior() const { |
---|
| 343 | return est; |
---|
| 344 | } |
---|
[737] | 345 | |
---|
[1064] | 346 | /*! configuration structure for basic PF |
---|
| 347 | \code |
---|
[1085] | 348 | particle = bdm::BootstrapParticle; % one bayes rule for each point in the empirical support |
---|
| 349 | - or - = bdm::MarginalizedParticle; % (in case of Marginalized Particle filtering |
---|
| 350 | prior = epdf_class; % prior probability density on the empirical variable |
---|
[1064] | 351 | --- optional --- |
---|
[1085] | 352 | n = 10; % number of particles |
---|
[1064] | 353 | resmethod = 'systematic', or 'multinomial', or 'stratified' |
---|
[1085] | 354 | % resampling method |
---|
| 355 | res_threshold = 0.5; % resample when active particles drop below 50% |
---|
[1064] | 356 | \endcode |
---|
| 357 | */ |
---|
| 358 | void from_setting ( const Setting &set ) { |
---|
| 359 | BM::from_setting ( set ); |
---|
| 360 | UI::get ( log_level, set, "log_level", UI::optional ); |
---|
[665] | 361 | |
---|
[1064] | 362 | shared_ptr<BM> bm0 = UI::build<BM>(set, "particle",UI::compulsory); |
---|
[737] | 363 | |
---|
[1064] | 364 | n =0; |
---|
| 365 | UI::get(n,set,"n",UI::optional);; |
---|
| 366 | if (n>0) { |
---|
| 367 | particles.set_length(n); |
---|
| 368 | for(int i=0; i<n; i++) { |
---|
| 369 | particles(i)=bm0->_copy(); |
---|
| 370 | } |
---|
| 371 | w = ones(n)/n; |
---|
| 372 | } |
---|
| 373 | shared_ptr<epdf> pri = UI::build<epdf>(set,"prior"); |
---|
| 374 | set_prior(pri.get()); |
---|
| 375 | // set resampling method |
---|
| 376 | resmethod_from_set ( set ); |
---|
| 377 | //set drv |
---|
[737] | 378 | |
---|
[1064] | 379 | rvc = bm0->_rvc(); |
---|
| 380 | dimc = bm0->dimensionc(); |
---|
| 381 | BM::set_rv(bm0->_rv()); |
---|
| 382 | yrv=bm0->_yrv(); |
---|
| 383 | dimy = bm0->dimensiony(); |
---|
| 384 | } |
---|
[1167] | 385 | |
---|
| 386 | void to_setting(Setting &set) const{ |
---|
| 387 | BM::to_setting(set); |
---|
| 388 | UI::save(particles, set,"particles"); |
---|
| 389 | UI::save(w,set,"w"); |
---|
| 390 | } |
---|
[766] | 391 | |
---|
[1064] | 392 | void log_register ( bdm::logger& L, const string& prefix ) { |
---|
| 393 | BM::log_register(L,prefix); |
---|
[1196] | 394 | if (log_level[logweights]) { |
---|
| 395 | L.add_vector( log_level, logweights, RV ( particles.length()), prefix); |
---|
| 396 | } |
---|
| 397 | if (log_level[logNeff]) { |
---|
| 398 | L.add_vector( log_level, logNeff, RV ( 1), prefix); |
---|
| 399 | } |
---|
| 400 | if (log_level[logmeans]) { |
---|
[1064] | 401 | for (int i=0; i<particles.length(); i++) { |
---|
| 402 | L.add_vector( log_level, logmeans, RV ( particles(i)->dimension() ), prefix , i); |
---|
| 403 | } |
---|
| 404 | } |
---|
| 405 | if (log_level[logvars]) { |
---|
| 406 | for (int i=0; i<particles.length(); i++) { |
---|
| 407 | L.add_vector( log_level, logvars, RV ( particles(i)->dimension() ), prefix , i); |
---|
| 408 | } |
---|
| 409 | } |
---|
| 410 | }; |
---|
| 411 | void log_write ( ) const { |
---|
| 412 | BM::log_write(); |
---|
[1187] | 413 | // weigths are before resamplign -- bayes |
---|
[1064] | 414 | if (log_level[logmeans]) { |
---|
| 415 | for (int i=0; i<particles.length(); i++) { |
---|
| 416 | log_level.store( logmeans, particles(i)->posterior().mean(), i); |
---|
| 417 | } |
---|
| 418 | } |
---|
| 419 | if (log_level[logvars]) { |
---|
| 420 | for (int i=0; i<particles.length(); i++) { |
---|
| 421 | log_level.store( logvars, particles(i)->posterior().variance(), i); |
---|
| 422 | } |
---|
| 423 | } |
---|
| 424 | |
---|
| 425 | } |
---|
| 426 | |
---|
| 427 | void set_prior(const epdf *pri) { |
---|
| 428 | const emix_base *emi=dynamic_cast<const emix_base*>(pri); |
---|
| 429 | if (emi) { |
---|
| 430 | bdm_assert(particles.length()>0, "initial particle is not assigned"); |
---|
| 431 | n = emi->_w().length(); |
---|
| 432 | int old_n = particles.length(); |
---|
| 433 | if (n!=old_n) { |
---|
| 434 | particles.set_length(n,true); |
---|
| 435 | } |
---|
| 436 | for(int i=old_n; i<n; i++) { |
---|
| 437 | particles(i)=particles(0)->_copy(); |
---|
| 438 | } |
---|
| 439 | |
---|
| 440 | for (int i =0; i<n; i++) { |
---|
| 441 | particles(i)->set_prior(emi->_com(i)); |
---|
| 442 | } |
---|
| 443 | } else { |
---|
| 444 | // try to find "n" |
---|
| 445 | bdm_assert(n>0, "Field 'n' must be filled when prior is not of type emix"); |
---|
| 446 | for (int i =0; i<n; i++) { |
---|
| 447 | particles(i)->set_prior(pri); |
---|
| 448 | } |
---|
| 449 | |
---|
| 450 | } |
---|
| 451 | } |
---|
| 452 | //! auxiliary function reading parameter 'resmethod' from configuration file |
---|
| 453 | void resmethod_from_set ( const Setting &set ) { |
---|
| 454 | string resmeth; |
---|
| 455 | if ( UI::get ( resmeth, set, "resmethod", UI::optional ) ) { |
---|
| 456 | if ( resmeth == "systematic" ) { |
---|
| 457 | resmethod = SYSTEMATIC; |
---|
| 458 | } else { |
---|
| 459 | if ( resmeth == "multinomial" ) { |
---|
| 460 | resmethod = MULTINOMIAL; |
---|
| 461 | } else { |
---|
| 462 | if ( resmeth == "stratified" ) { |
---|
| 463 | resmethod = STRATIFIED; |
---|
| 464 | } else { |
---|
| 465 | bdm_error ( "Unknown resampling method" ); |
---|
| 466 | } |
---|
| 467 | } |
---|
| 468 | } |
---|
| 469 | } else { |
---|
| 470 | resmethod = SYSTEMATIC; |
---|
| 471 | }; |
---|
| 472 | if ( !UI::get ( res_threshold, set, "res_threshold", UI::optional ) ) { |
---|
| 473 | res_threshold = 0.9; |
---|
| 474 | } |
---|
| 475 | //validate(); |
---|
| 476 | } |
---|
| 477 | |
---|
| 478 | void validate() { |
---|
| 479 | BM::validate(); |
---|
| 480 | est.validate(); |
---|
| 481 | bdm_assert ( n>0, "empty particle pool" ); |
---|
| 482 | n = w.length(); |
---|
| 483 | lls = zeros ( n ); |
---|
| 484 | |
---|
| 485 | if ( particles(0)->_rv()._dsize() > 0 ) { |
---|
| 486 | bdm_assert ( particles(0)->_rv()._dsize() == est.dimension(), "MPF:: Mismatch of RV " +particles(0)->_rv().to_string() + |
---|
| 487 | " of size (" +num2str(particles(0)->_rv()._dsize())+") and dimension of posterior ("+num2str(est.dimension()) + ")" ); |
---|
| 488 | } |
---|
| 489 | } |
---|
| 490 | //! resample posterior density (from outside - see MPF) |
---|
| 491 | void resample ( ) { |
---|
| 492 | ivec ind = zeros_i ( n ); |
---|
| 493 | bdm::resample(w,ind,resmethod); |
---|
| 494 | // copy the internals according to ind |
---|
| 495 | for (int i = 0; i < n; i++ ) { |
---|
| 496 | if ( ind ( i ) != i ) { |
---|
| 497 | delete particles(i); |
---|
| 498 | particles( i ) = particles( ind ( i ) )->_copy(); |
---|
| 499 | } |
---|
| 500 | w ( i ) = 1.0 / n; |
---|
| 501 | } |
---|
| 502 | } |
---|
| 503 | //! access function |
---|
| 504 | Array<BM*>& _particles() { |
---|
| 505 | return particles; |
---|
| 506 | } |
---|
| 507 | ~PF() { |
---|
| 508 | for (int i=0; i<particles.length(); i++) { |
---|
| 509 | delete particles(i); |
---|
| 510 | } |
---|
| 511 | } |
---|
| 512 | |
---|
[8] | 513 | }; |
---|
[737] | 514 | UIREGISTER ( PF ); |
---|
[8] | 515 | |
---|
[1064] | 516 | /*! Marginalized particle for state-space models with unknown parameters of distribuution of residues on \f$v_t\f$. |
---|
[8] | 517 | |
---|
[979] | 518 | \f{eqnarray*}{ |
---|
[1066] | 519 | x_t &=& g(x_{t-1}) + v_t,\\ |
---|
| 520 | y_t &\sim &fy(x_t), |
---|
| 521 | \f} |
---|
[1064] | 522 | |
---|
[1066] | 523 | This particle is a only a shell creating the residues calling internal estimator of their parameters. The internal estimator can be of any compatible type, e.g. ARX for Gaussian residues with unknown mean and variance. |
---|
[1064] | 524 | |
---|
[1066] | 525 | */ |
---|
[1064] | 526 | class NoiseParticleX : public MarginalizedParticleBase { |
---|
| 527 | protected: |
---|
| 528 | //! function transforming xt, ut -> x_t+1 |
---|
| 529 | shared_ptr<fnc> g; // pdf for non-linear part |
---|
| 530 | //! function transforming xt,ut -> yt |
---|
| 531 | shared_ptr<pdf> fy; // pdf for non-linear part |
---|
[8] | 532 | |
---|
[1064] | 533 | RV rvx; |
---|
| 534 | RV rvxc; |
---|
| 535 | RV rvyc; |
---|
| 536 | |
---|
| 537 | //!link from condition to f |
---|
| 538 | datalink_part cond2g; |
---|
| 539 | //!link from condition to h |
---|
| 540 | datalink_part cond2fy; |
---|
| 541 | //!link from xt to f |
---|
| 542 | datalink_part x2g; |
---|
| 543 | //!link from xt to h |
---|
| 544 | datalink_part x2fy; |
---|
| 545 | |
---|
| 546 | public: |
---|
| 547 | BM* _copy() const { |
---|
| 548 | return new NoiseParticleX(*this); |
---|
| 549 | }; |
---|
| 550 | void bayes(const vec &dt, const vec &cond) { |
---|
[1167] | 551 | //shared_ptr<epdf> pred_v=bm->epredictor(); |
---|
[1064] | 552 | |
---|
[1196] | 553 | vec xt; |
---|
| 554 | vec xthat; |
---|
| 555 | //vec vt=pred_v->sample(); |
---|
| 556 | // vec vt = bm->samplepred(); |
---|
| 557 | if (1){ |
---|
| 558 | //vec vt=pred_v->sample(); |
---|
| 559 | int dim_x = bm->dimensiony(); |
---|
| 560 | int dim_y = dimensiony(); |
---|
| 561 | |
---|
| 562 | vec post_mean =bm->posterior().mean(); |
---|
| 563 | mat SigV; |
---|
| 564 | if (post_mean.length()==dim_x){ |
---|
| 565 | SigV=diag(vec(post_mean._data(), dim_x)); |
---|
| 566 | } else { |
---|
| 567 | SigV=mat(post_mean._data(), dim_x, dim_x); |
---|
| 568 | } |
---|
[1064] | 569 | |
---|
[1196] | 570 | vec &xtm=est_emp.point; |
---|
| 571 | vec g_args(g->dimensionc()); |
---|
| 572 | x2g.filldown(xtm,g_args); |
---|
| 573 | cond2g.filldown(cond,g_args); |
---|
| 574 | xthat = g->eval(g_args); |
---|
| 575 | |
---|
| 576 | mat IC = eye(dim_x+dim_y); |
---|
| 577 | IC.set_submatrix(0,0,eye(dim_x)); |
---|
| 578 | mat SigJo=eye(dim_x + dim_y); |
---|
| 579 | SigJo.set_submatrix(0,0,SigV); |
---|
| 580 | mlnorm<ldmat>* mlfy=dynamic_cast<mlnorm<ldmat>* >(fy.get()); |
---|
| 581 | |
---|
| 582 | SigJo.set_submatrix(dim_x,dim_x, mlfy->_R()); |
---|
| 583 | IC.set_submatrix(dim_x,0,-mlfy->_A()); |
---|
| 584 | |
---|
| 585 | mat iIC=inv(IC); |
---|
| 586 | enorm<chmat> Jo; |
---|
| 587 | Jo._mu()=iIC*concat(xthat, zeros(dim_y)); |
---|
| 588 | Jo._R()=iIC*SigJo*iIC.T(); |
---|
| 589 | Jo.set_rv(concat(bm->_yrv(),yrv)); |
---|
| 590 | Jo.validate(); |
---|
| 591 | shared_ptr<pdf> Cond=Jo.condition(bm->_yrv()); |
---|
| 592 | //new sample |
---|
| 593 | |
---|
| 594 | xt = Cond->samplecond(dt);// + vt; |
---|
| 595 | } else{ |
---|
| 596 | |
---|
| 597 | |
---|
[1064] | 598 | //new sample |
---|
| 599 | vec &xtm=est_emp.point; |
---|
| 600 | vec g_args(g->dimensionc()); |
---|
| 601 | x2g.filldown(xtm,g_args); |
---|
| 602 | cond2g.filldown(cond,g_args); |
---|
[1196] | 603 | xthat = g->eval(g_args); |
---|
| 604 | |
---|
| 605 | xt = xthat + bm->samplepred(); |
---|
| 606 | |
---|
| 607 | } |
---|
[1064] | 608 | est_emp.point=xt; |
---|
| 609 | |
---|
| 610 | // the vector [v_t] updates bm, |
---|
[1196] | 611 | bm->bayes(xt-xthat); |
---|
[1064] | 612 | |
---|
| 613 | // residue of observation |
---|
| 614 | vec fy_args(fy->dimensionc()); |
---|
| 615 | x2fy.filldown(xt,fy_args); |
---|
| 616 | cond2fy.filldown(cond,fy_args); |
---|
| 617 | |
---|
[1167] | 618 | ll= fy->evallogcond(dt,fy_args); |
---|
[1064] | 619 | } |
---|
| 620 | void from_setting(const Setting &set) { |
---|
| 621 | MarginalizedParticleBase::from_setting(set); //reads bm, yrv,rvc, bm_rv, etc... |
---|
| 622 | |
---|
| 623 | g=UI::build<fnc>(set,"g",UI::compulsory); |
---|
| 624 | fy=UI::build<pdf>(set,"fy",UI::compulsory); |
---|
| 625 | UI::get(rvx,set,"rvx",UI::compulsory); |
---|
| 626 | est_emp.set_rv(rvx); |
---|
| 627 | |
---|
| 628 | UI::get(rvxc,set,"rvxc",UI::compulsory); |
---|
| 629 | UI::get(rvyc,set,"rvyc",UI::compulsory); |
---|
| 630 | |
---|
| 631 | } |
---|
[1167] | 632 | void to_setting (Setting &set) const { |
---|
| 633 | MarginalizedParticleBase::to_setting(set); //reads bm, yrv,rvc, bm_rv, etc... |
---|
| 634 | UI::save(g,set,"g"); |
---|
| 635 | UI::save(fy,set,"fy"); |
---|
| 636 | UI::save(bm,set,"bm"); |
---|
| 637 | } |
---|
[1064] | 638 | void validate() { |
---|
| 639 | MarginalizedParticleBase::validate(); |
---|
| 640 | |
---|
| 641 | dimy = fy->dimension(); |
---|
| 642 | bm->set_yrv(rvx); |
---|
| 643 | |
---|
| 644 | est_emp.set_rv(rvx); |
---|
| 645 | est_emp.set_dim(rvx._dsize()); |
---|
| 646 | est.validate(); |
---|
| 647 | // |
---|
| 648 | //check dimensions |
---|
| 649 | rvc = rvxc.subt(rvx.copy_t(-1)); |
---|
| 650 | rvc.add( rvyc); |
---|
| 651 | rvc=rvc.subt(rvx); |
---|
| 652 | |
---|
| 653 | bdm_assert(g->dimension()==rvx._dsize(),"rvx is not described"); |
---|
| 654 | bdm_assert(g->dimensionc()==rvxc._dsize(),"rvxc is not described"); |
---|
| 655 | bdm_assert(fy->dimensionc()==rvyc._dsize(),"rvyc is not described"); |
---|
| 656 | |
---|
| 657 | bdm_assert(bm->dimensiony()==g->dimension(), |
---|
| 658 | "Incompatible noise estimator of dimension " + |
---|
| 659 | num2str(bm->dimensiony()) + " does not match dimension of g , " + |
---|
| 660 | num2str(g->dimension())); |
---|
| 661 | |
---|
| 662 | dimc = rvc._dsize(); |
---|
| 663 | |
---|
| 664 | //establish datalinks |
---|
| 665 | x2g.set_connection(rvxc, rvx.copy_t(-1)); |
---|
| 666 | cond2g.set_connection(rvxc, rvc); |
---|
| 667 | |
---|
| 668 | x2fy.set_connection(rvyc, rvx); |
---|
| 669 | cond2fy.set_connection(rvyc, rvc); |
---|
| 670 | } |
---|
[979] | 671 | }; |
---|
| 672 | UIREGISTER(NoiseParticleX); |
---|
[870] | 673 | |
---|
[1192] | 674 | |
---|
[1170] | 675 | /*! Marginalized particle for state-space models with unknown parameters of distribuution of residues on \f$v_t\f$ and \f$ w_t \f$. |
---|
| 676 | |
---|
| 677 | \f{eqnarray*}{ |
---|
| 678 | x_t &=& g(x_{t-1}) + v_t,\\ |
---|
| 679 | y_t &= &h(x_t)+w_t, |
---|
| 680 | \f} |
---|
| 681 | |
---|
| 682 | This particle is a only a shell creating the residues calling internal estimator of their parameters. The internal estimator can be of any compatible type, e.g. ARX for Gaussian residues with unknown mean and variance. |
---|
| 683 | |
---|
| 684 | */ |
---|
| 685 | class NoiseParticleXY : public BM { |
---|
[1196] | 686 | protected: |
---|
[1170] | 687 | //! discrte particle |
---|
| 688 | dirac est_emp; |
---|
| 689 | //! internal Bayes Model |
---|
| 690 | shared_ptr<BM> bmx; |
---|
| 691 | shared_ptr<BM> bmy; |
---|
| 692 | |
---|
| 693 | //! \brief Internal class for custom posterior - product of empirical and exact part |
---|
| 694 | class eprod_3:public eprod_base { |
---|
| 695 | protected: |
---|
| 696 | NoiseParticleXY ∓ |
---|
| 697 | public: |
---|
| 698 | eprod_3(NoiseParticleXY &m):mp(m) {} |
---|
| 699 | const epdf* factor(int i) const { |
---|
| 700 | if (i==0) return &mp.bmx->posterior() ; |
---|
| 701 | if(i==1) return &mp.bmy->posterior(); |
---|
| 702 | return &mp.est_emp; |
---|
| 703 | } |
---|
| 704 | const int no_factors() const { |
---|
| 705 | return 3; |
---|
| 706 | } |
---|
| 707 | } est; |
---|
| 708 | |
---|
| 709 | protected: |
---|
| 710 | //! function transforming xt, ut -> x_t+1 |
---|
| 711 | shared_ptr<fnc> g; // pdf for non-linear part |
---|
| 712 | //! function transforming xt,ut -> yt |
---|
| 713 | shared_ptr<fnc> h; // pdf for non-linear part |
---|
| 714 | |
---|
| 715 | RV rvx; |
---|
| 716 | RV rvxc; |
---|
| 717 | RV rvyc; |
---|
| 718 | |
---|
| 719 | //!link from condition to f |
---|
| 720 | datalink_part cond2g; |
---|
| 721 | //!link from condition to h |
---|
| 722 | datalink_part cond2h; |
---|
| 723 | //!link from xt to f |
---|
| 724 | datalink_part x2g; |
---|
| 725 | //!link from xt to h |
---|
| 726 | datalink_part x2h; |
---|
| 727 | |
---|
| 728 | public: |
---|
| 729 | NoiseParticleXY():est(*this) {}; |
---|
[1177] | 730 | NoiseParticleXY(const NoiseParticleXY &m2):BM(m2),est(*this),g(m2.g),h(m2.h), rvx(m2.rvx),rvxc(m2.rvxc),rvyc(m2.rvyc) { |
---|
[1170] | 731 | bmx = m2.bmx->_copy(); |
---|
| 732 | bmy = m2.bmy->_copy(); |
---|
| 733 | est_emp = m2.est_emp; |
---|
| 734 | //est.validate(); |
---|
| 735 | validate(); |
---|
| 736 | }; |
---|
| 737 | |
---|
| 738 | const eprod_3& posterior() const { |
---|
| 739 | return est; |
---|
| 740 | } |
---|
| 741 | |
---|
| 742 | void set_prior(const epdf *pdf0) { |
---|
| 743 | const eprod *ep=dynamic_cast<const eprod*>(pdf0); |
---|
| 744 | if (ep) { // full prior |
---|
| 745 | bdm_assert(ep->no_factors()==2,"Incompatible prod"); |
---|
| 746 | bmx->set_prior(ep->factor(0)); |
---|
| 747 | bmy->set_prior(ep->factor(1)); |
---|
| 748 | est_emp.set_point(ep->factor(2)->sample()); |
---|
| 749 | } else { |
---|
| 750 | // assume prior is only for emp; |
---|
| 751 | est_emp.set_point(pdf0->sample()); |
---|
| 752 | } |
---|
| 753 | } |
---|
| 754 | |
---|
| 755 | BM* _copy() const { |
---|
| 756 | return new NoiseParticleXY(*this); |
---|
| 757 | }; |
---|
| 758 | |
---|
| 759 | void bayes(const vec &dt, const vec &cond) { |
---|
| 760 | //shared_ptr<epdf> pred_v=bm->epredictor(); |
---|
| 761 | |
---|
| 762 | //vec vt=pred_v->sample(); |
---|
| 763 | vec vt = bmx->samplepred(); |
---|
| 764 | |
---|
| 765 | //new sample |
---|
| 766 | vec &xtm=est_emp.point; |
---|
| 767 | vec g_args(g->dimensionc()); |
---|
| 768 | x2g.filldown(xtm,g_args); |
---|
| 769 | cond2g.filldown(cond,g_args); |
---|
| 770 | vec xt = g->eval(g_args) + vt; |
---|
| 771 | est_emp.point=xt; |
---|
| 772 | |
---|
| 773 | // the vector [v_t] updates bm, |
---|
| 774 | bmx->bayes(vt); |
---|
| 775 | |
---|
| 776 | // residue of observation |
---|
| 777 | vec h_args(h->dimensionc()); |
---|
| 778 | x2h.filldown(xt,h_args); |
---|
| 779 | cond2h.filldown(cond,h_args); |
---|
[1177] | 780 | |
---|
| 781 | vec z_y =h->eval(h_args)-dt; |
---|
[1187] | 782 | // ARX *abm = dynamic_cast<ARX*>(bmy.get()); |
---|
| 783 | /* double ll2; |
---|
[1177] | 784 | if (abm){ //ARX |
---|
| 785 | shared_ptr<epdf> pr_y(abm->epredictor_student(empty_vec)); |
---|
| 786 | ll2=pr_y->evallog(z_y); |
---|
| 787 | } else{ |
---|
| 788 | shared_ptr<epdf> pr_y(bmy->epredictor(empty_vec)); |
---|
| 789 | ll2=pr_y->evallog(z_y); |
---|
[1187] | 790 | }*/ |
---|
[1170] | 791 | |
---|
[1177] | 792 | bmy->bayes(z_y); |
---|
| 793 | // test _lls |
---|
[1187] | 794 | ll= bmy->_ll(); |
---|
[1170] | 795 | } |
---|
| 796 | void from_setting(const Setting &set) { |
---|
| 797 | BM::from_setting(set); //reads bm, yrv,rvc, bm_rv, etc... |
---|
| 798 | bmx = UI::build<BM>(set,"bmx",UI::compulsory); |
---|
| 799 | |
---|
| 800 | bmy = UI::build<BM>(set,"bmy",UI::compulsory); |
---|
| 801 | g=UI::build<fnc>(set,"g",UI::compulsory); |
---|
| 802 | h=UI::build<fnc>(set,"h",UI::compulsory); |
---|
| 803 | UI::get(rvx,set,"rvx",UI::compulsory); |
---|
| 804 | est_emp.set_rv(rvx); |
---|
| 805 | |
---|
| 806 | UI::get(rvxc,set,"rvxc",UI::compulsory); |
---|
| 807 | UI::get(rvyc,set,"rvyc",UI::compulsory); |
---|
| 808 | |
---|
| 809 | } |
---|
| 810 | void to_setting (Setting &set) const { |
---|
| 811 | BM::to_setting(set); //reads bm, yrv,rvc, bm_rv, etc... |
---|
| 812 | UI::save(g,set,"g"); |
---|
| 813 | UI::save(h,set,"h"); |
---|
| 814 | UI::save(bmx,set,"bmx"); |
---|
| 815 | UI::save(bmy,set,"bmy"); |
---|
| 816 | } |
---|
| 817 | void validate() { |
---|
| 818 | BM::validate(); |
---|
| 819 | |
---|
| 820 | dimy = h->dimension(); |
---|
| 821 | // bmx->set_yrv(rvx); |
---|
| 822 | // bmy-> |
---|
| 823 | |
---|
| 824 | est_emp.set_rv(rvx); |
---|
| 825 | est_emp.set_dim(rvx._dsize()); |
---|
| 826 | est.validate(); |
---|
| 827 | // |
---|
| 828 | //check dimensions |
---|
| 829 | rvc = rvxc.subt(rvx.copy_t(-1)); |
---|
| 830 | rvc.add( rvyc); |
---|
| 831 | rvc=rvc.subt(rvx); |
---|
| 832 | |
---|
| 833 | bdm_assert(g->dimension()==rvx._dsize(),"rvx is not described"); |
---|
| 834 | bdm_assert(g->dimensionc()==rvxc._dsize(),"rvxc is not described"); |
---|
| 835 | bdm_assert(h->dimensionc()==rvyc._dsize(),"rvyc is not described"); |
---|
| 836 | |
---|
| 837 | bdm_assert(bmx->dimensiony()==g->dimension(), |
---|
| 838 | "Incompatible noise estimator of dimension " + |
---|
| 839 | num2str(bmx->dimensiony()) + " does not match dimension of g , " + |
---|
| 840 | num2str(g->dimension()) |
---|
| 841 | ); |
---|
| 842 | |
---|
| 843 | dimc = rvc._dsize(); |
---|
| 844 | |
---|
| 845 | //establish datalinks |
---|
| 846 | x2g.set_connection(rvxc, rvx.copy_t(-1)); |
---|
| 847 | cond2g.set_connection(rvxc, rvc); |
---|
| 848 | |
---|
| 849 | x2h.set_connection(rvyc, rvx); |
---|
| 850 | cond2h.set_connection(rvyc, rvc); |
---|
| 851 | } |
---|
| 852 | |
---|
| 853 | }; |
---|
| 854 | UIREGISTER(NoiseParticleXY); |
---|
| 855 | |
---|
[1196] | 856 | class NoiseParticleXYprop: public NoiseParticleXY{ |
---|
| 857 | public: |
---|
| 858 | BM* _copy() const { |
---|
| 859 | return new NoiseParticleXYprop(*this); |
---|
| 860 | }; |
---|
| 861 | void bayes(const vec &dt, const vec &cond) { |
---|
| 862 | int dim_x = bmx->dimensiony(); |
---|
| 863 | int dim_y = bmy->dimensiony(); |
---|
| 864 | |
---|
| 865 | //vec vt=pred_v->sample(); |
---|
| 866 | mat SigV=mat(bmx->posterior().mean()._data(), dim_x, dim_x); |
---|
| 867 | mat SigW=mat(bmy->posterior().mean()._data(), dim_y, dim_y); |
---|
| 868 | |
---|
| 869 | vec &xtm=est_emp.point; |
---|
| 870 | vec g_args(g->dimensionc()); |
---|
| 871 | x2g.filldown(xtm,g_args); |
---|
| 872 | cond2g.filldown(cond,g_args); |
---|
| 873 | vec xthat = g->eval(g_args); |
---|
| 874 | |
---|
| 875 | mat IC = eye(dim_x+dim_y); |
---|
| 876 | IC.set_submatrix(0,0,eye(dim_x)); |
---|
| 877 | mat SigJo=eye(dim_x + dim_y); |
---|
| 878 | SigJo.set_submatrix(0,0,SigV); |
---|
| 879 | SigJo.set_submatrix(dim_x,dim_x, 10*SigW); |
---|
| 880 | |
---|
| 881 | diffbifn* hb=dynamic_cast<bilinfn*>(h.get()); |
---|
| 882 | mat C=zeros(dim_y,dim_x); |
---|
| 883 | if (hb){ |
---|
| 884 | hb->dfdx_cond(xtm, cond, C,true); |
---|
| 885 | IC.set_submatrix(dim_x,0,-C); |
---|
| 886 | } |
---|
| 887 | |
---|
| 888 | mat iIC=inv(IC); |
---|
| 889 | enorm<chmat> Jo; |
---|
| 890 | Jo._mu()=iIC*concat(xthat, zeros(dim_y)); |
---|
| 891 | Jo._R()=iIC*SigJo*iIC.T(); |
---|
| 892 | Jo.set_rv(concat(bmx->_yrv(), bmy->_yrv())); |
---|
| 893 | Jo.validate(); |
---|
| 894 | shared_ptr<pdf> Cond=Jo.condition(bmx->_yrv()); |
---|
| 895 | //new sample |
---|
| 896 | |
---|
| 897 | vec xt = Cond->samplecond(dt);// + vt; |
---|
| 898 | |
---|
| 899 | est_emp.point=xt; |
---|
| 900 | |
---|
| 901 | // the vector [v_t] updates bm, |
---|
| 902 | bmx->bayes(xt-xthat); |
---|
| 903 | |
---|
| 904 | // residue of observation |
---|
| 905 | vec h_args(h->dimensionc()); |
---|
| 906 | x2h.filldown(xt,h_args); |
---|
| 907 | cond2h.filldown(cond,h_args); |
---|
| 908 | |
---|
| 909 | vec z_y =h->eval(h_args)-dt; |
---|
| 910 | // ARX *abm = dynamic_cast<ARX*>(bmy.get()); |
---|
| 911 | /* double ll2; |
---|
| 912 | i f (abm){ //ARX* |
---|
| 913 | shared_ptr<epdf> pr_y(abm->epredictor_student(empty_vec)); |
---|
| 914 | ll2=pr_y->evallog(z_y); |
---|
| 915 | } else{ |
---|
| 916 | shared_ptr<epdf> pr_y(bmy->epredictor(empty_vec)); |
---|
| 917 | ll2=pr_y->evallog(z_y); |
---|
| 918 | }*/ |
---|
| 919 | |
---|
| 920 | bmy->bayes(z_y); |
---|
| 921 | // test _lls |
---|
| 922 | ll= bmy->_ll(); |
---|
| 923 | } |
---|
| 924 | }; |
---|
| 925 | UIREGISTER(NoiseParticleXYprop); |
---|
| 926 | |
---|
[979] | 927 | /*! Marginalized particle for state-space models with unknown parameters of residues distribution |
---|
| 928 | |
---|
| 929 | \f{eqnarray*}{ |
---|
[1066] | 930 | x_t &=& g(x_{t-1}) + v_t,\\ |
---|
| 931 | z_t &= &h(x_{t-1}) + w_t, |
---|
| 932 | \f} |
---|
[1064] | 933 | |
---|
[1066] | 934 | This particle is a only a shell creating the residues calling internal estimator of their parameters. The internal estimator can be of any compatible type, e.g. ARX for Gaussian residues with unknown mean and variance. |
---|
[1064] | 935 | |
---|
[1066] | 936 | */ |
---|
[1064] | 937 | class NoiseParticle : public MarginalizedParticleBase { |
---|
| 938 | protected: |
---|
| 939 | //! function transforming xt, ut -> x_t+1 |
---|
| 940 | shared_ptr<fnc> g; // pdf for non-linear part |
---|
| 941 | //! function transforming xt,ut -> yt |
---|
| 942 | shared_ptr<fnc> h; // pdf for non-linear part |
---|
| 943 | |
---|
| 944 | RV rvx; |
---|
| 945 | RV rvxc; |
---|
| 946 | RV rvyc; |
---|
| 947 | |
---|
| 948 | //!link from condition to f |
---|
| 949 | datalink_part cond2g; |
---|
| 950 | //!link from condition to h |
---|
| 951 | datalink_part cond2h; |
---|
| 952 | //!link from xt to f |
---|
| 953 | datalink_part x2g; |
---|
| 954 | //!link from xt to h |
---|
| 955 | datalink_part x2h; |
---|
| 956 | |
---|
| 957 | public: |
---|
| 958 | BM* _copy() const { |
---|
| 959 | return new NoiseParticle(*this); |
---|
| 960 | }; |
---|
| 961 | void bayes(const vec &dt, const vec &cond) { |
---|
| 962 | shared_ptr<epdf> pred_vw=bm->epredictor(); |
---|
| 963 | shared_ptr<epdf> pred_v = pred_vw->marginal(rvx); |
---|
| 964 | |
---|
| 965 | vec vt=pred_v->sample(); |
---|
| 966 | |
---|
| 967 | //new sample |
---|
| 968 | vec &xtm=est_emp.point; |
---|
| 969 | vec g_args(g->dimensionc()); |
---|
| 970 | x2g.filldown(xtm,g_args); |
---|
| 971 | cond2g.filldown(cond,g_args); |
---|
| 972 | vec xt = g->eval(g_args) + vt; |
---|
| 973 | est_emp.point=xt; |
---|
| 974 | |
---|
| 975 | // residue of observation |
---|
| 976 | vec h_args(h->dimensionc()); |
---|
| 977 | x2h.filldown(xt,h_args); |
---|
| 978 | cond2h.filldown(cond,h_args); |
---|
| 979 | vec wt = dt-h->eval(h_args); |
---|
| 980 | // the vector [v_t,w_t] is now complete |
---|
| 981 | bm->bayes(concat(vt,wt)); |
---|
| 982 | ll=bm->_ll(); |
---|
| 983 | } |
---|
| 984 | void from_setting(const Setting &set) { |
---|
| 985 | MarginalizedParticleBase::from_setting(set); //reads bm, yrv,rvc, bm_rv, etc... |
---|
| 986 | |
---|
| 987 | UI::get(g,set,"g",UI::compulsory); |
---|
| 988 | UI::get(h,set,"h",UI::compulsory); |
---|
| 989 | UI::get(rvx,set,"rvx",UI::compulsory); |
---|
| 990 | est_emp.set_rv(rvx); |
---|
| 991 | |
---|
| 992 | UI::get(rvxc,set,"rvxc",UI::compulsory); |
---|
| 993 | UI::get(rvyc,set,"rvyc",UI::compulsory); |
---|
| 994 | |
---|
| 995 | } |
---|
| 996 | void validate() { |
---|
| 997 | MarginalizedParticleBase::validate(); |
---|
| 998 | |
---|
| 999 | dimy = h->dimension(); |
---|
| 1000 | bm->set_yrv(concat(rvx,yrv)); |
---|
| 1001 | |
---|
| 1002 | est_emp.set_rv(rvx); |
---|
| 1003 | est_emp.set_dim(rvx._dsize()); |
---|
| 1004 | est.validate(); |
---|
| 1005 | // |
---|
| 1006 | //check dimensions |
---|
| 1007 | rvc = rvxc.subt(rvx.copy_t(-1)); |
---|
| 1008 | rvc.add( rvyc); |
---|
| 1009 | rvc=rvc.subt(rvx); |
---|
| 1010 | |
---|
| 1011 | bdm_assert(g->dimension()==rvx._dsize(),"rvx is not described"); |
---|
| 1012 | bdm_assert(g->dimensionc()==rvxc._dsize(),"rvxc is not described"); |
---|
| 1013 | bdm_assert(h->dimension()==rvyc._dsize(),"rvyc is not described"); |
---|
| 1014 | |
---|
| 1015 | bdm_assert(bm->dimensiony()==g->dimension()+h->dimension(), |
---|
| 1016 | "Incompatible noise estimator of dimension " + |
---|
| 1017 | num2str(bm->dimensiony()) + " does not match dimension of g and h, " + |
---|
| 1018 | num2str(g->dimension())+" and "+ num2str(h->dimension()) ); |
---|
| 1019 | |
---|
| 1020 | dimc = rvc._dsize(); |
---|
| 1021 | |
---|
| 1022 | //establish datalinks |
---|
| 1023 | x2g.set_connection(rvxc, rvx.copy_t(-1)); |
---|
| 1024 | cond2g.set_connection(rvxc, rvc); |
---|
| 1025 | |
---|
| 1026 | x2h.set_connection(rvyc, rvx); |
---|
| 1027 | cond2h.set_connection(rvyc, rvc); |
---|
| 1028 | } |
---|
[971] | 1029 | }; |
---|
| 1030 | UIREGISTER(NoiseParticle); |
---|
[811] | 1031 | |
---|
[971] | 1032 | |
---|
[32] | 1033 | } |
---|
[8] | 1034 | #endif // KF_H |
---|
| 1035 | |
---|