Edited at

簡単な粒子フィルタの実装

More than 5 years have passed since last update.


動機

今まで時系列データの予測にはカルマンフィルタを使う事が多かったのですが、そろそろ粒子フィルタにも手を出してみるか、ということでコードを書いてみました。

ただ今回の例は、季節周期などへの分解には対応していません。また、平滑化も行っていません。ご了承下さい。


参考文献

「予測にいかす統計モデリングの基本」だけでも、この記事の内容であれば実装可能だと思います。


粒子フィルタ概要

状態空間モデルの一種です。時系列をシステムモデルと観測モデルで表現出来るようなモノを状態空間モデルと呼びます。したがって、状態空間モデルは観測可能なデータの裏側にある状態(メカニズム)を組み込む事が出来るモデルになります。

その中で線形、ガウス型で記述されたモノがカルマンフィルタです。カルマンフィルタの欠点としては、ノイズに正規分布を仮定しているため、急激な値の変化への対応が難しい点が挙げられます。そのため、拡張カルマンフィルタなどが開発されているが、真の分布が多峰性の場合は良い近似が得られません。

一方、粒子フィルタはというと、分布を粒子で近似するため分布が多峰性の場合にも対応が可能です。そして、その近似は各粒子の時間発展とリサンプリングという2つの簡単な操作によって可能となります。

この2つの操作は、予測分布とフィルタ分布を発生させる事に対応します。

粒子フィルタが実際何をやっているかは、「予測にいかす統計モデリングの基本」の76ページの図をご覧いただくのがわかりやすいと思います。

では、以下にコードを記載します。


参考コード

# 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 # システムノイズと観測ノイズの分散比
sigma2 = pow(2,5) # システムノイズの分散
log_likelihood = 0.0 # 対数尤度
LSM = 0.0 # 2乗誤差
TIME = 1
PR=8

def __init__(self, PARTICLES_NUM):
self.PARTICLES_NUM = PARTICLES_NUM # 粒子の数
self.TEETH_OF_COMB = np.arange(0, 1, float(1.0)/self.PARTICLES_NUM) # リサンプリングに使用する横櫛(使い回す)
self.weights = np.zeros(self.PARTICLES_NUM) # 粒子の単位質量(観測データへの適合度)
self.particles = np.zeros(self.PARTICLES_NUM) # 粒子
self.predicted_particles = np.zeros(self.PARTICLES_NUM) # 予測分布(粒子)
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])を追記
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()


実験結果

粒子数5,100,1000で実験を行いました。その結果を貼り付けます。


粒子数5の場合

result_5.png


粒子数100の場合

resutl_100.png


粒子数1000の場合

resutlt_1000.png


コメント

システムノイズが決め打ちだったので、あまり良い結果は出ませんでした。

ですが、データのジャンプへの対応は出来ているのが確認出来るのと、粒子を増やしていくと徐々にフィットしていくのがわかると思います。

ちなみに2乗誤差は、粒子数5=366.39, 粒子数100=32.18, 粒子数1000=20.45 でした。


まとめ

枠組みはものすごく単純で、今回のような簡単なものであれば簡単に作れる事がわかりました。

ただ、今回の例ではあまりスゴさがわからなかったですね...スミマセン

ノイズの分布を決め打ちするのではなく実測データから作成したり、以前作っていたカルマンフィルタを粒子フィルタに移植したりしてみたり、実務で使用してみたいと思います。

お手数ですが間違いありましたら、ご指摘いただけますと助かります。


追記

上のコードで記載したモデルについて書くのを忘れていました。。。

この記事に書いた粒子フィルタは一回差分のトレンドモデルで、

  $x_t=x_{t-1}+v_t, v_t~N(0,α^2σ^2)$

  $y_t=x_t+w_t, w_t~N(0,σ^2)$

で表されます。

それと、記事内で粒子フィルタと記載していますが、コード内に書かれているのは粒子フィルタの中でも「モンテカルロフィルタ」と呼ばれるモノで、粒子フィルタの特殊系です。

この辺りは数理・計算の統計科学 (21世紀の統計科学)が参考になります。

また、図中の数値には期待値を使用しています(値と正規化された重みの積和を使っているので説明するまでもないと思うのですが...)。

期待値を使用しなければいけないということではありませんので、分布を見て決めて頂ければと思います。

以上追記でした。


追記2

フィルタ分布の生成の期待値の計算部分のコードが間違っておりました。

具体的には重みの総和で除するのを忘れていました。

訂正とお詫び申し上げますm(__)m