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