[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 | |
---|
[1064] | 45 | if ( do_resampling() ) { |
---|
| 46 | resample ( ); |
---|
| 47 | } |
---|
[8] | 48 | |
---|
| 49 | } |
---|
| 50 | |
---|
[283] | 51 | // void PF::set_est ( const epdf &epdf0 ) { |
---|
| 52 | // int i; |
---|
[477] | 53 | // |
---|
[283] | 54 | // for ( i=0;i<n;i++ ) { |
---|
| 55 | // _samples ( i ) = epdf0.sample(); |
---|
| 56 | // } |
---|
| 57 | // } |
---|
[8] | 58 | |
---|
[887] | 59 | // vec MPF::mpfepdf::mean() const { |
---|
| 60 | // const vec &w = pf->posterior()._w(); |
---|
| 61 | // vec pom = zeros ( BMs ( 0 )->posterior ().dimension() ); |
---|
| 62 | // //compute mean of BMs |
---|
| 63 | // for ( int i = 0; i < w.length(); i++ ) { |
---|
| 64 | // pom += BMs ( i )->posterior().mean() * w ( i ); |
---|
| 65 | // } |
---|
| 66 | // return concat ( pf->posterior().mean(), pom ); |
---|
| 67 | // } |
---|
[1064] | 68 | // |
---|
[887] | 69 | // vec MPF::mpfepdf::variance() const { |
---|
| 70 | // const vec &w = pf->posterior()._w(); |
---|
[1064] | 71 | // |
---|
[887] | 72 | // vec pom = zeros ( BMs ( 0 )->posterior ().dimension() ); |
---|
| 73 | // vec pom2 = zeros ( BMs ( 0 )->posterior ().dimension() ); |
---|
| 74 | // vec mea; |
---|
[1064] | 75 | // |
---|
[887] | 76 | // for ( int i = 0; i < w.length(); i++ ) { |
---|
| 77 | // // save current mean |
---|
| 78 | // mea = BMs ( i )->posterior().mean(); |
---|
| 79 | // pom += mea * w ( i ); |
---|
| 80 | // //compute variance |
---|
| 81 | // pom2 += ( BMs ( i )->posterior().variance() + pow ( mea, 2 ) ) * w ( i ); |
---|
| 82 | // } |
---|
| 83 | // return concat ( pf->posterior().variance(), pom2 - pow ( pom, 2 ) ); |
---|
| 84 | // } |
---|
[1064] | 85 | // |
---|
[887] | 86 | // void MPF::mpfepdf::qbounds ( vec &lb, vec &ub, double perc ) const { |
---|
| 87 | // //bounds on particles |
---|
| 88 | // vec lbp; |
---|
| 89 | // vec ubp; |
---|
| 90 | // pf->posterior().qbounds ( lbp, ubp ); |
---|
[1064] | 91 | // |
---|
[887] | 92 | // //bounds on Components |
---|
| 93 | // int dimC = BMs ( 0 )->posterior().dimension(); |
---|
| 94 | // int j; |
---|
| 95 | // // temporary |
---|
| 96 | // vec lbc ( dimC ); |
---|
| 97 | // vec ubc ( dimC ); |
---|
| 98 | // // minima and maxima |
---|
| 99 | // vec Lbc ( dimC ); |
---|
| 100 | // vec Ubc ( dimC ); |
---|
| 101 | // Lbc = std::numeric_limits<double>::infinity(); |
---|
| 102 | // Ubc = -std::numeric_limits<double>::infinity(); |
---|
[1064] | 103 | // |
---|
[887] | 104 | // for ( int i = 0; i < BMs.length(); i++ ) { |
---|
| 105 | // // check Coms |
---|
| 106 | // BMs ( i )->posterior().qbounds ( lbc, ubc ); |
---|
| 107 | // //save either minima or maxima |
---|
| 108 | // for ( j = 0; j < dimC; j++ ) { |
---|
| 109 | // if ( lbc ( j ) < Lbc ( j ) ) { |
---|
| 110 | // Lbc ( j ) = lbc ( j ); |
---|
| 111 | // } |
---|
| 112 | // if ( ubc ( j ) > Ubc ( j ) ) { |
---|
| 113 | // Ubc ( j ) = ubc ( j ); |
---|
| 114 | // } |
---|
| 115 | // } |
---|
| 116 | // } |
---|
| 117 | // lb = concat ( lbp, Lbc ); |
---|
| 118 | // ub = concat ( ubp, Ubc ); |
---|
| 119 | // } |
---|
[1064] | 120 | // |
---|
| 121 | // |
---|
| 122 | // |
---|
[887] | 123 | // void MPF::bayes ( const vec &yt, const vec &cond ) { |
---|
| 124 | // // follows PF::bayes in most places!!! |
---|
| 125 | // int i; |
---|
| 126 | // int n = pf->__w().length(); |
---|
| 127 | // vec &lls = pf->_lls(); |
---|
| 128 | // Array<vec> &samples = pf->__samples(); |
---|
[1064] | 129 | // |
---|
[887] | 130 | // // generate samples - time step |
---|
| 131 | // pf->bayes_gensmp ( vec ( 0 ) ); |
---|
| 132 | // // weight them - data step |
---|
| 133 | // #pragma parallel for |
---|
| 134 | // for ( i = 0; i < n; i++ ) { |
---|
| 135 | // vec bm_cond ( BMs ( i )->dimensionc() ); |
---|
| 136 | // this2bm.fill_cond ( yt, cond, bm_cond ); |
---|
| 137 | // pf2bm.filldown ( samples ( i ), bm_cond ); |
---|
| 138 | // BMs ( i ) -> bayes ( this2bm.pushdown ( yt ), bm_cond ); |
---|
| 139 | // lls ( i ) += BMs ( i )->_ll(); |
---|
| 140 | // } |
---|
[1064] | 141 | // |
---|
[887] | 142 | // pf->bayes_weights(); |
---|
[1064] | 143 | // |
---|
[887] | 144 | // ivec ind; |
---|
| 145 | // if ( pf->do_resampling() ) { |
---|
| 146 | // pf->resample ( ind ); |
---|
[1064] | 147 | // |
---|
[887] | 148 | // #pragma omp parallel for |
---|
| 149 | // for ( i = 0; i < n; i++ ) { |
---|
| 150 | // if ( ind ( i ) != i ) {//replace the current Bm by a new one |
---|
| 151 | // delete BMs ( i ); |
---|
| 152 | // BMs ( i ) = (BM*) BMs ( ind ( i ) )->_copy(); //copy constructor |
---|
| 153 | // } |
---|
| 154 | // }; |
---|
| 155 | // } |
---|
| 156 | // }; |
---|
[1064] | 157 | // |
---|
| 158 | // |
---|
[887] | 159 | // void MPF_ARXg::bayes ( const vec &yt, const vec &cond ) { |
---|
| 160 | // // follows PF::bayes in most places!!! |
---|
| 161 | // int i; |
---|
| 162 | // int n = pf->__w().length(); |
---|
| 163 | // vec &lls = pf->_lls(); |
---|
| 164 | // Array<vec> &samples = pf->__samples(); |
---|
[1064] | 165 | // |
---|
[887] | 166 | // // generate samples - time step |
---|
| 167 | // for (int i =0;i<n; i++){ |
---|
| 168 | // mat M; chmat R; |
---|
| 169 | // BMsp(i)->posterior().sample_mat(M,R); |
---|
| 170 | // vec diff=R._Ch().T()*randn(samples(i).length()); |
---|
| 171 | // samples(i) = g->eval(samples(i)) + diff; |
---|
| 172 | // //////////// |
---|
| 173 | // // samples(i) = cond; |
---|
| 174 | // ///////// |
---|
| 175 | // BMsp(i)->bayes(diff); |
---|
| 176 | // } |
---|
| 177 | // // weight them - data step |
---|
| 178 | // enorm<ldmat> ooo; |
---|
| 179 | // ooo.set_parameters(zeros(2),0.1*eye(2)); |
---|
| 180 | // ooo.validate(); |
---|
[1064] | 181 | // |
---|
[887] | 182 | // #pragma parallel for |
---|
| 183 | // for ( i = 0; i < n; i++ ) { |
---|
| 184 | // vec tmp=yt - h->eval(samples(i)); |
---|
| 185 | // BMso ( i ) -> bayes ( tmp ); |
---|
| 186 | // lls ( i ) += BMso ( i )->_ll(); |
---|
| 187 | // ///// |
---|
| 188 | // ooo._mu()=h->eval(samples(i)); |
---|
| 189 | // lls(i) = ooo.evallog_nn(yt); |
---|
| 190 | // } |
---|
[1064] | 191 | // |
---|
[887] | 192 | // pf->bayes_weights(); |
---|
[1064] | 193 | // |
---|
[887] | 194 | // ivec ind; |
---|
| 195 | // if ( pf->do_resampling() ) { |
---|
| 196 | // pf->resample ( ind ); |
---|
[1064] | 197 | // |
---|
[887] | 198 | // #pragma omp parallel for |
---|
| 199 | // for ( i = 0; i < n; i++ ) { |
---|
| 200 | // if ( ind ( i ) != i ) {//replace the current Bm by a new one |
---|
| 201 | // delete BMso ( i ); |
---|
| 202 | // BMso ( i ) = (ARX*) BMso ( ind ( i ) )->_copy(); //copy constructor |
---|
| 203 | // delete BMsp ( i ); |
---|
| 204 | // BMsp ( i ) = (ARX*) BMsp ( ind ( i ) )->_copy(); //copy constructor |
---|
| 205 | // } |
---|
| 206 | // }; |
---|
| 207 | // } |
---|
| 208 | // }; |
---|
[1064] | 209 | // |
---|
[887] | 210 | // void MPF_ARXg::from_setting(const libconfig::Setting& set) { |
---|
| 211 | // BM::from_setting(set); |
---|
[1064] | 212 | // |
---|
[887] | 213 | // pf = new PF; |
---|
| 214 | // // prior must be set before BM |
---|
| 215 | // pf->prior_from_set ( set ); |
---|
| 216 | // pf->resmethod_from_set ( set ); |
---|
[1064] | 217 | // |
---|
[887] | 218 | // // read functions g and h |
---|
| 219 | // g=UI::build<fnc>(set,"g"); |
---|
| 220 | // h=UI::build<fnc>(set,"h"); |
---|
[1064] | 221 | // |
---|
[887] | 222 | // set_dim( g->dimension()); |
---|
| 223 | // dimy = h->dimension(); |
---|
[1064] | 224 | // |
---|
[887] | 225 | // shared_ptr<ARX> arxo=UI::build<ARX>(set,"arxo"); |
---|
| 226 | // shared_ptr<ARX> arxp=UI::build<ARX>(set,"arxp"); |
---|
| 227 | // int n = pf->__samples().length(); |
---|
| 228 | // BMso.set_length(n); |
---|
| 229 | // BMsp.set_length(n); |
---|
| 230 | // for(int i=0; i<n; i++){ |
---|
| 231 | // BMso(i)=arxo->_copy(); |
---|
| 232 | // BMsp(i)=arxp->_copy(); |
---|
| 233 | // } |
---|
[1064] | 234 | // |
---|
[887] | 235 | // yrv = arxo->_yrv(); |
---|
| 236 | // //rvc = arxo->_rvc(); |
---|
[1064] | 237 | // |
---|
[887] | 238 | // validate(); |
---|
[1064] | 239 | // |
---|
[887] | 240 | // } |
---|
[739] | 241 | |
---|
| 242 | |
---|
| 243 | } |
---|
[32] | 244 | //MPF::MPF:{} |
---|