Top Tutorial 3

Tutorial 4: How to create an application using a Particle Filter

In this example, we use a Gaussian random walk (with a constant shift of 1) as a toy state model: p(x_{i+1}|x_i)=N(x_i+1, sigma), where N(x, s) is a Gaussian distribution with mean x and std-deviation s. This means that, given the current state x_i, the next state x_{i+1} is likely to be x_i+1. This simple linear model is used in the state transition distribution as well as in the importance distribution. We use also a Gaussian for the observation model, i.e., p(y_i|x_i)=N(y_i-x_i, sigma2).

Prior distribution

We start by defining the prior distribution p(x)=N(0, sigma3). This method is used to initially set all states, including the previous and the current states, according to the prior distribution.
#include <sl/monte_carlo.h>

template< typename T >
struct prior_sampler {
    std::normal_distribution< T > _sampler;   
    prior_sampler(const T &stddev = 1) : _sampler(0, stddev) {}

    template< typename T_SIR >
    void sample(T_SIR &sir) {
        for (int i1 = -sir.sample_sets().size() + 1; i1 <= 0; i1++) {
            dotimes (i0, sir.sample_set().size()) {
                sir.sample_set_at(i1) [i0].position() [0] = _sampler(sir.rn_engine());
            }
        }
    }
};

Observation distribution

We now define the likelihood, e.g., the observation distribution p(y_i|x_i)=N(y_i-x_i, sigma2).
template< typename T >
struct observation_dist_sampler {
    std::normal_distribution< T > _sampler;
    
    observation_dist_sampler(const T &stddev) : _sampler(0, stddev) {}
    template< typename T_SIR >
    typename T_SIR::value_type pdf(T_SIR &sir, unsigned i0) {
        return sl::pdf(_sampler, sir.observation() [0] - sir.sample_set() [i0].position() [0]);
    }
};

Transition distribution

The following code implements the transition distribution p(x_{i+1}|x_i)=N(x_{i+1}-1-x_i, sigma1).
template< typename T >
struct transition_dist_sampler {
    std::normal_distribution< T > _sampler;
    
    transition_dist_sampler(const T &stddev) : _sampler(0, stddev) {}
    template< typename T_SIR >
    typename T_SIR::value_type pdf(T_SIR &sir, unsigned i0) {
        return sl::pdf(_sampler, sir.sample_set() [i0].position()[0] - sir.prev_sample_set() [i0].position()[0] - 1);
    }
};

Importance distribution

The following code implements the importance distribution p(x_{i+1}|x_i)=N(x_{i+1}-1-x_i, sigma), which is used to generate the current samples. Ideally, it should be very close to the current posterior distribution.
template< typename T >
struct importance_dist_sampler {
    std::normal_distribution< T > _sampler;
    
    importance_dist_sampler(const T &stddev) : _sampler(0,stddev) {}
    template< typename T_SIR >
    void sample(T_SIR &sir) {
        dotimes (i0, sir.sample_set().size()) {
            sir.sample_set() [i0].position() [0] = sir.prev_sample_set() [i0].position()[0] + 1 + _sampler(sir.rn_engine());
        }
    }
    template< typename T_SIR >
    typename T_SIR::value_type pdf(T_SIR &sir, unsigned i0) {
        return sl::pdf(_sampler, sir.sample_set() [i0].position() [0] - sir.prev_sample_set() [i0].position()[0] - 1);
    }
};

Setting up and running the particle filter

int main() {
    std::mt19937 engine(time(NULL));

    prior_sampler< double > prior_sampler(0.05);
    observation_dist_sampler< double > observation_op(0.025);
    transition_dist_sampler< double > state_transition_op(0.025);
    importance_dist_sampler< double > importance_op(0.05);
    
    sl::generic_SIR< observation_dist_sampler< double >,
                     transition_dist_sampler< double >,
                     importance_dist_sampler< double > > p_filter(engine,
                                                               observation_op,
                                                               state_transition_op,
                                                               importance_op);
    int nsamples = 50, dim = 1, history_size = 2;
    p_filter.init_construction(nsamples, dim, history_size);
    p_filter.posterior_processor()._interactive = true;
    p_filter.init(prior_sampler);
        
    std::normal_distribution< double > sampler(0, 0.025);
    sl::vector< double > obs_v(1, 0);
    dotimes (i0, 10) {
        obs_v [0] += 1 + sampler(engine);
        std::cout << "sample: " << i0 << ",  obs = " << obs_v [0] << std::endl;
        p_filter(obs_v);
    }
    std::cout << "done!\n";

    return 1;
}
This code is from the file 'particle-filter-example.cpp' in the mcmll/examples directory. Switch to the mcmll/examples directory and compile all examples by the command 'scons', or alternatively, compile it with:
  g++ particle-filter-example.cpp  -I path-to-mcmll/lib -std=c++0x -D__SL_USE_CXX0X__ -o particle-filter-example
where path-to-mcmll is the installation path of the MCMLL library (for example: /home/me/cpp/mcmll). Run example with:
  ./particle-filter-example