Changeset 60 for bdm/estim/libPF.h

Timestamp:

03/31/08 13:53:36 (16 years ago)

Author:

smidl

Message:

nove rozlozeni Q

Files:

• bdm/estim/libPF.h (modified) (9 diffs)

bdm/estim/libPF.h

@@ -19 +19 @@
 
 using namespace itpp;
+
+//UGLY HACK
+extern double PF_SSAT;//used for StrSim:06 test... if length>0 the value is written.
 
 /*!
@@ -60 +63 @@
 
 class MPF : public PF {
-        BM_T* Bms[1000];
+        BM_T* Bms[10000];
 
         //! internal class for MPDF providing composition of eEmp with external components
@@ -105 +108 @@
                 //
 
-                if ( n>1000 ) it_error ( "increase 1000 here!" );
+                if ( n>10000 ) it_error ( "increase 10000 here!" );
 
                 for ( int i=0;i<n;i++ ) {
@@ -126 +129 @@
                 for ( int i=0;i<n;i++ ) {Bms[i]->condition ( _samples ( i ) );}
         }
+
+//SimStr:
+        double SSAT;    
 };
 
@@ -132 +138 @@
         int i;
         vec lls ( n );
+        vec llsP ( n );
         ivec ind;
         double mlls=-std::numeric_limits<double>::infinity();
 
+        // StrSim:06
+        double sumLWL=0.0;
+        double sumL2WL=0.0;
+        double WL = 0.0;
+
         for ( i=0;i<n;i++ ) {
                 //generate new samples from paramater evolution model;
-                _samples ( i ) = par.samplecond ( _samples ( i ), lls ( i ) );
+                _samples ( i ) = par.samplecond ( _samples ( i ), llsP ( i ) );
                 Bms[i]->condition ( _samples ( i ) );
                 Bms[i]->bayes ( dt );
@@ -144 +156 @@
                 if ( lls ( i ) >mlls ) mlls=lls ( i ); //find maximum likelihood (for numerical stability)
         }
+        
+        if ( true ) {
+                for ( i=0;i<n;i++ ) {
+                        WL = _w ( i ) *exp ( llsP ( i ) ); //using old weights!
+                        sumLWL += exp ( lls ( i ) ) *WL;
+                        sumL2WL += exp ( 2*lls ( i ) ) *WL;
+                }
+                SSAT  = sumL2WL/(sumLWL*sumLWL);
+        }
 
         // compute weights
@@ -150 +171 @@
         }
 
-        _w /=sum(_w); //?
-
-        double eff = 1.0/( _w*_w );
+        _w /=sum ( _w ); //?
+
+
+
+        double eff = 1.0/ ( _w*_w );
         if ( eff < ( 0.1*n ) ) {
                 ind = est.resample();
@@ -161 +184 @@
                                 //fixme this would require new assignment operator
                                 // *Bms[i] = *Bms[ind ( i ) ];
-                                
+
                                 // poor-man's solution: replicate constructor here
                                 // copied from MPF::MPF
@@ -170 +193 @@
                         }
                 };
-                cout << '.'; 
+                cout << '.';
         }
 }