00001
00013 #ifndef PARTICLES_H
00014 #define PARTICLES_H
00015
00016
00017 #include "../stat/exp_family.h"
00018
00019 namespace bdm {
00020
00027 class PF : public BM {
00028 protected:
00030 int n;
00032 eEmp est;
00034 vec &_w;
00036 Array<vec> &_samples;
00038 shared_ptr<mpdf> par;
00040 shared_ptr<mpdf> obs;
00042 vec lls;
00043
00045 RESAMPLING_METHOD resmethod;
00048 double res_threshold;
00049
00052
00054 bool opt_L_smp;
00056 bool opt_L_wei;
00058
00059 public:
00062 PF ( ) : est(), _w ( est._w() ), _samples ( est._samples() ), opt_L_smp ( false ), opt_L_wei ( false ) {
00063 LIDs.set_size ( 5 );
00064 };
00065
00066 void set_parameters (int n0, double res_th0=0.5, RESAMPLING_METHOD rm = SYSTEMATIC ) {
00067 n = n0;
00068 res_threshold = res_th0;
00069 resmethod = rm;
00070 };
00071 void set_model ( shared_ptr<mpdf> par0, shared_ptr<mpdf> obs0) {
00072 par = par0;
00073 obs = obs0;
00074
00075 est.set_rv(par->_rv());
00076 };
00077 void set_statistics ( const vec w0, const epdf &epdf0 ) {
00078 est.set_statistics ( w0, epdf0 );
00079 };
00080 void set_statistics ( const eEmp &epdf0 ) {
00081 bdm_assert_debug(epdf0._rv().equal(par->_rv()),"Incompatibel input");
00082 est=epdf0;
00083 };
00089 void set_options ( const string &opt ) {
00090 BM::set_options ( opt );
00091 opt_L_wei = ( opt.find ( "logweights" ) != string::npos );
00092 opt_L_smp = ( opt.find ( "logsamples" ) != string::npos );
00093 }
00095 virtual void bayes_gensmp();
00097 virtual void bayes_weights();
00099 virtual bool do_resampling(){
00100 double eff = 1.0 / ( _w * _w );
00101 return eff < ( res_threshold*n );
00102 }
00103 void bayes ( const vec &dt );
00105 vec& __w() { return _w; }
00107 vec& _lls() { return lls; }
00109 RESAMPLING_METHOD _resmethod() const { return resmethod; }
00111 const eEmp& posterior() const {return est;}
00112
00125 void from_setting(const Setting &set){
00126 par = UI::build<mpdf>(set,"parameter_pdf",UI::compulsory);
00127 obs = UI::build<mpdf>(set,"observation_pdf",UI::compulsory);
00128
00129 prior_from_set(set);
00130 resmethod_from_set(set);
00131
00132
00133
00134 RV u = par->_rvc().remove_time().subt( par->_rv() );
00135
00136 RV obs_u = obs->_rvc().remove_time().subt( par->_rv() );
00137
00138 u.add(obs_u);
00139
00140 set_drv(concat(obs->_rv(),u) );
00141 }
00143 void resmethod_from_set(const Setting &set){
00144 string resmeth;
00145 if (UI::get(resmeth,set,"resmethod",UI::optional)){
00146 if (resmeth=="systematic") {
00147 resmethod= SYSTEMATIC;
00148 } else {
00149 if (resmeth=="multinomial"){
00150 resmethod=MULTINOMIAL;
00151 } else {
00152 if (resmeth=="stratified"){
00153 resmethod= STRATIFIED;
00154 } else {
00155 bdm_error("Unknown resampling method");
00156 }
00157 }
00158 }
00159 } else {
00160 resmethod=SYSTEMATIC;
00161 };
00162 if(!UI::get(res_threshold, set, "res_threshold", UI::optional)){
00163 res_threshold=0.5;
00164 }
00165 }
00167 void prior_from_set(const Setting & set){
00168 shared_ptr<epdf> pri = UI::build<epdf>(set,"prior",UI::compulsory);
00169
00170 eEmp *test_emp=dynamic_cast<eEmp*>(&(*pri));
00171 if (test_emp) {
00172 est=*test_emp;
00173 } else {
00174 int n;
00175 if (!UI::get(n,set,"n",UI::optional)){n=10;}
00176
00177 set_statistics(ones(n)/n, *pri);
00178 }
00179
00180 }
00181
00182 void validate(){
00183 n=_w.length();
00184 lls=zeros(n);
00185 if (par->_rv()._dsize()>0) {
00186 bdm_assert(par->_rv()._dsize()==est.dimension(),"Mismatch of RV and dimension of posterior" );
00187 }
00188 }
00190 void resample(ivec &ind){
00191 est.resample(ind,resmethod);
00192 }
00194 Array<vec>& __samples(){return _samples;}
00195 };
00196 UIREGISTER(PF);
00197
00205 class MPF : public BM {
00206 protected:
00208 shared_ptr<PF> pf;
00210 Array<BM*> BMs;
00211
00213
00214 class mpfepdf : public epdf {
00216 shared_ptr<PF> &pf;
00218 Array<BM*> &BMs;
00219 public:
00221 mpfepdf (shared_ptr<PF> &pf0, Array<BM*> &BMs0): epdf(), pf(pf0), BMs(BMs0) { };
00223 void read_parameters(){
00224 rv = concat(pf->posterior()._rv(), BMs(0)->posterior()._rv());
00225 dim = pf->posterior().dimension() + BMs(0)->posterior().dimension();
00226 bdm_assert_debug(dim == rv._dsize(), "Wrong name ");
00227 }
00228 vec mean() const {
00229 const vec &w = pf->posterior()._w();
00230 vec pom = zeros ( BMs(0)->posterior ().dimension() );
00231
00232 for ( int i = 0; i < w.length(); i++ ) {
00233 pom += BMs ( i )->posterior().mean() * w ( i );
00234 }
00235 return concat ( pf->posterior().mean(), pom );
00236 }
00237 vec variance() const {
00238 const vec &w = pf->posterior()._w();
00239
00240 vec pom = zeros ( BMs(0)->posterior ().dimension() );
00241 vec pom2 = zeros ( BMs(0)->posterior ().dimension() );
00242 vec mea;
00243
00244 for ( int i = 0; i < w.length(); i++ ) {
00245
00246 mea = BMs ( i )->posterior().mean();
00247 pom += mea * w ( i );
00248
00249 pom2 += ( BMs ( i )->posterior().variance() + pow ( mea, 2 ) ) * w ( i );
00250 }
00251 return concat ( pf->posterior().variance(), pom2 - pow ( pom, 2 ) );
00252 }
00253
00254 void qbounds ( vec &lb, vec &ub, double perc = 0.95 ) const {
00255
00256 vec lbp;
00257 vec ubp;
00258 pf->posterior().qbounds ( lbp, ubp );
00259
00260
00261 int dimC = BMs ( 0 )->posterior().dimension();
00262 int j;
00263
00264 vec lbc ( dimC );
00265 vec ubc ( dimC );
00266
00267 vec Lbc ( dimC );
00268 vec Ubc ( dimC );
00269 Lbc = std::numeric_limits<double>::infinity();
00270 Ubc = -std::numeric_limits<double>::infinity();
00271
00272 for ( int i = 0; i < BMs.length(); i++ ) {
00273
00274 BMs ( i )->posterior().qbounds ( lbc, ubc );
00275
00276 for ( j = 0; j < dimC; j++ ) {
00277 if ( lbc ( j ) < Lbc ( j ) ) {
00278 Lbc ( j ) = lbc ( j );
00279 }
00280 if ( ubc ( j ) > Ubc ( j ) ) {
00281 Ubc ( j ) = ubc ( j );
00282 }
00283 }
00284 }
00285 lb = concat ( lbp, Lbc );
00286 ub = concat ( ubp, Ubc );
00287 }
00288
00289 vec sample() const {
00290 bdm_error ( "Not implemented" );
00291 return vec();
00292 }
00293
00294 double evallog ( const vec &val ) const {
00295 bdm_error ( "not implemented" );
00296 return 0.0;
00297 }
00298 };
00299
00301 mpfepdf jest;
00302
00304 bool opt_L_mea;
00305
00306 public:
00308 MPF () : jest (pf,BMs) {};
00310 void set_parameters ( shared_ptr<mpdf> par0, shared_ptr<mpdf> obs0, int n0, RESAMPLING_METHOD rm = SYSTEMATIC ) {
00311 pf->set_model ( par0, obs0);
00312 pf->set_parameters(n0, rm );
00313 BMs.set_length ( n0 );
00314 }
00316 void set_BM ( const BM &BMcond0 ) {
00317
00318 int n=pf->__w().length();
00319 BMs.set_length(n);
00320
00321
00322 for ( int i = 0; i < n; i++ ) {
00323 BMs ( i ) = BMcond0._copy_();
00324 }
00325 };
00326
00327 void bayes ( const vec &dt );
00328 const epdf& posterior() const {
00329 return jest;
00330 }
00333 void set_options ( const string &opt ) {
00334 BM::set_options(opt);
00335 opt_L_mea = ( opt.find ( "logmeans" ) != string::npos );
00336 }
00337
00339 const BM* _BM ( int i ) {
00340 return BMs ( i );
00341 }
00342
00354 void from_setting(const Setting &set){
00355 shared_ptr<mpdf> par = UI::build<mpdf>(set,"parameter_pdf",UI::compulsory);
00356 shared_ptr<mpdf> obs= new mpdf();
00357
00358 pf = new PF;
00359
00360 pf->prior_from_set(set);
00361 pf->resmethod_from_set(set);
00362 pf->set_model(par,obs);
00363
00364 shared_ptr<BM> BM0 =UI::build<BM>(set,"BM",UI::compulsory);
00365 set_BM(*BM0);
00366
00367 string opt;
00368 if (UI::get(opt,set,"options",UI::optional)){
00369 set_options(opt);
00370 }
00371
00372
00373 RV u = par->_rvc().remove_time().subt( par->_rv() );
00374 set_drv(concat(BM0->_drv(),u) );
00375 validate();
00376 }
00377 void validate(){
00378 try{
00379 pf->validate();
00380 } catch (std::exception &e){
00381 throw UIException("Error in PF part of MPF:");
00382 }
00383 jest.read_parameters();
00384 for ( int i = 0; i < pf->__w().length(); i++ ) {
00385 BMs ( i )->condition ( pf->posterior()._sample ( i ) );
00386 }
00387 }
00388
00389 };
00390 UIREGISTER(MPF);
00391
00392 }
00393 #endif // KF_H
00394
00395
