Implementation of a simple particle filter

Motivation

Until now, I used to use the Kalman filter to predict time-series data, but I decided to write a code to see if I could start working on the particle filter. However, this example does not correspond to the decomposition into seasonal cycles. Also, it is not smoothed. please note that.

References

I think that "the basics of statistical modeling that can be used for prediction" can be implemented if it is the content of this article.

Particle filter overview

A type of state space model. A thing that can express a time series with a system model and an observation model is called a state space model. Therefore, the state space model is a model that can incorporate the state (mechanism) behind the observable data. Among them, the thing described in linear and Gaussian type is the Kalman filter. The disadvantage of the Kalman filter is that it is difficult to respond to sudden changes in value because it assumes a normal distribution for noise. For this reason, extended Kalman filters have been developed, but good approximations cannot be obtained if the true distribution is multimodal. On the other hand, the particle filter approximates the distribution with particles, so it can be used even when the distribution is multimodal. And the approximation is possible by two simple operations, time evolution and resampling of each particle. These two operations correspond to generating a predicted distribution and a filter distribution. You can easily see what the particle filter actually does by looking at the diagram on page 76 of "Basics of Statistical Modeling for Prediction".

Then, I will write the code below.

Reference code

python


# coding: utf-8

from math import log, pow, sqrt
import numpy as np
from scipy.stats import norm
from numpy.random import uniform
from multiprocessing import Pool
import matplotlib.pyplot as plt


class ParticleFilter:
    alpha2 = 0.15 #Dispersion ratio of system noise and observed noise
    sigma2 = pow(2,5) #Dispersion of system noise
    log_likelihood = 0.0 #Log likelihood
    LSM = 0.0 #Square error
    TIME = 1
    PR=8 
    
    def __init__(self, PARTICLES_NUM):
        self.PARTICLES_NUM = PARTICLES_NUM #Number of particles
        self.TEETH_OF_COMB = np.arange(0, 1, float(1.0)/self.PARTICLES_NUM) #Horizontal comb used for resampling (reuse)
        self.weights = np.zeros(self.PARTICLES_NUM) #Unit mass of particles (goodness of fit to observation data)
        self.particles = np.zeros(self.PARTICLES_NUM) #particle
        self.predicted_particles = np.zeros(self.PARTICLES_NUM) #Predicted distribution (particles)
        np.random.seed(555)
        self.predicted_value = []
        self.filtered_value = []
    
    def init_praticles_distribution(self):
        """initialize particles
        x_0|0
        """
        self.particles = norm.rvs(0,1,size=self.PARTICLES_NUM)
        
    def get_system_noise(self):
        """v_t"""
        return norm.rvs(0, self.alpha2*self.sigma2, size=self.PARTICLES_NUM)
        
    def calc_pred_particles(self):
        """calculate system function
        x_t|t-1
        """
        return self.particles + self.get_system_noise()  
        
    def calc_particles_weight(self,y):
        """calculate fitness probabilities between observation value and predicted value
        w_t
        """
        locs = self.calc_pred_particles()
        self.predicted_particles = locs
                  
        self.weights = norm.pdf([y]*self.PARTICLES_NUM, loc=locs,
                                scale=[sqrt(self.sigma2)]*self.PARTICLES_NUM)
                  
    def calc_likelihood(self):
        """alculate likelihood at that point
        p(y_t|y_1:t-1)
        """
        res = sum(self.weights)/self.PARTICLES_NUM
        self.log_likelihood += log(res)
#        return res
      
    def normalize_weights(self):
        """wtilda_t"""
        self.weights = self.weights/sum(self.weights)
      
    def resample(self,y):
        """x_t|t"""
        self.normalize_weights()

        self.memorize_predicted_value()

        # accumulate weight
        cum = np.cumsum(self.weights)
        
        # create roulette pointer 
        base = uniform(0,float(1.0)/self.PARTICLES_NUM)
        pointers = self.TEETH_OF_COMB + base
        
        # select particles
        selected_idx = [np.where(cum>=p)[0][0] for p in pointers]
        """
        pool = Pool(processes=self.PR)
        selected_idx = pool.map(get_slected_particles_idx, ((cum,p) for p in pointers))
        pool.close()
        pool.join()     
        """

#         print "select",selected_idx
        self.particles = self.predicted_particles[selected_idx]                
        self.memorize_filtered_value(selected_idx, y)
        
    
    def memorize_predicted_value(self):
        predicted_value = sum(self.predicted_particles*self.weights)
        self.predicted_value.append(predicted_value)

    def memorize_filtered_value(self, selected_idx, y):
        filtered_value = sum(self.particles*self.weights[selected_idx])/sum(self.weights[selected_idx]) # /sum(self.weights[selected_idx])Added
        self.filtered_value.append(filtered_value)
        self.calculate_LSM(y,filtered_value)

    def calculate_LSM(self,y,filterd_value):
        self.LSM += pow(y-filterd_value,2)

    def ahead(self,y):
        """compute system model and observation model"""
        print 'calculating time at %d' % self.TIME
        self.calc_pred_particles()
        self.calc_particles_weight(y)
        self.calc_likelihood()
        self.resample(y)
        self.TIME += 1

def get_slected_particles_idx((cum,p)):
    """multiprocessing function"""
    try:
        return np.where(cum>=p)[0][0]
    
    except Exception:
        import sys
        import traceback
        sys.stderr.write(traceback.format_exc())    

