13
Tôi quan tâm đến thuật toán đơn giản cho bộ lọc hạt được đưa ra ở đây: http://www.aiqus.com/upfiles/PFAlgo.png Có vẻ như rất đơn giản nhưng tôi không có ý tưởng về cách thực hiện nó. Bất kỳ ý tưởng nào về cách thực hiện nó (chỉ để hiểu rõ hơn cách hoạt động)?Thực hiện phương pháp monte carlo tuần tự (bộ lọc hạt)
Edit: Đây là một ví dụ đơn giản tuyệt vời mà giải thích cách hoạt động: http://www.aiqus.com/questions/39942/very-simple-particle-filters-algorithm-sequential-monte-carlo-method-implementation?page=1#39950
Tôi đã cố gắng để thực hiện nó trong C++: http://pastebin.com/M1q1HcN4 nhưng tôi lưu ý chắc chắn nếu tôi làm điều đó một cách đúng đắn . Bạn có thể vui lòng kiểm tra xem tôi có hiểu rõ không, hoặc có một số hiểu lầm theo mã của tôi không?
#include <iostream>
#include <vector>
#include <boost/random/mersenne_twister.hpp>
#include <boost/random/uniform_01.hpp>
#include <boost/random/uniform_int_distribution.hpp>
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. <0.5, 0.5>
char state; // e.g. 'A' or 'B'
float weight; // e.g. 0.8
}
Particle;
/// ===========================================================================
int main()
{
vector<char> 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< vector<Particle> > Xall; // vector of all particles from time 0 to t
/// Step (1) Initialisation
vector<Particle> X; // a vector of N particles
for(int i = 0; i < 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 <= x.distribution.PA) x.state = 'A'; // r <= 0.5
if(x.distribution.PA < r && r <= x.distribution.PA + x.distribution.PB) x.state = 'B'; // 0.5 < r <= 1
X.push_back(x);
}
Xall.push_back(X);
X.clear();
/// Observing data
for(int t = 1; t <= 18; ++t)
{
char y = Y[t-1]; // current observation
/// Step (2) Importance sampling
float sumWeights = 0;
vector<Particle> X; // a vector of N particles
for(int i = 0; i < 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 <= x.distribution.PA) x.state = 'A';
if(x.distribution.PA < r && r <= 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 < N; ++i)
X[i].weight /= sumWeights;
/// Step (3) resampling N particles according to weights
float PA = 0, PB = 0;
for(int i = 0; i < N; ++i)
{
if(X[i].state == 'A') PA += X[i].weight;
else if(X[i].state == 'B') PB += X[i].weight;
}
vector<Particle> reX; // new vector of particles
for(int i = 0; i < N; ++i)
{
Particle x;
x.distribution.PA = PA;
x.distribution.PB = PB;
float r = uniform_generator();
if(r <= x.distribution.PA) x.state = 'A';
if(x.distribution.PA < r && r <= 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< mt19937, double > uniform_gen(gen);
return uniform_gen();
}
Khi nào bạn có thể lọc bộ lọc này trong thế giới thực? Bạn có thể chạy thử nghiệm đối với một vấn đề với giải pháp phân tích không? Nếu bạn đã thực hiện nó một cách chính xác, bạn sẽ nhận được cùng một số. Rất khó để có được kết quả đúng với việc triển khai sai! –
@AlessandroTeruzzi Đây chỉ là một ví dụ đơn giản được đưa ra [ở đây] (http://www.aiqus.com/questions/39942/very-simple-particle-filters-algorithm-sequential-monte-carlo-method-implementation/39950). Tôi đã thực hiện nó chỉ để hiểu rõ hơn về khái niệm [bộ lọc hạt được đưa ra bởi thuật toán này] (http://www.aiqus.com/upfiles/PFAlgo.png), nhưng tôi không chắc chắn nếu tôi đã triển khai nó một cách chính xác, vì Tôi không hiểu thuật toán rất rõ. Tôi không biết làm thế nào để kiểm tra nếu nó hoạt động, vì thuật toán và đầu ra của nó vẫn không rõ ràng với tôi (ngay cả khi thuật toán dường như rất đơn giản). – shn
Đề nghị đầu tiên của tôi cho một thuật toán chung: không cố gắng thực hiện một cái gì đó bạn không hoàn toàn hiểu. Đầu tiên cam kết, sau đó thực hiện. Nếu không, bạn sẽ không thể nói những gì đang xảy ra. –