Help us understand the problem. What is going on with this article?

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

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

shima_x
分析とか雑用とかやってます
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
No comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
ユーザーは見つかりませんでした