[384] | 1 | #include "particles.h" |
---|
[8] | 2 | |
---|
[283] | 3 | namespace bdm { |
---|
[8] | 4 | |
---|
| 5 | using std::endl; |
---|
| 6 | |
---|
[11] | 7 | |
---|
[737] | 8 | void PF::bayes_weights() { |
---|
| 9 | // |
---|
[1064] | 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 | } |
---|
[11] | 15 | |
---|
[1064] | 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; |
---|
[638] | 30 | } |
---|
| 31 | |
---|
[737] | 32 | void PF::bayes ( const vec &yt, const vec &cond ) { |
---|
| 33 | |
---|
[1064] | 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 | } |
---|
[11] | 42 | |
---|
[1064] | 43 | bayes_weights(); |
---|
[11] | 44 | |
---|
[1187] | 45 | if (log_level[logweights]) { |
---|
| 46 | log_level.store( logweights, w); |
---|
| 47 | } |
---|
[1196] | 48 | if (log_level[logNeff]) { |
---|
| 49 | log_level.store( logNeff, 1.0 / ( w * w )); |
---|
| 50 | } |
---|
[1187] | 51 | |
---|
[1064] | 52 | if ( do_resampling() ) { |
---|
| 53 | resample ( ); |
---|
| 54 | } |
---|
[8] | 55 | |
---|
| 56 | } |
---|
| 57 | |
---|
[283] | 58 | // void PF::set_est ( const epdf &epdf0 ) { |
---|
| 59 | // int i; |
---|
[477] | 60 | // |
---|
[283] | 61 | // for ( i=0;i<n;i++ ) { |
---|
| 62 | // _samples ( i ) = epdf0.sample(); |
---|
| 63 | // } |
---|
| 64 | // } |
---|
[8] | 65 | |
---|
[887] | 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 | // } |
---|
[1064] | 75 | // |
---|
[887] | 76 | // vec MPF::mpfepdf::variance() const { |
---|
| 77 | // const vec &w = pf->posterior()._w(); |
---|
[1064] | 78 | // |
---|
[887] | 79 | // vec pom = zeros ( BMs ( 0 )->posterior ().dimension() ); |
---|
| 80 | // vec pom2 = zeros ( BMs ( 0 )->posterior ().dimension() ); |
---|
| 81 | // vec mea; |
---|
[1064] | 82 | // |
---|
[887] | 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 | // } |
---|
[1064] | 92 | // |
---|
[887] | 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 ); |
---|
[1064] | 98 | // |
---|
[887] | 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(); |
---|
[1064] | 110 | // |
---|
[887] | 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 | // } |
---|
[1064] | 127 | // |
---|
| 128 | // |
---|
| 129 | // |
---|
[887] | 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(); |
---|
[1064] | 136 | // |
---|
[887] | 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 | // } |
---|
[1064] | 148 | // |
---|
[887] | 149 | // pf->bayes_weights(); |
---|
[1064] | 150 | // |
---|
[887] | 151 | // ivec ind; |
---|
| 152 | // if ( pf->do_resampling() ) { |
---|
| 153 | // pf->resample ( ind ); |
---|
[1064] | 154 | // |
---|
[887] | 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 | // }; |
---|
[1064] | 164 | // |
---|
| 165 | // |
---|
[887] | 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(); |
---|
[1064] | 172 | // |
---|
[887] | 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(); |
---|
[1064] | 188 | // |
---|
[887] | 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 | // } |
---|
[1064] | 198 | // |
---|
[887] | 199 | // pf->bayes_weights(); |
---|
[1064] | 200 | // |
---|
[887] | 201 | // ivec ind; |
---|
| 202 | // if ( pf->do_resampling() ) { |
---|
| 203 | // pf->resample ( ind ); |
---|
[1064] | 204 | // |
---|
[887] | 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 | // }; |
---|
[1064] | 216 | // |
---|
[887] | 217 | // void MPF_ARXg::from_setting(const libconfig::Setting& set) { |
---|
| 218 | // BM::from_setting(set); |
---|
[1064] | 219 | // |
---|
[887] | 220 | // pf = new PF; |
---|
| 221 | // // prior must be set before BM |
---|
| 222 | // pf->prior_from_set ( set ); |
---|
| 223 | // pf->resmethod_from_set ( set ); |
---|
[1064] | 224 | // |
---|
[887] | 225 | // // read functions g and h |
---|
| 226 | // g=UI::build<fnc>(set,"g"); |
---|
| 227 | // h=UI::build<fnc>(set,"h"); |
---|
[1064] | 228 | // |
---|
[887] | 229 | // set_dim( g->dimension()); |
---|
| 230 | // dimy = h->dimension(); |
---|
[1064] | 231 | // |
---|
[887] | 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 | // } |
---|
[1064] | 241 | // |
---|
[887] | 242 | // yrv = arxo->_yrv(); |
---|
| 243 | // //rvc = arxo->_rvc(); |
---|
[1064] | 244 | // |
---|
[887] | 245 | // validate(); |
---|
[1064] | 246 | // |
---|
[887] | 247 | // } |
---|
[739] | 248 | |
---|
| 249 | |
---|
| 250 | } |
---|
[32] | 251 | //MPF::MPF:{} |
---|