[PYTHON] 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
- [Basics of Statistical Modeling for Prediction-From Introduction to Bayesian Statistics to Applications](http://www.amazon.co.jp/gp/product/4061557955/ref=as_li_qf_sp_asin_tl?ie=UTF8&camp=247&creative=1211&creativeASIN=4061557955&linkCode= as2 & tag = shimashimao06-22) <img src = "http://ir-jp.amazon-adsystem.com/e/ir?t=shimashimao06-22&l=as2&o=9&a=4061557955" width="1" height="1 "border =" 0 "alt =" "style =" border: none! Important; margin: 0px! Important; "/>
- Introduction to Time Series Analysis
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
When the number of particles is 100
When the number of particles is 1000
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 ...
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.
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