Changeset 887 for library/bdm/stat

Timestamp:

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

Author:

smidl

Message:

new base for particle filtering

Location:

library/bdm/stat

Files:

• exp_family.cpp (modified) (4 diffs)
• exp_family.h (modified) (5 diffs)

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;}
+        };
+
 ////
 ///////