title: Implementation of sequential monte carlo method (particle filters)
body: <p>I'm interested in the simple algorithm for particles filter given here: <a href="http://www.aiqus.com/upfiles/PFAlgo.png">http://www.aiqus.com/upfiles/PFAlgo.png</a> It seems very simple but I have no idea on how to do it practically.
  Any idea on how to implement it (just to better understand how it works) ?</p>
  
  <p><strong>Edit:</strong>
  This is a great simple example that explain how it works: <a href="http://www.aiqus.com/questions/39942/very-simple-particle-filters-algorithm-sequential-monte-carlo-method-implementation?page=1#39950">http://www.aiqus.com/questions/39942/very-simple-particle-filters-algorithm-sequential-monte-carlo-method-implementation?page=1#39950</a></p>
  
  <p>I've tried to implement it in C++: <a href="http://pastebin.com/M1q1HcN4">http://pastebin.com/M1q1HcN4</a> but I'm note sure if I do it the right way. Can you please check if I understood it well, or there are some misunderstanding according to my code ?</p>
  
  <pre><code>#include &lt;iostream&gt;
  #include &lt;vector&gt;
  #include &lt;boost/random/mersenne_twister.hpp&gt;
  #include &lt;boost/random/uniform_01.hpp&gt;
  #include &lt;boost/random/uniform_int_distribution.hpp&gt;
  
  using namespace std;
  using namespace boost;
  
  double uniform_generator(void);
  
  #define N 4 // number of particles
  
  #define evolutionProba_A_A 1.0/3.0 // P(X_t = A | X_t-1 = A)
  #define evolutionProba_A_B 1.0/3.0 // P(X_t = A | X_t-1 = B)
  #define evolutionProba_B_B 2.0/3.0 // P(X_t = B | X_t-1 = B)
  #define evolutionProba_B_A 2.0/3.0 // P(X_t = B | X_t-1 = A)
  
  #define observationProba_A_A 4.0/5.0 // P(Y_t = A | X_t = A)
  #define observationProba_A_B 1.0/5.0 // P(Y_t = A | X_t = B)
  #define observationProba_B_B 4.0/5.0 // P(Y_t = B | X_t = B)
  #define observationProba_B_A 1.0/5.0 // P(Y_t = A | X_t = A)
  
  /// ===========================================================================
  
  typedef struct distrib { float PA; float PB; } Distribution;
  
  typedef struct particle
  {
      Distribution distribution; // e.g. &lt;0.5, 0.5&gt;
      char state; // e.g. 'A' or 'B'
      float weight; // e.g. 0.8
  }
  Particle;
  
  /// ===========================================================================
  
  int main()
  {
      vector&lt;char&gt; Y; // data observations
      Y.push_back('A'); Y.push_back('B'); Y.push_back('A'); Y.push_back('A'); Y.push_back('A'); Y.push_back('B');
      Y.push_back('A'); Y.push_back('A'); Y.push_back('B'); Y.push_back('A'); Y.push_back('B'); Y.push_back('A');
      Y.push_back('A'); Y.push_back('B'); Y.push_back('B'); Y.push_back('A'); Y.push_back('A'); Y.push_back('B');
  
      vector&lt; vector&lt;Particle&gt; &gt; Xall; // vector of all particles from time 0 to t
  
      /// Step (1) Initialisation
      vector&lt;Particle&gt; X; // a vector of N particles
      for(int i = 0; i &lt; N; ++i)
      {
          Particle x;
  
          // sample particle Xi from initial distribution
          x.distribution.PA = 0.5; x.distribution.PB = 0.5;
          float r = uniform_generator();
          if( r &lt;= x.distribution.PA ) x.state = 'A'; // r &lt;= 0.5
          if( x.distribution.PA &lt; r &amp;&amp; r &lt;= x.distribution.PA + x.distribution.PB ) x.state = 'B'; // 0.5 &lt; r &lt;= 1
  
          X.push_back(x);
      }
  
      Xall.push_back(X);
      X.clear();
  
      /// Observing data
      for(int t = 1; t &lt;= 18; ++t)
      {
          char y = Y[t-1]; // current observation
  
          /// Step (2) Importance sampling
          float sumWeights = 0;
          vector&lt;Particle&gt; X; // a vector of N particles
          for(int i = 0; i &lt; N; ++i)
          {
              Particle x;
  
              // P(X^i_t = A) = P(X^i_t = A | X^i_t-1 = A) * P(X^i_t-1 = A) + P(X^i_t = A | X^i_t-1 = B) * P(X^i_t-1 = B)
              x.distribution.PA = evolutionProba_A_A * Xall[t-1][i].distribution.PA + evolutionProba_A_B * Xall[t-1][i].distribution.PB;
  
              // P(X^i_t = B) = P(X^i_t = B | X^i_t-1 = A) * P(X^i_t-1 = A) + P(X^i_t = B | X^i_t-1 = B) * P(X^i_t-1 = B)
              x.distribution.PB = evolutionProba_B_A * Xall[t-1][i].distribution.PA + evolutionProba_B_B * Xall[t-1][i].distribution.PB;
  
              // sample the a particle from this distribution
              float r = uniform_generator();
              if( r &lt;= x.distribution.PA ) x.state = 'A';
              if( x.distribution.PA &lt; r &amp;&amp; r &lt;= x.distribution.PA + x.distribution.PB ) x.state = 'B';
  
              // compute weight of this particle according to the observation y
              if( y == 'A' )
              {
                  if( x.state == 'A' ) x.weight = observationProba_A_A; // P(y = A | X^i_t = A)
                  else if( x.state == 'B' ) x.weight = observationProba_A_B; // P(y = A | X^i_t = B)
              }
              else if( y == 'B' )
              {
                  if( x.state == 'A' ) x.weight = observationProba_B_A; // P(y = B | X^i_t = A)
                  else if( x.state == 'B' ) x.weight = observationProba_B_B; // P(y = B | X^i_t = B)
              }
  
              sumWeights += x.weight;
  
              X.push_back(x);
          }
  
          // normalise weights
          for(int i = 0; i &lt; N; ++i)
              X[i].weight /= sumWeights;
  
          /// Step (3) resampling N particles according to weights
          float PA = 0, PB = 0;
          for(int i = 0; i &lt; N; ++i)
          {
              if( X[i].state == 'A' ) PA += X[i].weight;
              else if( X[i].state == 'B' ) PB += X[i].weight;
          }
  
          vector&lt;Particle&gt; reX; // new vector of particles
          for(int i = 0; i &lt; N; ++i)
          {
              Particle x;
  
              x.distribution.PA = PA;
              x.distribution.PB = PB;
  
              float r = uniform_generator();
              if( r &lt;= x.distribution.PA ) x.state = 'A';
              if( x.distribution.PA &lt; r &amp;&amp; r &lt;= x.distribution.PA + x.distribution.PB ) x.state = 'B';
  
              reX.push_back(x);
          }
  
          Xall.push_back(reX);
      }
  
      return 0;
  }
  
  /// ===========================================================================
  
  double uniform_generator(void)
  {
      mt19937 gen(55);
      static uniform_01&lt; mt19937, double &gt; uniform_gen(gen);
      return uniform_gen();
  }
  </code></pre>