Particle++

a C++ class template for particle filter
source at github

Particle filters or sequential Monte Carlo methods require sampling a large amount of particles at each iteration. This makes it particularly slow in MATLAB simulation. Therefore C++ implementation is desirable for the sake of speed. This template provides some useful C++ classes for users to simulate particle filters. We try to use the STL as much as possible to provide user with an intuitive and consistant interface. Unlike other particle filter softwares, we leave all the definition work of numerical functions to the user, therefore make the template highly extensible. Due to this design, the template does not depend on any other package or library. Only C++11 compiler is required. The following simple example shows how to use the template.

Example: Stochastic Volatility model

Consider the following stochastic volatility model,

 	X_n = alpha X_{n-1} + V_n
 	Y_n = beta exp(X_n / 2) W_n

where V_n and W_n are independent standand normal variables. In this case, the state pdf is f(x'|x)=mathcal{N}(x';alpha x,1) and the observation pdf is g(y|x)=mathcal{N}(y;0,beta^2 exp(x/2)). We simply use standard Gaussian as proposal distribution, q(x'|x,y)=mathcal{N}(x';x,1). Here we let alpha = 0.91 and beta = 1.

To use the particle filter template, you need to follow the following steps:

  1. Define statetype and obsvtype type

  2. Define f, g, q, q_sam function as ( q_sam is the proposal sampler )

    1. long double f( state_type x_n, state_type x_{n-1} )

    2. long double g( state_type x_n, obsv_type y )

    3. long double q( state_type x_n, state_type x_{n-1}, obsv_type y )

    4. state_type q_sam( state_type x_{n-1}, obsv_type y ) // return value is x_n

  3. construct a particle filter object with the above four functions

The code and the simulation result are shown below.



tracking result of the example 






/// main.cpp
/// this example shows how to use the C++ particle filter template

# include <iostream>
# include <fstream>
# include <cmath>
# include <chrono>
# include <random>

# include "pfilter.h"  // include the template

// initialize the random seed
unsigned seed = std::chrono::system_clock::now().time_since_epoch().count();
std::default_random_engine generator (seed);

const long double PI = 3.14159265359;
const long double alpha = 0.91;
const long double beta = 1.0;

typedef long double statetype;
typedef long double obsvtype;

std::normal_distribution <statetype> distribution(0.0,1.0);

// state pdf
long double f(statetype x1, statetype x2){
    return exp(-0.5*pow((x1-alpha*x2),2));
}

// observation pdf
long double g(statetype x, obsvtype y){
    return 1/exp(x/2)*exp(-0.5*pow(y/beta/exp(x/2),2));
}

// proposal distribution
long double q(statetype x1, statetype x2, obsvtype y){
    return exp(-0.5*pow((x1-alpha*x2),2));
}

// proposal sampling function
statetype q_sam(statetype x, obsvtype y){
    return distribution(generator)+alpha*x;
}


int main(){
	// construct the particle filter object
    pfilter <statetype,obsvtype> A(f,g,q,q_sam);  
	
    std::ifstream in("data_y");   	// data input
    std::ofstream on("data_xhat");	// data output
    in >> A;            // input data
    A.initialize(200);  // initialize with 
						// the number of particles we want to use
    A.iterate(); 	    // run
    on << A; 			  // output data
    return 0;    	    //  we are done
}