Changeset 887 for library/bdm

Timestamp:

03/29/10 23:02:03 (14 years ago)

Author:

smidl

Message:

new base for particle filtering

Location:

library/bdm

Files:

• estim/particles.cpp (modified) (3 diffs)
• estim/particles.h (modified) (8 diffs)
• stat/exp_family.cpp (modified) (4 diffs)
• stat/exp_family.h (modified) (5 diffs)

library/bdm/estim/particles.cpp

@@ -5 +5 @@
 using std::endl;
 
-void PF::bayes_gensmp ( const vec &ut ) {
-        if ( ut.length() > 0 ) {
-                vec cond ( par->dimensionc() );
-                cond.set_subvector ( par->dimension(), ut );
-#pragma parallel for
-                for ( int i = 0; i < n; i++ ) {
-                        cond.set_subvector ( 0, _samples ( i ) );
-                        _samples ( i ) = par->samplecond ( cond );
-                        lls ( i ) = 0;
-                }
-        } else {
-#pragma parallel for
-                for ( int i = 0; i < n; i++ ) {
-                        _samples ( i ) = par->samplecond ( _samples ( i ) );
-                        lls ( i ) = 0;
-                }
-        }
-}
 
 void PF::bayes_weights() {
@@ -29 +11 @@
         // compute weights
         for ( int  i = 0; i < n; i++ ) {
-                _w ( i ) *= exp ( lls ( i ) - mlls ); // multiply w by likelihood
+                w ( i ) *= exp ( lls ( i ) - mlls ); // multiply w by likelihood
         }
 
         //renormalize
-        double sw = sum ( _w );
+        double sw = sum ( w );
         if ( !std::isfinite ( sw ) ) {
                 for ( int i = 0; i < n; i++ ) {
-                        if ( !std::isfinite ( _w ( i ) ) ) {
-                                _w ( i ) = 0;
+                        if ( !std::isfinite ( w ( i ) ) ) {
+                                w ( i ) = 0;
                         }
                 }
-                sw = sum ( _w );
+                sw = sum ( w );
                 if ( !std::isfinite ( sw ) || sw == 0.0 ) {
                         bdm_error ( "Particle filter is lost; no particle is good enough." );
                 }
         }
-        _w /= sw;
+        w /= sw;
 }
 
 void PF::bayes ( const vec &yt, const vec &cond ) {
-        const vec &ut = cond; //todo
-
-        int i;
-        // generate samples - time step
-        bayes_gensmp ( ut );
+
+/*      if (auxiliary){
+                ...
+        }*/
         // weight them - data step
-        for ( i = 0; i < n; i++ ) {
-                lls ( i ) += obs->evallogcond ( yt, _samples ( i ) ); //+= because lls may have something from gensmp!
+        for (int i = 0; i < n; i++ ) {
+                particles(i)->bayes(yt,cond);
+                lls ( i ) = particles(i)->_ll(); //+= because lls may have something from gensmp!
         }
 
         bayes_weights();
 
-        if ( do_resampling() ) {
-                est.resample ( resmethod );
-        }
+        if ( do_resampling() ) { resample ( ); }
 
 }
@@ -75 +55 @@
 // }
 
-vec MPF::mpfepdf::mean() const {
-        const vec &w = pf->posterior()._w();
-        vec pom = zeros ( BMs ( 0 )->posterior ().dimension() );
-        //compute mean of BMs
-        for ( int i = 0; i < w.length(); i++ ) {
-                pom += BMs ( i )->posterior().mean() * w ( i );
-        }
-        return concat ( pf->posterior().mean(), pom );
-}
-
-vec MPF::mpfepdf::variance() const {
-        const vec &w = pf->posterior()._w();
-
-        vec pom = zeros ( BMs ( 0 )->posterior ().dimension() );
-        vec pom2 = zeros ( BMs ( 0 )->posterior ().dimension() );
-        vec mea;
-
-        for ( int i = 0; i < w.length(); i++ ) {
-                // save current mean
-                mea = BMs ( i )->posterior().mean();
-                pom += mea * w ( i );
-                //compute variance
-                pom2 += ( BMs ( i )->posterior().variance() + pow ( mea, 2 ) ) * w ( i );
-        }
-        return concat ( pf->posterior().variance(), pom2 - pow ( pom, 2 ) );
-}
-
-void MPF::mpfepdf::qbounds ( vec &lb, vec &ub, double perc ) const {
-        //bounds on particles
-        vec lbp;
-        vec ubp;
-        pf->posterior().qbounds ( lbp, ubp );
-
-        //bounds on Components
-        int dimC = BMs ( 0 )->posterior().dimension();
-        int j;
-        // temporary
-        vec lbc ( dimC );
-        vec ubc ( dimC );
-        // minima and maxima
-        vec Lbc ( dimC );
-        vec Ubc ( dimC );
-        Lbc = std::numeric_limits<double>::infinity();
-        Ubc = -std::numeric_limits<double>::infinity();
-
-        for ( int i = 0; i < BMs.length(); i++ ) {
-                // check Coms
-                BMs ( i )->posterior().qbounds ( lbc, ubc );
-                //save either minima or maxima
-                for ( j = 0; j < dimC; j++ ) {
-                        if ( lbc ( j ) < Lbc ( j ) ) {
-                                Lbc ( j ) = lbc ( j );
-                        }
-                        if ( ubc ( j ) > Ubc ( j ) ) {
-                                Ubc ( j ) = ubc ( j );
-                        }
-                }
-        }
-        lb = concat ( lbp, Lbc );
-        ub = concat ( ubp, Ubc );
-}
-
-
-
-void MPF::bayes ( const vec &yt, const vec &cond ) {
-        // follows PF::bayes in most places!!!
-        int i;
-        int n = pf->__w().length();
-        vec &lls = pf->_lls();
-        Array<vec> &samples = pf->__samples();
-
-        // generate samples - time step
-        pf->bayes_gensmp ( vec ( 0 ) );
-        // weight them - data step
-#pragma parallel for
-        for ( i = 0; i < n; i++ ) {
-                vec bm_cond ( BMs ( i )->dimensionc() );
-                this2bm.fill_cond ( yt, cond, bm_cond );
-                pf2bm.filldown ( samples ( i ), bm_cond );
-                BMs ( i ) -> bayes ( this2bm.pushdown ( yt ), bm_cond );
-                lls ( i ) += BMs ( i )->_ll();
-        }
-
-        pf->bayes_weights();
-
-        ivec ind;
-        if ( pf->do_resampling() ) {
-                pf->resample ( ind );
-
-#pragma omp parallel for
-                for ( i = 0; i < n; i++ ) {
-                        if ( ind ( i ) != i ) {//replace the current Bm by a new one
-                                delete BMs ( i );
-                                BMs ( i ) = (BM*) BMs ( ind ( i ) )->_copy(); //copy constructor
-                        }
-                };
-        }
-};
-
-
-void MPF_ARXg::bayes ( const vec &yt, const vec &cond ) {
-        // follows PF::bayes in most places!!!
-        int i;
-        int n = pf->__w().length();
-        vec &lls = pf->_lls();
-        Array<vec> &samples = pf->__samples();
-        
-        // generate samples - time step
-        for (int i =0;i<n; i++){
-                mat M; chmat R;
-                BMsp(i)->posterior().sample_mat(M,R);
-                vec diff=R._Ch().T()*randn(samples(i).length());
-                samples(i) = g->eval(samples(i)) + diff;
-                ////////////
-//              samples(i) = cond;
-                /////////
-                BMsp(i)->bayes(diff);
-        }
-        // weight them - data step
-        enorm<ldmat> ooo;
-        ooo.set_parameters(zeros(2),0.1*eye(2));
-        ooo.validate();
-        
-        #pragma parallel for
-        for ( i = 0; i < n; i++ ) {
-                vec tmp=yt - h->eval(samples(i));
-                BMso ( i ) -> bayes ( tmp );
-                lls ( i ) += BMso ( i )->_ll();
-                /////
-                ooo._mu()=h->eval(samples(i));
-                lls(i) = ooo.evallog_nn(yt);
-        }
-        
-        pf->bayes_weights();
-        
-        ivec ind;
-        if ( pf->do_resampling() ) {
-                pf->resample ( ind );
-                
-                #pragma omp parallel for
-                for ( i = 0; i < n; i++ ) {
-                        if ( ind ( i ) != i ) {//replace the current Bm by a new one
-                                delete BMso ( i );
-                                BMso ( i ) = (ARX*) BMso ( ind ( i ) )->_copy(); //copy constructor
-                                delete BMsp ( i );
-                                BMsp ( i ) = (ARX*) BMsp ( ind ( i ) )->_copy(); //copy constructor
-                        }
-                };
-        }
-};
-
-void MPF_ARXg::from_setting(const libconfig::Setting& set) {
-        BM::from_setting(set);
-        
-        pf = new PF;
-        // prior must be set before BM
-        pf->prior_from_set ( set );
-        pf->resmethod_from_set ( set );
-                
-        // read functions g and h
-        g=UI::build<fnc>(set,"g");
-        h=UI::build<fnc>(set,"h");
-        
-        set_dim( g->dimension());
-        dimy = h->dimension();
-        
-        shared_ptr<ARX> arxo=UI::build<ARX>(set,"arxo");
-        shared_ptr<ARX> arxp=UI::build<ARX>(set,"arxp");
-        int n = pf->__samples().length();
-        BMso.set_length(n);
-        BMsp.set_length(n);
-        for(int i=0; i<n; i++){
-                BMso(i)=arxo->_copy();
-                BMsp(i)=arxp->_copy();
-        }
-        
-        yrv = arxo->_yrv();
-        //rvc = arxo->_rvc();
-        
-        validate();
-        
-}
+// vec MPF::mpfepdf::mean() const {
+//      const vec &w = pf->posterior()._w();
+//      vec pom = zeros ( BMs ( 0 )->posterior ().dimension() );
+//      //compute mean of BMs
+//      for ( int i = 0; i < w.length(); i++ ) {
+//              pom += BMs ( i )->posterior().mean() * w ( i );
+//      }
+//      return concat ( pf->posterior().mean(), pom );
+// }
+// 
+// vec MPF::mpfepdf::variance() const {
+//      const vec &w = pf->posterior()._w();
+// 
+//      vec pom = zeros ( BMs ( 0 )->posterior ().dimension() );
+//      vec pom2 = zeros ( BMs ( 0 )->posterior ().dimension() );
+//      vec mea;
+// 
+//      for ( int i = 0; i < w.length(); i++ ) {
+//              // save current mean
+//              mea = BMs ( i )->posterior().mean();
+//              pom += mea * w ( i );
+//              //compute variance
+//              pom2 += ( BMs ( i )->posterior().variance() + pow ( mea, 2 ) ) * w ( i );
+//      }
+//      return concat ( pf->posterior().variance(), pom2 - pow ( pom, 2 ) );
+// }
+// 
+// void MPF::mpfepdf::qbounds ( vec &lb, vec &ub, double perc ) const {
+//      //bounds on particles
+//      vec lbp;
+//      vec ubp;
+//      pf->posterior().qbounds ( lbp, ubp );
+// 
+//      //bounds on Components
+//      int dimC = BMs ( 0 )->posterior().dimension();
+//      int j;
+//      // temporary
+//      vec lbc ( dimC );
+//      vec ubc ( dimC );
+//      // minima and maxima
+//      vec Lbc ( dimC );
+//      vec Ubc ( dimC );
+//      Lbc = std::numeric_limits<double>::infinity();
+//      Ubc = -std::numeric_limits<double>::infinity();
+// 
+//      for ( int i = 0; i < BMs.length(); i++ ) {
+//              // check Coms
+//              BMs ( i )->posterior().qbounds ( lbc, ubc );
+//              //save either minima or maxima
+//              for ( j = 0; j < dimC; j++ ) {
+//                      if ( lbc ( j ) < Lbc ( j ) ) {
+//                              Lbc ( j ) = lbc ( j );
+//                      }
+//                      if ( ubc ( j ) > Ubc ( j ) ) {
+//                              Ubc ( j ) = ubc ( j );
+//                      }
+//              }
+//      }
+//      lb = concat ( lbp, Lbc );
+//      ub = concat ( ubp, Ubc );
+// }
+// 
+// 
+// 
+// void MPF::bayes ( const vec &yt, const vec &cond ) {
+//      // follows PF::bayes in most places!!!
+//      int i;
+//      int n = pf->__w().length();
+//      vec &lls = pf->_lls();
+//      Array<vec> &samples = pf->__samples();
+// 
+//      // generate samples - time step
+//      pf->bayes_gensmp ( vec ( 0 ) );
+//      // weight them - data step
+// #pragma parallel for
+//      for ( i = 0; i < n; i++ ) {
+//              vec bm_cond ( BMs ( i )->dimensionc() );
+//              this2bm.fill_cond ( yt, cond, bm_cond );
+//              pf2bm.filldown ( samples ( i ), bm_cond );
+//              BMs ( i ) -> bayes ( this2bm.pushdown ( yt ), bm_cond );
+//              lls ( i ) += BMs ( i )->_ll();
+//      }
+// 
+//      pf->bayes_weights();
+// 
+//      ivec ind;
+//      if ( pf->do_resampling() ) {
+//              pf->resample ( ind );
+// 
+// #pragma omp parallel for
+//              for ( i = 0; i < n; i++ ) {
+//                      if ( ind ( i ) != i ) {//replace the current Bm by a new one
+//                              delete BMs ( i );
+//                              BMs ( i ) = (BM*) BMs ( ind ( i ) )->_copy(); //copy constructor
+//                      }
+//              };
+//      }
+// };
+// 
+// 
+// void MPF_ARXg::bayes ( const vec &yt, const vec &cond ) {
+//      // follows PF::bayes in most places!!!
+//      int i;
+//      int n = pf->__w().length();
+//      vec &lls = pf->_lls();
+//      Array<vec> &samples = pf->__samples();
+//      
+//      // generate samples - time step
+//      for (int i =0;i<n; i++){
+//              mat M; chmat R;
+//              BMsp(i)->posterior().sample_mat(M,R);
+//              vec diff=R._Ch().T()*randn(samples(i).length());
+//              samples(i) = g->eval(samples(i)) + diff;
+//              ////////////
+// //           samples(i) = cond;
+//              /////////
+//              BMsp(i)->bayes(diff);
+//      }
+//      // weight them - data step
+//      enorm<ldmat> ooo;
+//      ooo.set_parameters(zeros(2),0.1*eye(2));
+//      ooo.validate();
+//      
+//      #pragma parallel for
+//      for ( i = 0; i < n; i++ ) {
+//              vec tmp=yt - h->eval(samples(i));
+//              BMso ( i ) -> bayes ( tmp );
+//              lls ( i ) += BMso ( i )->_ll();
+//              /////
+//              ooo._mu()=h->eval(samples(i));
+//              lls(i) = ooo.evallog_nn(yt);
+//      }
+//      
+//      pf->bayes_weights();
+//      
+//      ivec ind;
+//      if ( pf->do_resampling() ) {
+//              pf->resample ( ind );
+//              
+//              #pragma omp parallel for
+//              for ( i = 0; i < n; i++ ) {
+//                      if ( ind ( i ) != i ) {//replace the current Bm by a new one
+//                              delete BMso ( i );
+//                              BMso ( i ) = (ARX*) BMso ( ind ( i ) )->_copy(); //copy constructor
+//                              delete BMsp ( i );
+//                              BMsp ( i ) = (ARX*) BMsp ( ind ( i ) )->_copy(); //copy constructor
+//                      }
+//              };
+//      }
+// };
+// 
+// void MPF_ARXg::from_setting(const libconfig::Setting& set) {
+//      BM::from_setting(set);
+//      
+//      pf = new PF;
+//      // prior must be set before BM
+//      pf->prior_from_set ( set );
+//      pf->resmethod_from_set ( set );
+//              
+//      // read functions g and h
+//      g=UI::build<fnc>(set,"g");
+//      h=UI::build<fnc>(set,"h");
+//      
+//      set_dim( g->dimension());
+//      dimy = h->dimension();
+//      
+//      shared_ptr<ARX> arxo=UI::build<ARX>(set,"arxo");
+//      shared_ptr<ARX> arxp=UI::build<ARX>(set,"arxp");
+//      int n = pf->__samples().length();
+//      BMso.set_length(n);
+//      BMsp.set_length(n);
+//      for(int i=0; i<n; i++){
+//              BMso(i)=arxo->_copy();
+//              BMsp(i)=arxp->_copy();
+//      }
+//      
+//      yrv = arxo->_yrv();
+//      //rvc = arxo->_rvc();
+//      
+//      validate();
+//      
+// }
 
 

library/bdm/estim/particles.h

@@ -16 +16 @@
 
 #include "../estim/arx_ext.h"
+#include "../stat/emix.h"
 
 namespace bdm {
+        
+//! class used in PF
+class MarginalizedParticle : public BM{
+        protected:
+        //! discrte particle
+        dirac est_emp;
+        //! internal Bayes Model
+        shared_ptr<BM> bm; 
+        //! pdf with for transitional par
+        shared_ptr<pdf> par; // pdf for non-linear part
+        //! link from this to bm
+        shared_ptr<datalink_part> cond2bm;
+        //! link from cond to par
+        shared_ptr<datalink_part> cond2par;
+        //! link from emp 2 par
+        shared_ptr<datalink_part> emp2bm;
+        //! link from emp 2 par
+        shared_ptr<datalink_part> emp2par;
+        
+        //! custom posterior - product
+        class eprod_2:public eprod_base {
+                protected:
+                MarginalizedParticle &mp;
+                public:
+                eprod_2(MarginalizedParticle &m):mp(m){}
+                const epdf* factor(int i) const {return (i==0) ? &mp.bm->posterior() : &mp.est_emp;}
+                const int no_factors() const {return 2;}
+        } est;
+        
+        public:
+                MarginalizedParticle():est(*this){};
+                BM* _copy() const{return new MarginalizedParticle(*this);};
+                void bayes(const vec &dt, const vec &cond){
+                        vec par_cond(par->dimensionc());
+                        cond2par->filldown(cond,par_cond); // copy ut
+                        emp2par->filldown(est_emp._point(),par_cond); // copy xt-1
+                        
+                        //sample new particle
+                        est_emp.set_point(par->samplecond(par_cond));
+                        //if (evalll)
+                        vec bm_cond(bm->dimensionc());
+                        cond2bm->filldown(cond, bm_cond);// set e.g. ut
+                        emp2bm->filldown(est_emp._point(), bm_cond);// set e.g. ut
+                        bm->bayes(dt, bm_cond);
+                        ll=bm->_ll();
+                }
+                const eprod_2& posterior() const {return est;}
+        
+                void set_prior(const epdf *pdf0){
+                        const eprod *ep=dynamic_cast<const eprod*>(pdf0);
+                        if (ep){ // full prior
+                                bdm_assert(ep->no_factors()==2,"Incompatible prod");
+                                bm->set_prior(ep->factor(0));
+                                est_emp.set_point(ep->factor(1)->sample());
+                        } else {
+                                // assume prior is only for emp;
+                                est_emp.set_point(pdf0->sample());
+                        }
+                }
+
+                /*! parse structure
+                \code
+                class = "BootstrapParticle";
+                parameter_pdf = {class = 'epdf_offspring', ...};
+                observation_pdf = {class = 'epdf_offspring',...};
+                \endcode
+                If rvs are set, then it checks for compatibility.
+                */
+                void from_setting(const Setting &set){
+                        BM::from_setting ( set );
+                        par = UI::build<pdf> ( set, "parameter_pdf", UI::compulsory );
+                        bm = UI::build<BM> ( set, "bm", UI::compulsory );
+                }
+                void validate(){
+                        yrv = bm->_rv();
+                        dimy = bm->dimension();
+                        set_rv( par->_rv());
+                        set_dim( par->dimension());
+
+                        rvc = par->_rvc().subt(par->_rv().copy_t(-1));
+                        rvc.add(bm->_rvc()); //
+                        
+                        cond2bm=new datalink_part;
+                        cond2par=new datalink_part;
+                        emp2bm  =new datalink_part;
+                        emp2par =new datalink_part;
+                        cond2bm->set_connection(bm->_rvc(), rvc);
+                        cond2par->set_connection(par->_rvc(), rvc);
+                        emp2bm->set_connection(bm->_rvc(), _rv());
+                        emp2par->set_connection(par->_rvc(), _rv().copy_t(-1));
+                        
+                        dimc = rvc._dsize();
+                };
+};
+UIREGISTER(MarginalizedParticle);
+
+//! class used in PF
+class BootstrapParticle : public BM{
+        dirac est;
+        shared_ptr<pdf> par;
+        shared_ptr<pdf> obs;
+        shared_ptr<datalink_part> cond2par;
+        shared_ptr<datalink_part> cond2obs;
+        shared_ptr<datalink_part> xt2obs;
+        shared_ptr<datalink_part> xtm2par;
+        public:
+                BM* _copy() const{return new BootstrapParticle(*this);};
+                void bayes(const vec &dt, const vec &cond){
+                        vec par_cond(par->dimensionc());
+                        cond2par->filldown(cond,par_cond); // copy ut
+                        xtm2par->filldown(est._point(),par_cond); // copy xt-1
+                        
+                        //sample new particle
+                        est.set_point(par->samplecond(par_cond));
+                        //if (evalll)
+                        vec obs_cond(obs->dimensionc());
+                        cond2obs->filldown(cond, obs_cond);// set e.g. ut
+                        xt2obs->filldown(est._point(), obs_cond);// set e.g. ut
+                        ll=obs->evallogcond(dt,obs_cond);
+                }
+                const dirac& posterior() const {return est;}
+                
+                void set_prior(const epdf *pdf0){est.set_point(pdf0->sample());}
+                
+                /*! parse structure
+                \code
+                class = "BootstrapParticle";
+                parameter_pdf = {class = 'epdf_offspring', ...};
+                observation_pdf = {class = 'epdf_offspring',...};
+                \endcode
+                If rvs are set, then it checks for compatibility.
+                */
+                void from_setting(const Setting &set){
+                        BM::from_setting ( set );
+                        par = UI::build<pdf> ( set, "parameter_pdf", UI::compulsory );
+                        obs = UI::build<pdf> ( set, "observation_pdf", UI::compulsory );
+                }
+                void validate(){
+                        yrv = obs->_rv();
+                        dimy = obs->dimension();
+                        set_rv( par->_rv());
+                        set_dim( par->dimension());
+                        
+                        rvc = par->_rvc().subt(par->_rv().copy_t(-1));
+                        rvc.add(obs->_rvc()); //
+                        
+                        cond2obs=new datalink_part;
+                        cond2par=new datalink_part;
+                        xt2obs  =new datalink_part;
+                        xtm2par =new datalink_part;
+                        cond2obs->set_connection(obs->_rvc(), rvc);
+                        cond2par->set_connection(par->_rvc(), rvc);
+                        xt2obs->set_connection(obs->_rvc(), _rv());
+                        xtm2par->set_connection(par->_rvc(), _rv().copy_t(-1));
+                        
+                        dimc = rvc._dsize();
+                };
+};
+UIREGISTER(BootstrapParticle);
+
 
 /*!
@@ -33 +194 @@
         LOG_LEVEL(PF,weights,samples);
         
+        class pf_mix: public emix_base{
+                Array<BM*> &bms;
+                public:
+                        pf_mix(vec &w0, Array<BM*> &bms0):emix_base(w0),bms(bms0){}
+                        const epdf* component(const int &i)const{return &(bms(i)->posterior());}
+                        int no_coms() const {return bms.length();}
+        };
 protected:
         //!number of particles;
         int n;
         //!posterior density
-        eEmp est;
-        //! pointer into \c eEmp
-        vec &_w;
-        //! pointer into \c eEmp
-        Array<vec> &_samples;
-        //! Parameter evolution model
-        shared_ptr<pdf> par;
-        //! Observation model
-        shared_ptr<pdf> obs;
+        pf_mix est;
+        //! weights;
+        vec w;
+        //! particles
+        Array<BM*> particles;
         //! internal structure storing loglikelihood of predictions
         vec lls;
@@ -62 +226 @@
         //! \name Constructors
         //!@{
-        PF ( ) : est(), _w ( est._w() ), _samples ( est._samples() ) {
-        };
+        PF ( ) : est(w,particles) { };
 
         void set_parameters ( int n0, double res_th0 = 0.5, RESAMPLING_METHOD rm = SYSTEMATIC ) {
@@ -70 +233 @@
                 resmethod = rm;
         };
-        void set_model ( shared_ptr<pdf> par0, shared_ptr<pdf> obs0 ) {
-                par = par0;
-                obs = obs0;
+        void set_model ( const BM *particle0, const epdf *prior) {
+                if (n>0){
+                        particles.set_length(n);
+                        for (int i=0; i<n;i++){
+                                particles(i) = particle0->_copy();
+                                particles(i)->set_prior(prior);
+                        }
+                }
                 // set values for posterior
-                est.set_rv ( par->_rv() );
+                est.set_rv ( particle0->posterior()._rv() );
         };
         void set_statistics ( const vec w0, const epdf &epdf0 ) {
-                est.set_statistics ( w0, epdf0 );
+                //est.set_statistics ( w0, epdf0 );
         };
-        void set_statistics ( const eEmp &epdf0 ) {
+/*      void set_statistics ( const eEmp &epdf0 ) {
                 bdm_assert_debug ( epdf0._rv().equal ( par->_rv() ), "Incompatible input" );
                 est = epdf0;
-        };
+        };*/
         //!@}
 
-        //! Set posterior density by sampling from epdf0
-        //! bayes I - generate samples and add their weights to lls
-        virtual void bayes_gensmp ( const vec &ut );
-        //! bayes II - compute weights of the
+        //! bayes compute weights of the
         virtual void bayes_weights();
         //! important part of particle filtering - decide if it is time to perform resampling
         virtual bool do_resampling() {
-                double eff = 1.0 / ( _w * _w );
+                double eff = 1.0 / ( w * w );
                 return eff < ( res_threshold*n );
         }
         void bayes ( const vec &yt, const vec &cond );
-        //!access function
-        vec& __w() {
-                return _w;
-        }
         //!access function
         vec& _lls() {
@@ -109 +270 @@
         }
         //! return correctly typed posterior (covariant return)
-        const eEmp& posterior() const {
+        const pf_mix& posterior() const {
                 return est;
         }
@@ -127 +288 @@
         void from_setting ( const Setting &set ) {
                 BM::from_setting ( set );
-                par = UI::build<pdf> ( set, "parameter_pdf", UI::compulsory );
-                obs = UI::build<pdf> ( set, "observation_pdf", UI::compulsory );
-
-                prior_from_set ( set );
+                
+                shared_ptr<BM> bm0 = UI::build<BM>(set, "particle",UI::compulsory);
+                
+                shared_ptr<epdf> pri = UI::build<epdf> ( set, "prior", UI::compulsory );
+                n =0;
+                UI::get(n,set,"n",UI::optional);;
+                if (n>0){
+                        particles.set_length(n);
+                        for(int i=0;i<n;i++){particles(i)=bm0->_copy();}
+                        w = ones(n)/n;
+                }
+                set_prior(pri.get());
+                // set resampling method
                 resmethod_from_set ( set );
-                // set resampling method
                 //set drv
-                //find potential input - what remains in rvc when we subtract rv
-                RV u = par->_rvc().remove_time().subt ( par->_rv() );
-                //find potential input - what remains in rvc when we subtract x_t
-                RV obs_u = obs->_rvc().remove_time().subt ( par->_rv() );
-
-                u.add ( obs_u ); // join both u, and check if they do not overlap
-
-                set_yrv ( obs->_rv() );
-                rvc = u;
+
+                set_yrv ( bm0->_rv() );
+                rvc = bm0->_rvc();
+                BM::set_rv(bm0->_rv());
+                yrv=bm0->_yrv();
+        }
+        
+        void set_prior(const epdf *pri){
+                const emix_base *emi=dynamic_cast<const emix_base*>(pri);
+                if (emi) {
+                        bdm_assert(particles.length()>0, "initial particle is not assigned");
+                        n = emi->_w().length();
+                        int old_n = particles.length();
+                        if (n!=old_n){
+                                particles.set_length(n,true);
+                        } 
+                        for(int i=old_n;i<n;i++){particles(i)=particles(0)->_copy();}
+                        
+                        for (int i =0; i<n; i++){
+                                particles(i)->set_prior(emi->_com(i));
+                        }
+                } else {
+                        // try to find "n"
+                        bdm_assert(n>0, "Field 'n' must be filled when prior is not of type emix");
+                        for (int i =0; i<n; i++){
+                                particles(i)->set_prior(pri);
+                        }
+                        
+                }
         }
         //! auxiliary function reading parameter 'resmethod' from configuration file
@@ -165 +354 @@
                 };
                 if ( !UI::get ( res_threshold, set, "res_threshold", UI::optional ) ) {
-                        res_threshold = 0.5;
+                        res_threshold = 0.9;
                 }
                 //validate();
-        }
-        //! load prior information from set and set internal structures accordingly
-        void prior_from_set ( const Setting & set ) {
-                shared_ptr<epdf> pri = UI::build<epdf> ( set, "prior", UI::compulsory );
-
-                eEmp *test_emp = dynamic_cast<eEmp*> ( & ( *pri ) );
-                if ( test_emp ) { // given pdf is sampled
-                        est = *test_emp;
-                } else {
-                        //int n;
-                        if ( !UI::get ( n, set, "n", UI::optional ) ) {
-                                n = 10;
-                        }
-                        // sample from prior
-                        set_statistics ( ones ( n ) / n, *pri );
-                }
-                n = est._w().length();
-                lls=zeros(n);
         }
 
         void validate() {
                 BM::validate();
-                bdm_assert ( par, "PF::parameter_pdf not specified." );
-                n = _w.length();
+                est.validate();
+                bdm_assert ( n>0, "empty particle pool" );
+                n = w.length();
                 lls = zeros ( n );
 
-                if ( par->_rv()._dsize() > 0 ) {
-                        bdm_assert ( par->_rv()._dsize() == est.dimension(), "Mismatch of RV and dimension of posterior" );
+                if ( particles(0)->_rv()._dsize() > 0 ) {
+                        bdm_assert (  particles(0)->_rv()._dsize() == est.dimension(), "Mismatch of RV and dimension of posterior" );
                 }
         }
         //! resample posterior density (from outside - see MPF)
-        void resample ( ivec &ind ) {
-                est.resample ( ind, resmethod );
+        void resample ( ) {
+                ivec ind = zeros_i ( n );
+                bdm::resample(w,ind,resmethod);
+                // copy the internals according to ind
+                for (int i = 0; i < n; i++ ) {
+                        if ( ind ( i ) != i ) {
+                                particles( i ) = particles( ind ( i ) )->_copy();
+                        }
+                        w ( i ) = 1.0 / n;
+                }
         }
         //! access function
-        Array<vec>& __samples() {
-                return _samples;
-        }
-
-        virtual double logpred ( const vec &yt ) const NOT_IMPLEMENTED(0);
-
-        virtual epdf* epredictor() const NOT_IMPLEMENTED(NULL);
-        
-        virtual pdf* predictor() const NOT_IMPLEMENTED(NULL);
+        Array<BM*>& _particles() {
+                return particles;
+        }
+
 };
 UIREGISTER ( PF );
@@ -222 +397 @@
 */
 
-class MPF : public BM  {
-        //! \var log_level_enums means 
-        //! TODO DOPLNIT
-        LOG_LEVEL(MPF,means);
-
-protected:
-        //! particle filter on non-linear variable
-        shared_ptr<PF> pf;
-        //! Array of Bayesian models
-        Array<BM*> BMs;
-
-        //! internal class for pdf providing composition of eEmp with external components
-
-        class mpfepdf : public epdf  {
-                //! pointer to particle filter
-                shared_ptr<PF> &pf;
-                //! pointer to Array of BMs
-                Array<BM*> &BMs;
-        public:
-                //! constructor
-                mpfepdf ( shared_ptr<PF> &pf0, Array<BM*> &BMs0 ) : epdf(), pf ( pf0 ), BMs ( BMs0 ) { };
-                //! a variant of set parameters - this time, parameters are read from BMs and pf
-                void read_parameters() {
-                        rv = concat ( pf->posterior()._rv(), BMs ( 0 )->posterior()._rv() );
-                        dim = pf->posterior().dimension() + BMs ( 0 )->posterior().dimension();
-                        bdm_assert_debug ( dim == rv._dsize(), "Wrong name " );
-                }
-                vec mean() const;
-
-                vec variance() const;
-
-                void qbounds ( vec &lb, vec &ub, double perc = 0.95 ) const;
-
-                vec sample() const NOT_IMPLEMENTED(0);
-
-                double evallog ( const vec &val ) const NOT_IMPLEMENTED(0);             
-        };
-
-        //! Density joining PF.est with conditional parts
-        mpfepdf jest;
-
-        //! datalink from global yt and cond (Up) to BMs yt and cond (Down)
-        datalink_m2m this2bm;
-        //! datalink from global yt and cond (Up) to PFs yt and cond (Down)
-        datalink_m2m this2pf;
-        //!datalink from PF part to BM
-        datalink_part pf2bm;
-
-public:
-        //! Default constructor.
-        MPF () :  jest ( pf, BMs ) {};
-        //! set all parameters at once
-        void set_pf ( shared_ptr<pdf> par0, int n0, RESAMPLING_METHOD rm = SYSTEMATIC ) {
-                if (!pf) pf=new PF;
-                pf->set_model ( par0, par0 ); // <=== nasty!!!
-                pf->set_parameters ( n0, rm );
-                pf->set_rv(par0->_rv());
-                BMs.set_length ( n0 );
-        }
-        //! set a prototype of BM, copy it to as many times as there is particles in pf
-        void set_BM ( const BM &BMcond0 ) {
-
-                int n = pf->__w().length();
-                BMs.set_length ( n );
-                // copy
-                //BMcond0 .condition ( pf->posterior()._sample ( 0 ) );
-                for ( int i = 0; i < n; i++ ) {
-                        BMs ( i ) = (BM*) BMcond0._copy();
-                }
-        };
-
-        void bayes ( const vec &yt, const vec &cond );
-
-        const epdf& posterior() const {
-                return jest;
-        }
-
-        //!Access function
-        const BM* _BM ( int i ) {
-                return BMs ( i );
-        }
-        PF& _pf() {return *pf;}
-
-
-        virtual double logpred ( const vec &yt ) const NOT_IMPLEMENTED(0); 
-                
-        virtual epdf* epredictor() const NOT_IMPLEMENTED(NULL); 
-        
-        virtual pdf* predictor() const NOT_IMPLEMENTED(NULL); 
-
-
-        /*! configuration structure for basic PF
-        \code
-        BM              = BM_class;           // Bayesian filtr for analytical part of the model
-        parameter_pdf   = pdf_class;         // transitional pdf for non-parametric part of the model
-        prior           = epdf_class;         // prior probability density
-        --- optional ---
-        n               = 10;                 // number of particles
-        resmethod       = 'systematic', or 'multinomial', or 'stratified'
-                                                                                  // resampling method
-        \endcode
-        */
-        void from_setting ( const Setting &set ) {
-                BM::from_setting( set );
-
-                shared_ptr<pdf> par = UI::build<pdf> ( set, "parameter_pdf", UI::compulsory );
-
-                pf = new PF;
-                // prior must be set before BM
-                pf->prior_from_set ( set );
-                pf->resmethod_from_set ( set );
-                pf->set_model ( par, par ); // too hackish!
-
-                shared_ptr<BM> BM0 = UI::build<BM> ( set, "BM", UI::compulsory );
-                set_BM ( *BM0 );
-
-                //set drv
-                //??set_yrv(concat(BM0->_yrv(),u) );
-                set_yrv ( BM0->_yrv() );
-                rvc = BM0->_rvc().subt ( par->_rv() );
-                //find potential input - what remains in rvc when we subtract rv
-                RV u = par->_rvc().subt ( par->_rv().copy_t ( -1 ) );
-                rvc.add ( u );
-                dimc = rvc._dsize();
-                validate();
-        }
-
-        void validate() {
-                BM::validate();
-                try {
-                        pf->validate();
-                } catch ( std::exception ) {
-                        throw UIException ( "Error in PF part of MPF:" );
-                }
-                jest.read_parameters();
-                this2bm.set_connection ( BMs ( 0 )->_yrv(), BMs ( 0 )->_rvc(), yrv, rvc );
-                this2pf.set_connection ( pf->_yrv(), pf->_rvc(), yrv, rvc );
-                pf2bm.set_connection ( BMs ( 0 )->_rvc(), pf->posterior()._rv() );
-        }
-};
-UIREGISTER ( MPF );
+// class MPF : public BM  {
+//      //! Introduces new option 
+//      //! \li means - meaning TODO
+//      LOG_LEVEL(MPF,means);
+// protected:
+//      //! particle filter on non-linear variable
+//      shared_ptr<PF> pf;
+//      //! Array of Bayesian models
+//      Array<BM*> BMs;
+// 
+//      //! internal class for pdf providing composition of eEmp with external components
+// 
+//      class mpfepdf : public epdf  {
+//              //! pointer to particle filter
+//              shared_ptr<PF> &pf;
+//              //! pointer to Array of BMs
+//              Array<BM*> &BMs;
+//      public:
+//              //! constructor
+//              mpfepdf ( shared_ptr<PF> &pf0, Array<BM*> &BMs0 ) : epdf(), pf ( pf0 ), BMs ( BMs0 ) { };
+//              //! a variant of set parameters - this time, parameters are read from BMs and pf
+//              void read_parameters() {
+//                      rv = concat ( pf->posterior()._rv(), BMs ( 0 )->posterior()._rv() );
+//                      dim = pf->posterior().dimension() + BMs ( 0 )->posterior().dimension();
+//                      bdm_assert_debug ( dim == rv._dsize(), "Wrong name " );
+//              }
+//              vec mean() const;
+// 
+//              vec variance() const;
+// 
+//              void qbounds ( vec &lb, vec &ub, double perc = 0.95 ) const;
+// 
+//              vec sample() const NOT_IMPLEMENTED(0);
+// 
+//              double evallog ( const vec &val ) const NOT_IMPLEMENTED(0);             
+//      };
+// 
+//      //! Density joining PF.est with conditional parts
+//      mpfepdf jest;
+// 
+//      //! datalink from global yt and cond (Up) to BMs yt and cond (Down)
+//      datalink_m2m this2bm;
+//      //! datalink from global yt and cond (Up) to PFs yt and cond (Down)
+//      datalink_m2m this2pf;
+//      //!datalink from PF part to BM
+//      datalink_part pf2bm;
+// 
+// public:
+//      //! Default constructor.
+//      MPF () :  jest ( pf, BMs ) {};
+//      //! set all parameters at once
+//      void set_pf ( shared_ptr<pdf> par0, int n0, RESAMPLING_METHOD rm = SYSTEMATIC ) {
+//              if (!pf) pf=new PF;
+//              pf->set_model ( par0, par0 ); // <=== nasty!!!
+//              pf->set_parameters ( n0, rm );
+//              pf->set_rv(par0->_rv());
+//              BMs.set_length ( n0 );
+//      }
+//      //! set a prototype of BM, copy it to as many times as there is particles in pf
+//      void set_BM ( const BM &BMcond0 ) {
+// 
+//              int n = pf->__w().length();
+//              BMs.set_length ( n );
+//              // copy
+//              //BMcond0 .condition ( pf->posterior()._sample ( 0 ) );
+//              for ( int i = 0; i < n; i++ ) {
+//                      BMs ( i ) = (BM*) BMcond0._copy();
+//              }
+//      };
+// 
+//      void bayes ( const vec &yt, const vec &cond );
+// 
+//      const epdf& posterior() const {
+//              return jest;
+//      }
+// 
+//      //!Access function
+//      const BM* _BM ( int i ) {
+//              return BMs ( i );
+//      }
+//      PF& _pf() {return *pf;}
+// 
+// 
+//      virtual double logpred ( const vec &yt ) const NOT_IMPLEMENTED(0); 
+//              
+//      virtual epdf* epredictor() const NOT_IMPLEMENTED(NULL); 
+//      
+//      virtual pdf* predictor() const NOT_IMPLEMENTED(NULL); 
+// 
+// 
+//      /*! configuration structure for basic PF
+//      \code
+//      BM              = BM_class;           // Bayesian filtr for analytical part of the model
+//      parameter_pdf   = pdf_class;         // transitional pdf for non-parametric part of the model
+//      prior           = epdf_class;         // prior probability density
+//      --- optional ---
+//      n               = 10;                 // number of particles
+//      resmethod       = 'systematic', or 'multinomial', or 'stratified'
+//                                                                                // resampling method
+//      \endcode
+//      */
+//      void from_setting ( const Setting &set ) {
+//              BM::from_setting( set );
+// 
+//              shared_ptr<pdf> par = UI::build<pdf> ( set, "parameter_pdf", UI::compulsory );
+// 
+//              pf = new PF;
+//              // prior must be set before BM
+//              pf->prior_from_set ( set );
+//              pf->resmethod_from_set ( set );
+//              pf->set_model ( par, par ); // too hackish!
+// 
+//              shared_ptr<BM> BM0 = UI::build<BM> ( set, "BM", UI::compulsory );
+//              set_BM ( *BM0 );
+// 
+//              //set drv
+//              //??set_yrv(concat(BM0->_yrv(),u) );
+//              set_yrv ( BM0->_yrv() );
+//              rvc = BM0->_rvc().subt ( par->_rv() );
+//              //find potential input - what remains in rvc when we subtract rv
+//              RV u = par->_rvc().subt ( par->_rv().copy_t ( -1 ) );
+//              rvc.add ( u );
+//              dimc = rvc._dsize();
+//              validate();
+//      }
+// 
+//      void validate() {
+//              BM::validate();
+//              try {
+//                      pf->validate();
+//              } catch ( std::exception ) {
+//                      throw UIException ( "Error in PF part of MPF:" );
+//              }
+//              jest.read_parameters();
+//              this2bm.set_connection ( BMs ( 0 )->_yrv(), BMs ( 0 )->_rvc(), yrv, rvc );
+//              this2pf.set_connection ( pf->_yrv(), pf->_rvc(), yrv, rvc );
+//              pf2bm.set_connection ( BMs ( 0 )->_rvc(), pf->posterior()._rv() );
+//      }
+// };
+// UIREGISTER ( MPF );
 
 /*! ARXg for estimation of state-space variances
 */
-class MPF_ARXg :public BM{
-        protected:
-        shared_ptr<PF> pf;
-        //! pointer to Array of BMs
-        Array<ARX*> BMso;
-        //! pointer to Array of BMs
-        Array<ARX*> BMsp;
-        //!parameter evolution
-        shared_ptr<fnc> g;
-        //!observation function
-        shared_ptr<fnc> h;
-        
-        public:
-                void bayes(const vec &yt, const vec &cond );
-                void from_setting(const Setting &set) ;
-                void validate() {
-                        bdm_assert(g->dimensionc()==g->dimension(),"not supported yet");
-                        bdm_assert(h->dimensionc()==g->dimension(),"not supported yet");                        
-                }
-
-                double logpred(const vec &cond) const NOT_IMPLEMENTED(0.0);
-                epdf* epredictor() const NOT_IMPLEMENTED(NULL);
-                pdf* predictor() const NOT_IMPLEMENTED(NULL);
-                const epdf& posterior() const {return pf->posterior();};
-                
-                void log_register( logger &L, const string &prefix ){
-                        BM::log_register(L,prefix);
-                        logrec->ids.set_size ( 3 );
-                        logrec->ids(1)= L.add_vector(RV("Q",dimension()*dimension()), prefix+L.prefix_sep()+"Q");
-                        logrec->ids(2)= L.add_vector(RV("R",dimensiony()*dimensiony()), prefix+L.prefix_sep()+"R");
-                        
-                };
-                void log_write() const {
-                        BM::log_write();
-                        mat mQ=zeros(dimension(),dimension());
-                        mat pom=zeros(dimension(),dimension());
-                        mat mR=zeros(dimensiony(),dimensiony());
-                        mat pom2=zeros(dimensiony(),dimensiony());
-                        mat dum;
-                        const vec w=pf->posterior()._w();
-                        for (int i=0; i<w.length(); i++){
-                                BMsp(i)->posterior().mean_mat(dum,pom);
-                                mQ += w(i) * pom;
-                                BMso(i)->posterior().mean_mat(dum,pom2);
-                                mR += w(i) * pom2;
-                                
-                        }
-                        logrec->L.log_vector ( logrec->ids ( 1 ), cvectorize(mQ) );
-                        logrec->L.log_vector ( logrec->ids ( 2 ), cvectorize(mR) );
-                        
-                }
-};
-UIREGISTER(MPF_ARXg);
+// class MPF_ARXg :public BM{
+//      protected:
+//      shared_ptr<PF> pf;
+//      //! pointer to Array of BMs
+//      Array<ARX*> BMso;
+//      //! pointer to Array of BMs
+//      Array<ARX*> BMsp;
+//      //!parameter evolution
+//      shared_ptr<fnc> g;
+//      //!observation function
+//      shared_ptr<fnc> h;
+//      
+//      public:
+//              void bayes(const vec &yt, const vec &cond );
+//              void from_setting(const Setting &set) ;
+//              void validate() {
+//                      bdm_assert(g->dimensionc()==g->dimension(),"not supported yet");
+//                      bdm_assert(h->dimensionc()==g->dimension(),"not supported yet");                        
+//              }
+// 
+//              double logpred(const vec &cond) const NOT_IMPLEMENTED(0.0);
+//              epdf* epredictor() const NOT_IMPLEMENTED(NULL);
+//              pdf* predictor() const NOT_IMPLEMENTED(NULL);
+//              const epdf& posterior() const {return pf->posterior();};
+//              
+//              void log_register( logger &L, const string &prefix ){
+//                      BM::log_register(L,prefix);
+//                      logrec->ids.set_size ( 3 );
+//                      logrec->ids(1)= L.add_vector(RV("Q",dimension()*dimension()), prefix+L.prefix_sep()+"Q");
+//                      logrec->ids(2)= L.add_vector(RV("R",dimensiony()*dimensiony()), prefix+L.prefix_sep()+"R");
+//                      
+//              };
+//              void log_write() const {
+//                      BM::log_write();
+//                      mat mQ=zeros(dimension(),dimension());
+//                      mat pom=zeros(dimension(),dimension());
+//                      mat mR=zeros(dimensiony(),dimensiony());
+//                      mat pom2=zeros(dimensiony(),dimensiony());
+//                      mat dum;
+//                      const vec w=pf->posterior()._w();
+//                      for (int i=0; i<w.length(); i++){
+//                              BMsp(i)->posterior().mean_mat(dum,pom);
+//                              mQ += w(i) * pom;
+//                              BMso(i)->posterior().mean_mat(dum,pom2);
+//                              mR += w(i) * pom2;
+//                              
+//                      }
+//                      logrec->L.log_vector ( logrec->ids ( 1 ), cvectorize(mQ) );
+//                      logrec->L.log_vector ( logrec->ids ( 2 ), cvectorize(mR) );
+//                      
+//              }
+// };
+// UIREGISTER(MPF_ARXg);
 
 

library/bdm/stat/exp_family.cpp

@@ -423 +423 @@
 }
 
-void eEmp::resample ( ivec &ind, RESAMPLING_METHOD method ) {
+void eEmp::resample ( RESAMPLING_METHOD method ) {
+        ivec ind = zeros_i ( n );
+        bdm::resample(w,ind,method);
+        // copy the internals according to ind
+        for (int i = 0; i < n; i++ ) {
+                if ( ind ( i ) != i ) {
+                        samples ( i ) = samples ( ind ( i ) );
+                }
+                w ( i ) = 1.0 / n;
+        }
+}
+
+void resample ( const vec &w, ivec &ind, RESAMPLING_METHOD method ) {
+        int n = w.length();
         ind = zeros_i ( n );
         ivec N_babies = zeros_i ( n );
@@ -430 +443 @@
         int i, j, parent;
         double u0;
-
+        
         switch ( method ) {
-        case MULTINOMIAL:
-                u ( n - 1 ) = pow ( UniRNG.sample(), 1.0 / n );
-
-                for ( i = n - 2; i >= 0; i-- ) {
-                        u ( i ) = u ( i + 1 ) * pow ( UniRNG.sample(), 1.0 / ( i + 1 ) );
-                }
-
-                break;
-
-        case STRATIFIED:
-
-                for ( i = 0; i < n; i++ ) {
-                        u ( i ) = ( i + UniRNG.sample() ) / n;
-                }
-
-                break;
-
-        case SYSTEMATIC:
-                u0 = UniRNG.sample();
-
-                for ( i = 0; i < n; i++ ) {
-                        u ( i ) = ( i + u0 ) / n;
-                }
-
-                break;
-
-        default:
-                bdm_error ( "PF::resample(): Unknown resampling method" );
-        }
-
+                case MULTINOMIAL:
+                        u ( n - 1 ) = pow ( UniRNG.sample(), 1.0 / n );
+                        
+                        for ( i = n - 2; i >= 0; i-- ) {
+                                u ( i ) = u ( i + 1 ) * pow ( UniRNG.sample(), 1.0 / ( i + 1 ) );
+                        }
+                        
+                        break;
+                        
+                case STRATIFIED:
+                        
+                        for ( i = 0; i < n; i++ ) {
+                                u ( i ) = ( i + UniRNG.sample() ) / n;
+                        }
+                        
+                        break;
+                        
+                case SYSTEMATIC:
+                        u0 = UniRNG.sample();
+                        
+                        for ( i = 0; i < n; i++ ) {
+                                u ( i ) = ( i + u0 ) / n;
+                        }
+                        
+                        break;
+                        
+                default:
+                        bdm_error ( "PF::resample(): Unknown resampling method" );
+        }
+        
         // U is now full
         j = 0;
-
+        
         for ( i = 0; i < n; i++ ) {
                 while ( u ( i ) > cumDist ( j ) ) j++;
-
+                
                 N_babies ( j ) ++;
         }
@@ -474 +487 @@
         // * particles with at least one baby should not move *
         // This assures that reassignment can be done inplace;
-
+        
         // find the first parent;
         parent = 0;
         while ( N_babies ( parent ) == 0 ) parent++;
-
+        
         // Build index
         for ( i = 0; i < n; i++ ) {
@@ -488 +501 @@
                         // if yes, find the next one
                         while ( ( N_babies ( parent ) == 0 ) || ( N_babies ( parent ) == 1 && parent > i ) ) parent++;
-
+                        
                         // Replicate parent
                         ind ( i ) = parent;
-
+                        
                         N_babies ( parent ) --; //this index was now replicated;
                 }
-
-        }
-
-// copy the internals according to ind
-        for ( i = 0; i < n; i++ ) {
-                if ( ind ( i ) != i ) {
-                        samples ( i ) = samples ( ind ( i ) );
-                }
-                w ( i ) = 1.0 / n;
+                
         }
 }

library/bdm/stat/exp_family.h

@@ -97 +97 @@
         virtual void flatten ( const BMEF * B ) NOT_IMPLEMENTED_VOID;
 
-        double logpred ( const vec &yt ) const NOT_IMPLEMENTED(0);
-        
-        virtual epdf* epredictor() const NOT_IMPLEMENTED(NULL);
-
-        virtual pdf* predictor() const NOT_IMPLEMENTED(NULL);
 
         void to_setting ( Setting &set ) const
@@ -267 +262 @@
 UIREGISTER2 ( enorm, fsqmat );
 SHAREDPTR2 ( enorm, fsqmat );
+
+typedef enorm<ldmat> egauss;
+UIREGISTER(egauss);
+
 
 //forward declaration
@@ -1596 +1595 @@
 //! Switch between various resampling methods.
 enum RESAMPLING_METHOD { MULTINOMIAL = 0, STRATIFIED = 1, SYSTEMATIC = 3 };
+
+//! Shortcut for multinomial.sample(int n). Various simplifications are allowed see RESAMPLING_METHOD 
+void resample(const vec &w, ivec &ind, RESAMPLING_METHOD=SYSTEMATIC);
+
 /*!
 \brief Weighted empirical density
@@ -1664 +1667 @@
         };
         //! Function performs resampling, i.e. removal of low-weight samples and duplication of high-weight samples such that the new samples represent the same density.
-        //! The vector with indeces of new samples is returned in variable \c index.
-        void resample ( ivec &index, RESAMPLING_METHOD method = SYSTEMATIC );
-
-        //! Resampling without returning index of new particles.
-        void resample ( RESAMPLING_METHOD method = SYSTEMATIC ) {
-                ivec ind;
-                resample ( ind, method );
-        };
+        void resample ( RESAMPLING_METHOD method = SYSTEMATIC );
 
         //! inherited operation : NOT implemented
@@ -1874 +1870 @@
 }
 
+/*! Dirac delta function distribution */
+class dirac: public epdf{
+        protected: 
+                vec point;
+        public:
+                double evallog (const vec &dt) const {return -inf;}
+                vec mean () const {return point;}
+                vec variance () const {return pow(point,2);}
+                void qbounds ( vec &lb, vec &ub, double percentage = 0.95 ) const { lb = point; ub = point;}
+                //! access
+                const vec& _point() {return point;}
+                void set_point(const vec& p){point =p; dim=p.length();}
+                vec sample() const {return point;}
+        };
+
 ////
 ///////