if __name__=='__main__':
    pf = ParticleFilter(1000)
    pf.init_praticles_distribution()
    
    data = np.hstack((norm.rvs(0,1,size=20),norm.rvs(2,1,size=60),norm.rvs(-1,0.5,size=20)))
    
    for d in data:
        pf.ahead(d)
    print 'log likelihood:', pf.log_likelihood
    print 'LSM:', pf.LSM
    
    rng = range(100)
    plt.plot(rng,data,label=u"training data")
    plt.plot(rng,pf.predicted_value,label=u"predicted data")
    plt.plot(rng,pf.filtered_value,label=u"filtered data")
    plt.xlabel('TIME',fontsize=18)
    plt.ylabel('Value',fontsize=18)    
    plt.legend() 
    plt.show()

Experimental result

The experiment was conducted with a particle number of 5,100,1000. Paste the result.

When the number of particles is 5

result_5.png

When the number of particles is 100

resutl_100.png

When the number of particles is 1000

resutlt_1000.png

comment

The system noise was decisive, so I didn't get very good results. However, you can see that it is able to respond to data jumps, and that it gradually fits as you increase the number of particles. By the way, the squared error was 5 = 366.39 particles, 100 = 32.18 particles, 1000 = 20.45 particles.

Summary

The framework is extremely simple, and I found that it was easy to create something as simple as this time. However, I didn't really understand how amazing it was in this example ... Sumimasen </ l> I would like to create it from actual measurement data instead of fixing the noise distribution, or try to transplant the Kalman filter that I made before to a particle filter, and use it in practice.

We apologize for the inconvenience, but if you make a mistake, we would appreciate it if you could point it out.

Postscript

I forgot to write about the model described in the code above. .. .. The particle filter described in this article is a one-time difference trend model.   x_t=x_{t-1}+v_t, v_t~N(0,α^2σ^2)   y_t=x_t+w_t, w_t~N(0,σ^2) It is represented by. Also, although it is described as a particle filter in the article, what is written in the code is a thing called "Monte Carlo filter" among the particle filters, which is a special system of particle filters. This area is [Statistical science of mathematics and calculation (Statistical science of 21st century)](http://www.amazon.co.jp/gp/product/4130440837/ref=as_li_qf_sp_asin_tl?ie=UTF8&camp=247&creative=1211&creativeASIN=4130440837&linkCode = as2 & tag = shimashimao06-22) <img src = "http://ir-jp.amazon-adsystem.com/e/ir?t=shimashimao06-22&l=as2&o=9&a=4130440837" width="1" height=" 1 "border =" 0 "alt =" "style =" border: none! Important; margin: 0px! Important; "/> will be helpful. In addition, the expected value is used for the numerical value in the figure (I think it is not necessary to explain because it uses the product sum of the value and the normalized weight ...). It does not mean that you have to use the expected value, so I hope you can decide by looking at the distribution. That was the postscript.

Addendum 2

The code for the calculation part of the expected value for generating the filter distribution was incorrect. Specifically, I forgot to divide by the sum of the weights. We apologize for the correction m (__) m

Recommended Posts

Implementation of a simple particle filter
A simple sample of pivot_table.
Python implementation of particle filters
A simple Python implementation of the k-nearest neighbor method (k-NN)
Python implementation of self-organizing particle filters
Implementation of a two-layer neural network 2
Explanation and implementation of simple perceptron
Pandas: A very simple example of DataFrame.rolling ()
A simple example of how to use ArgumentParser
Proposal of Kuwahara filter as a photographic expression
Simple implementation example of one kind of data augmentation
[Python] Implementation of clustering using a mixed Gaussian model
Implementation example of simple LISP processing system (Python version)
A python implementation of the Bayesian linear regression class
A very simple example of an ortoolpy optimization problem
Implementation of Fibonacci sequence
Super simple: A collection of shells that output dates
Implementation of a convolutional neural network using only Numpy
A note of trying a simple MCMC tutorial on PyMC3
A reminder about the implementation of recommendations in Python
Implementation of VGG16 using Keras created without using a trained model
[With simple explanation] Scratch implementation of deep Boltzmann machine with Python ②
[With simple explanation] Scratch implementation of deep Boltzmann machine with Python ①
Quantum computer implementation of quantum walk 2
Implementation of TF-IDF using gensim
Implementation of MathJax on Sphinx
A memorandum of kernel compilation
Create a (simple) REST server
Explanation and implementation of SocialFoceModel
Operation of filter (None, list)
A small memorandum of openpyxl
Implementation of game theory-Prisoner's dilemma-
Implementation of independent component analysis
Simple FPS measurement of python
Simple simulation of virus infection
Quantum computer implementation of quantum walk 3
Create a simple textlint server
Filter the output of tracemalloc
A brief summary of Linux
Implementation of quicksort in Python
Quantum computer implementation of quantum walk 1
Deep reinforcement learning 2 Implementation of reinforcement learning
Implementation of Scale-space for SIFT
A memorandum of using eigen3
A simple data analysis of Bitcoin provided by CoinMetrics in Python
Evaluate the performance of a simple regression model using LeaveOneOut cross-validation
Try to make a blackjack strategy by reinforcement learning ((1) Implementation of blackjack)
[Python] A simple function to find the center coordinates of a circle