Pastebin launched a little side project called VERYVIRAL.com, check it out ;-) Want more features on Pastebin? Sign Up, it's FREE!
Guest

Untitled

By: a guest on May 10th, 2012  |  syntax: C++  |  size: 5.02 KB  |  views: 1,491  |  expires: Never
download  |  raw  |  embed  |  report abuse  |  print
Text below is selected. Please press Ctrl+C to copy to your clipboard. (⌘+C on Mac)
  1. #include <iostream>
  2. #include <vector>
  3. #include <boost/random/mersenne_twister.hpp>
  4. #include <boost/random/uniform_01.hpp>
  5. #include <boost/random/uniform_int_distribution.hpp>
  6.  
  7. using namespace std;
  8. using namespace boost;
  9.  
  10. double uniform_generator(void);
  11.  
  12. #define N 4 // number of particles
  13.  
  14. #define evolutionProba_A_A 1.0/3.0 // P(X_t = A | X_t-1 = A)
  15. #define evolutionProba_A_B 1.0/3.0 // P(X_t = A | X_t-1 = B)
  16. #define evolutionProba_B_B 2.0/3.0 // P(X_t = B | X_t-1 = B)
  17. #define evolutionProba_B_A 2.0/3.0 // P(X_t = B | X_t-1 = A)
  18.  
  19. #define observationProba_A_A 4.0/5.0 // P(Y_t = A | X_t = A)
  20. #define observationProba_A_B 1.0/5.0 // P(Y_t = A | X_t = B)
  21. #define observationProba_B_B 4.0/5.0 // P(Y_t = B | X_t = B)
  22. #define observationProba_B_A 1.0/5.0 // P(Y_t = A | X_t = A)
  23.  
  24. /// ===========================================================================
  25.  
  26. typedef struct distrib { float PA; float PB; } Distribution;
  27.  
  28. typedef struct particle
  29. {
  30.     Distribution distribution; // e.g. <0.5, 0.5>
  31.     char state; // e.g. 'A' or 'B'
  32.     float weight; // e.g. 0.8
  33. }
  34. Particle;
  35.  
  36. /// ===========================================================================
  37.  
  38. int main()
  39. {
  40.     vector<char> Y; // data observations
  41.     Y.push_back('A'); Y.push_back('B'); Y.push_back('A'); Y.push_back('A'); Y.push_back('A'); Y.push_back('B');
  42.     Y.push_back('A'); Y.push_back('A'); Y.push_back('B'); Y.push_back('A'); Y.push_back('B'); Y.push_back('A');
  43.     Y.push_back('A'); Y.push_back('B'); Y.push_back('B'); Y.push_back('A'); Y.push_back('A'); Y.push_back('B');
  44.  
  45.     vector< vector<Particle> > Xall; // vector of all particles from time 0 to t
  46.  
  47.     /// Step (1) Initialisation
  48.     vector<Particle> X; // a vector of N particles
  49.     for(int i = 0; i < N; ++i)
  50.     {
  51.         Particle x;
  52.  
  53.         // sample particle Xi from initial distribution
  54.         x.distribution.PA = 0.5; x.distribution.PB = 0.5;
  55.         float r = uniform_generator();
  56.         if( r <= x.distribution.PA ) x.state = 'A'; // r <= 0.5
  57.         if( x.distribution.PA < r && r <= x.distribution.PA + x.distribution.PB ) x.state = 'B'; // 0.5 < r <= 1
  58.  
  59.         X.push_back(x);
  60.     }
  61.  
  62.     Xall.push_back(X);
  63.     X.clear();
  64.  
  65.     /// Observing data
  66.     for(int t = 1; t <= 18; ++t)
  67.     {
  68.         char y = Y[t-1]; // current observation
  69.  
  70.         /// Step (2) Importance sampling
  71.         float sumWeights = 0;
  72.         vector<Particle> X; // a vector of N particles
  73.         for(int i = 0; i < N; ++i)
  74.         {
  75.             Particle x;
  76.  
  77.             // 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)
  78.             x.distribution.PA = evolutionProba_A_A * Xall[t-1][i].distribution.PA + evolutionProba_A_B * Xall[t-1][i].distribution.PB;
  79.  
  80.             // 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)
  81.             x.distribution.PB = evolutionProba_B_A * Xall[t-1][i].distribution.PA + evolutionProba_B_B * Xall[t-1][i].distribution.PB;
  82.  
  83.             // sample the a particle from this distribution
  84.             float r = uniform_generator();
  85.             if( r <= x.distribution.PA ) x.state = 'A';
  86.             if( x.distribution.PA < r && r <= x.distribution.PA + x.distribution.PB ) x.state = 'B';
  87.  
  88.             // compute weight of this particle according to the observation y
  89.             if( y == 'A' )
  90.             {
  91.                 if( x.state == 'A' ) x.weight = observationProba_A_A; // P(y = A | X^i_t = A)
  92.                 else if( x.state == 'B' ) x.weight = observationProba_A_B; // P(y = A | X^i_t = B)
  93.             }
  94.             else if( y == 'B' )
  95.             {
  96.                 if( x.state == 'A' ) x.weight = observationProba_B_A; // P(y = B | X^i_t = A)
  97.                 else if( x.state == 'B' ) x.weight = observationProba_B_B; // P(y = B | X^i_t = B)
  98.             }
  99.  
  100.             sumWeights += x.weight;
  101.  
  102.             X.push_back(x);
  103.         }
  104.  
  105.         // normalise weights
  106.         for(int i = 0; i < N; ++i)
  107.             X[i].weight /= sumWeights;
  108.  
  109.         /// Step (3) resampling N particles according to weights
  110.         float PA = 0, PB = 0;
  111.         for(int i = 0; i < N; ++i)
  112.         {
  113.             if( X[i].state == 'A' ) PA += X[i].weight;
  114.             else if( X[i].state == 'B' ) PB += X[i].weight;
  115.         }
  116.  
  117.         vector<Particle> reX; // new vector of particles
  118.         for(int i = 0; i < N; ++i)
  119.         {
  120.             Particle x;
  121.  
  122.             x.distribution.PA = PA;
  123.             x.distribution.PB = PB;
  124.  
  125.             float r = uniform_generator();
  126.             if( r <= x.distribution.PA ) x.state = 'A';
  127.             if( x.distribution.PA < r && r <= x.distribution.PA + x.distribution.PB ) x.state = 'B';
  128.  
  129.             reX.push_back(x);
  130.         }
  131.  
  132.         Xall.push_back(reX);
  133.     }
  134.  
  135.     return 0;
  136. }
  137.  
  138. /// ===========================================================================
  139.  
  140. double uniform_generator(void)
  141. {
  142.     mt19937 gen(55);
  143.     static uniform_01< mt19937, double > uniform_gen(gen);
  144.     return uniform_gen();
  145. }