[PYTHON] Implementierung eines einfachen Partikelfilters
Motivation
Bisher habe ich den Kalman-Filter verwendet, um Zeitreihendaten vorherzusagen, aber ich habe beschlossen, einen Code zu schreiben, um auch den Partikelfilter auszuprobieren. Dieses Beispiel entspricht jedoch nicht der Zerlegung in saisonale Zyklen. Auch wird es nicht geglättet. bitte beachte, dass.
Verweise
- Grundlagen der statistischen Modellierung für die Vorhersage - von der Einführung in die Bayes'sche Statistik bis zur Anwendung 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! Wichtig; Rand: 0px! Wichtig; "/>
- [Einführung in die Zeitreihenanalyse](http://www.amazon.co.jp/gp/product/4000054554/ref=as_li_tf_tl?ie=UTF8&camp=247&creative=1211&creativeASIN=4000054554&linkCode=as2&tag=shimashimao0g6g = "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 = "Rahmen: keine! wichtig; Rand: 0px! wichtig;" />
Ich denke, dass "die Grundlagen der statistischen Modellierung, die für die Vorhersage verwendet werden können", implementiert werden können, wenn es sich um den Inhalt dieses Artikels handelt.
Übersicht über Partikelfilter
Eine Art Zustandsraummodell. Eine Sache, die die Zeitreihen mit einem Systemmodell und einem Beobachtungsmodell ausdrücken kann, wird als Zustandsraummodell bezeichnet. Daher ist das Zustandsraummodell ein Modell, das den Zustand (Mechanismus) hinter den beobachtbaren Daten einbeziehen kann. Unter diesen ist das im linearen und Gaußschen Typ beschriebene Ding der Kalman-Filter. Der Nachteil des Kalman-Filters besteht darin, dass es schwierig ist, auf plötzliche Wertänderungen zu reagieren, da es eine Normalverteilung für Rauschen annimmt. Aus diesem Grund wurden erweiterte Kalman-Filter entwickelt, aber gute Näherungen können nicht erhalten werden, wenn die wahre Verteilung multimodal ist. Andererseits approximiert der Partikelfilter die Verteilung mit Partikeln, so dass er auch dann verwendet werden kann, wenn die Verteilung multimodal ist. Die Annäherung ist durch zwei einfache Operationen möglich: Zeitentwicklung und Resampling jedes Partikels. Diese beiden Operationen entsprechen dem Erzeugen einer vorhergesagten Verteilung und einer Filterverteilung. Sie können leicht sehen, was der Partikelfilter tatsächlich tut, indem Sie sich das Diagramm auf Seite 76 von "Grundlagen der statistischen Modellierung für die Vorhersage" ansehen.
Dann werde ich den folgenden Code schreiben.
Referenzcode
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 #Dispersionsverhältnis von Systemrauschen und beobachtetem Rauschen
sigma2 = pow(2,5) #Streuung des Systemrauschens
log_likelihood = 0.0 #Protokollwahrscheinlichkeit
LSM = 0.0 #Quadratischer Fehler
TIME = 1
PR=8
def __init__(self, PARTICLES_NUM):
self.PARTICLES_NUM = PARTICLES_NUM #Anzahl der Partikel
self.TEETH_OF_COMB = np.arange(0, 1, float(1.0)/self.PARTICLES_NUM) #Horizontaler Kamm zum erneuten Abtasten (Wiederverwendung)
self.weights = np.zeros(self.PARTICLES_NUM) #Einheitsmasse der Partikel (Eignung für Beobachtungsdaten)
self.particles = np.zeros(self.PARTICLES_NUM) #Partikel
self.predicted_particles = np.zeros(self.PARTICLES_NUM) #Voraussichtliche Verteilung (Partikel)
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])Hinzugefügt
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()
Versuchsergebnis
Das Experiment wurde mit 5.100.1000 Partikeln durchgeführt. Fügen Sie das Ergebnis ein.
Wenn die Anzahl der Partikel 5 beträgt
Wenn die Anzahl der Partikel 100 beträgt
Wenn die Anzahl der Partikel 1000 beträgt
Kommentar
Das Systemrauschen war entscheidend, so dass ich keine sehr guten Ergebnisse erzielte. Sie können jedoch sehen, dass es auf Datensprünge reagieren kann und dass es allmählich passt, wenn Sie die Anzahl der Partikel erhöhen. Der quadratische Fehler betrug übrigens 5 = 366,39 Partikel, 100 = 32,18 Partikel, 1000 = 20,45 Partikel.
Zusammenfassung
Das Framework ist extrem einfach und ich fand, dass es einfach war, etwas so Einfaches wie dieses Mal zu erstellen.
Ich habe jedoch nicht wirklich verstanden, wie erstaunlich es in diesem Beispiel war ...
Wir entschuldigen uns für die Unannehmlichkeiten, aber wenn Sie einen Fehler machen, würden wir uns freuen, wenn Sie darauf hinweisen könnten.
Nachtrag
Ich habe vergessen, über das im obigen Code beschriebene Modell zu schreiben. .. ..
Der in diesem Artikel beschriebene Partikelfilter ist ein einmaliges Differenztrendmodell.
Nachtrag 2
Der Code für den Berechnungsteil des erwarteten Werts zum Generieren der Filterverteilung war falsch. Insbesondere habe ich vergessen, durch die Summe der Gewichte zu teilen. Wir entschuldigen uns für die Korrektur m (__) m
Recommended Posts