[PYTHON] Mise en place d'un filtre à particules simple
Motivation
Jusqu'à présent, j'utilisais le filtre de Kalman pour prédire les données de séries chronologiques, mais j'ai décidé d'écrire un code pour voir si je pouvais commencer à travailler sur le filtre à particules. Cependant, cet exemple ne correspond pas à une décomposition en cycles saisonniers. De plus, il n'est pas lissé. veuillez noter que.
Les références
- [Bases de la modélisation statistique pour la prédiction - de l'introduction aux statistiques bayésiennes à l'application](http://www.amazon.co.jp/gp/product/4061557955/ref=as_li_qf_sp_asin_tl?ie=UTF8&camp=247&creative=1211&creativeASIN=4061557955&linkCreativeASIN=4061557955 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 à l'analyse des séries temporelles "http://ir-jp.amazon-adsystem.com/e/ir?t=shimashimao06-22&l=as2&o=9&a=4000054554" width = "1" height = "1" border = "0" alt = "" style = "border: none! important; margin: 0px! important;" />
Je pense que "les bases de la modélisation statistique qui peuvent être utilisées pour la prédiction" peuvent être implémentées si c'est le contenu de cet article.
Vue d'ensemble du filtre à particules
Un type de modèle d'espace d'états. Une chose qui peut exprimer la série temporelle avec un modèle système et un modèle d'observation est appelée un modèle d'espace d'états. Par conséquent, le modèle d'espace d'états est un modèle qui peut incorporer l'état (mécanisme) derrière les données observables. Parmi eux, la chose décrite en type linéaire et gaussien est le filtre de Kalman. L'inconvénient du filtre de Kalman est qu'il est difficile de répondre à des changements brusques de valeur car il suppose une distribution normale du bruit. Pour cette raison, des filtres de Kalman étendus ont été développés, mais de bonnes approximations ne peuvent être obtenues si la vraie distribution est multimodale. D'autre part, le filtre à particules se rapproche de la distribution avec des particules, de sorte qu'il peut être utilisé même lorsque la distribution est multimodale. Et l'approximation est possible par deux opérations simples: évolution temporelle et rééchantillonnage de chaque particule. Ces deux opérations correspondent à la génération d'une distribution prédite et d'une distribution de filtre. Vous pouvez facilement voir ce que fait réellement le filtre à particules en regardant le diagramme à la page 76 de «Bases de la modélisation statistique pour la prédiction».
Ensuite, j'écrirai le code ci-dessous.
Code de référence
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 #Rapport de dispersion du bruit du système et du bruit observé
sigma2 = pow(2,5) #Dispersion du bruit du système
log_likelihood = 0.0 #Probabilité du journal
LSM = 0.0 #Erreur carrée
TIME = 1
PR=8
def __init__(self, PARTICLES_NUM):
self.PARTICLES_NUM = PARTICLES_NUM #Nombre de particules
self.TEETH_OF_COMB = np.arange(0, 1, float(1.0)/self.PARTICLES_NUM) #Peigne horizontal utilisé pour le rééchantillonnage (réutilisation)
self.weights = np.zeros(self.PARTICLES_NUM) #Masse unitaire des particules (aptitude aux données d'observation)
self.particles = np.zeros(self.PARTICLES_NUM) #particule
self.predicted_particles = np.zeros(self.PARTICLES_NUM) #Distribution prédite (particule)
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])Ajoutée
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()
Résultat expérimental
L'expérience a été menée avec 5 100 000 particules. Collez le résultat.
Lorsque le nombre de particules est de 5
Lorsque le nombre de particules est de 100
Lorsque le nombre de particules est de 1000
commentaire
Le bruit du système a été déterminant, je n'ai donc pas obtenu de très bons résultats. Cependant, vous pouvez voir qu'il est capable de répondre aux sauts de données et qu'il s'adapte progressivement à mesure que vous augmentez le nombre de particules. À propos, l'erreur quadratique était de 5 = 366,39 particules, 100 = 32,18 particules, 1000 = 20,45 particules.
Résumé
Le cadre est extrêmement simple et j'ai trouvé qu'il était facile de créer quelque chose d'aussi simple que cette fois.
Cependant, je n'ai pas vraiment compris à quel point c'était incroyable dans cet exemple ...
Nous vous prions de nous excuser pour la gêne occasionnée, mais si vous faites une erreur, veuillez nous en informer.
Postscript
J'ai oublié d'écrire sur le modèle décrit dans le code ci-dessus. .. ..
Le filtre à particules décrit dans cet article est un modèle de tendance à différence unique.
Addendum 2
Le code de la partie calcul de la valeur attendue pour générer la distribution du filtre était incorrect. Plus précisément, j'ai oublié de diviser par la somme des poids. Nous nous excusons pour la correction m (__) m
Recommended Posts