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

- In line 1, we include the header file containing the classes required by the SMC framework. We assume that all following code is inside one singel file, so that this header file is valid also for all the remaining code
- In line 5, we create a 1-D Gaussian distribution object
- In line 6, the constructor is provided with a std-deviation for the Gaussian distribution
- In line 8-15, the
**required**'sample' method is defined, which will be called by the SMC-framework - In the following lines, all samples, including the samples which make up the history, are drawn from the prior distribution. The method 'sir.sample_sets()' returns an array of all available posterior approximations (sample sets) over time (some preceding and the current one), the method 'sir.sample_set()' returns an array to the samples which make up the current posterior approximation

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]); } };

- In line 3, we create a 1-D Gaussian distribution object
- In line 5, the constructor is provided with a std-deviation for the Gaussian distribution
- In line 6-9, the
**required**'pdf' method for calculating the**p**robability**d**ensity**f**unction is defined, which will be called by the SMC-framework - In line 8, the pdf p(y_i|x_i)=N(y_i-x_i, sigma2) is calculated using sample no. i0 and returned

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); } };

- In line 3, we create a 1-D Gaussian distribution object
- In line 5, the constructor is provided with a std-deviation for the Gaussian distribution
- In line 6-9, the
**required**'pdf' method for calculating the**p**robability**d**ensity**f**unction is defined, which will be called by the SMC-framework - In line 8, the pdf p(x_{i+1}|x_i)=N(x_{i+1}-1-x_i, sigma1) is calculated using sample no. i0 and returned

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); } };

- In line 3, we create a 1-D Gaussian distribution object
- In line 5, the constructor is provided with a std-deviation for the Gaussian distribution
- In line 6-11, the
**required**'sample' method is defined, which will be called by the SMC-framework. This method generates samples according to the distribution p(x_{i+1}|x_i)=N(x_{i+1}-1-x_i, sigma) - In line 8-10, new samples are drawn according to the importance distribution
- In line 12-15, the
**required**'pdf' method for calculating the**p**robability**d**ensity**f**unction is defined, which will be called by the SMC-framework - In line 14, the pdf pi(x_{i+1}|x_i)=N(x_{i+1}-1-x_i, sigma) is calculated using sample no. i0 and returned

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; }

- Line 2: create the mersenne twister random number engine
- Line 4-7: create instances of the prior, observation, transition and importance distributions
- Line 9: create the particle filter 'p_filter' based on the Sequential Importance Resampling (SIR) method, provide the random number engine and instances of all required distributions
- Line 15: define the number of samples to use for the approximation of the posterior distribution nsamples, the state dimension 'dim' and the total size of tracked state history 'history_size'
- Line 16: initialize the filter internal structures which includes the allocation of memory according to the provided settings
- Line 17: enable wait for key after each posterior update
- Line 18: apply the prior sampler to initially populate all previous and current states
- Line 20: create a Gaussian distribution object
- Line 21: create the observation vector 'obs'
- Line 22: start the loop to iteratively generate 10 artificial observations, which are sequentially fed to the particle filter
- Line 23: we apply a simple linear state update
- Line 24: print out the current sample number and the current observation
- Line 25: feed the filter with the current observation and let the filter update the posterior distribution approximation (samples and corresponding weights). At the same time, the expectation E (ML-estimate) and the Maximum a Posteriori (MAP) estimate of the state is printed

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