Check out the Pastebin Gadgets Shop. We have thousands of fun, geeky & affordable gadgets on sale :-)Want more features on Pastebin? Sign Up, it's FREE!

# Untitled

By: a guest on May 10th, 2012  |  syntax: C++  |  size: 5.02 KB  |  views: 2,307  |  expires: Never
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. }
clone this paste RAW Paste Data
Top