14

ここにあるパーティクル フィルタの単純なアルゴリズムに興味があります。それを実装する方法についてのアイデアはありますか（それがどのように機能するかをよりよく理解するためだけに）？

編集: これは、それがどのように機能するかを説明する非常に単純な例です: 1#39950

C++ で実装しようとしました: http://pastebin.com/M1q1HcN4しかし、正しい方法で実装しているかどうかはわかりません。私がよく理解しているかどうか、または私のコードによる誤解があるかどうかを確認していただけますか?

#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();
}
4

2 に答える 2

21

このオンライン コースは非常に理解しやすく、粒子フィルターについて非常によく説明されています。

それは「ロボット車のプログラミング」と呼ばれ、位置特定の 3 つの方法について説明しています。モンテカルロ位置特定、カルマン フィルター、および粒子フィルターです。

コースは完全に無料です (現在は終了しているため、積極的に参加することはできませんが、講義を視聴することはできます)、スタンフォード大学の教授が教えます。「クラス」は驚くほど簡単に理解できました。小さな演習もありました。一部は論理的なものでしたが、かなりの部分がプログラミングでした。また、あなたが遊ぶことができる宿題。

実際には、すべてのフィルターに対して Python で独自のコードを作成し、テストも行います。コースの終わりまでに、3 つのフィルターすべてを Python に実装する必要があります。

チェックアウトすることを強くお勧めします。

于 2012-05-15T12:43:42.763 に答える
3

ソースへのいくつかの良いリンクは次のとおりです。

Matlab：http ://wiki.esl.fim.uni-passau.de/index.php/Filtering_III:_Kalman_and_Particle_Filter

C：http ://www.robots.ox.ac.uk/~misard/condensation.html （凝縮アルゴリズム）

C ++： http: //robotics.usc.edu/~boyoon/particle.html

Java：http ://www.oursland.net/projects/particlefilter/

于 2012-05-04T20:48:07.837 に答える