logo标签列表

顺序蒙特卡罗方法的实现（粒子滤波器）

c++probabilitymontecarloinferenceparticle-filter
14
我对这里提供的粒子滤波器简单算法很感兴趣：http://www.aiqus.com/upfiles/PFAlgo.png 它看起来非常简单，但我不知道如何在实践中使用它。
你有没有关于如何实现它的任何想法（只是为了更好地理解它是如何工作的）？ 编辑： 这是一个很好的简单示例，解释了它是如何工作的：http://www.aiqus.com/questions/39942/very-simple-particle-filters-algorithm-sequential-monte-carlo-method-implementation?page=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();
}
你什么时候能在现实世界中使用这个过滤器？你能对具有分析解的问题运行测试吗？如果你已经正确地实现了它，你将得到相同的数字。用错误的实现得出正确结果的可能性非常小！ - Alessandro Teruzzi
@AlessandroTeruzzi 这只是一个简单示例的实现，示例可以在这里找到。我实现它只是为了更好地理解该算法提供的粒子滤波器概念，但我不确定我是否正确地实现了它，因为我并没有很好地理解该算法。我不知道如何测试它是否有效，因为该算法及其输出对我来说仍然不是很清晰（即使该算法似乎非常简单）。 - shn
我对通用算法的第一个建议是：不要试图实现你没有完全理解的东西。先理解，再实现。否则你将无法确定出了什么问题。 - Jorge Leitao
1
@J.C.Leitão，您说得对，但不幸的是，粒子滤波器就是一个例子，除非您尝试实现它，否则您无法真正理解它。 - shn
2个回答
21

这个在线课程非常易懂，讲解粒子滤波器十分清晰。

它被称为“编程机器人车”，介绍了三种本地化方法：蒙特卡罗本地化、卡尔曼滤波器和粒子滤波器。

这门课程完全免费（已经结束了，所以您不能积极参与，但仍然可以观看讲座），由斯坦福教授授课。 "课程" 对我来说令人惊讶的容易跟进，并配有一些小练习--有些仅仅是逻辑上的，但很多是编程方面的。此外还有家庭作业，你可以试着完成。

实际上，它要求你使用Python编写所有滤波器的代码，并为你测试。通过这门课程，你应该能在Python中完成所有3个滤波器的实现。

强烈建议去看看。

好的，我会去查看一下。但是既然我只对粒子滤波器（也许还有其他两种滤波器）感兴趣，那么我是否需要从头开始检查所有单元的所有课程？ - shn
如果你对所有的过滤器都感兴趣，那么你一定要查看整个课程。但即使你不是，我建议你浏览其他课程——它将为你提供非常好的本地化方法概述，并且你将更容易理解每个过滤器最适合的位置。我认为一节课通常可以在4小时左右完成——如果你想轻松一点，可以分1或2天完成。;) - penelope
顺便说一句，我不会用粒子滤波器用于机器人技术方面（如定位等），而是会用于与在线时序数据聚类相关的另一个问题。 - shn
不太喜欢任何基于在线的东西 :( 但我认为它可以帮助你，因为解释非常易懂，而且不太关注机器人的用途。顺便说一句，应该尽量减少闲聊评论 ;) - penelope
1
旧链接失效了，现在更新为：https://www.udacity.com/course/artificial-intelligence-for-robotics--cs373 - Bill Kotsias
3
它们与我的预期有些不同。是否有方法实现 www.cs.ubc.ca/~arnaud/doucet_defreitas_gordon_smcbookintro.ps 第9页中介绍的简单方法？ - shn
当然，有办法 :) 那么链接与您的期望有何不同？您卡在哪里了？我认为第11页上的算法相当简单明了。 - Throwback1986
我想从头开始实现粒子滤波器，以便理解它。我已编辑了我的第一篇帖子，添加了一个简单的C++实现，该实现解释了这里的一个简单示例：http://www.aiqus.com/questions/39942/very-simple-particle-filters-algorithm-sequential-monte-carlo-method-implementation?page =1＃39950您能否请检查我是否理解得好，或者是否存在误解？ - shn
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接