Changeset 887 for library/bdm/stat/exp_family.cpp

Timestamp:

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

Author:

smidl

Message:

new base for particle filtering

Files:

• library/bdm/stat/exp_family.cpp (modified) (4 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;
+                
         }
 }