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

自己組織化型状態空間モデルのコード(未完成)

More than 5 years have passed since last update.

動機

変化点検知などで精度の高い時系列予測モデルが必要になったため、コード作成を試みてみた。結果を先に書きます。「スゲー精度のいいのが出来たゼ!(ドヤァッ」ってしたかったのですが、挙動がおかしい中途半端なコードが出来上がりました...orz
今後、修正でき次第上げ直します。。

参考文献

「初期分布探索付き自己組織化状態空間モデルによる金融時系列解析の最前線:t分布付き確率的ボラティリティ変動モデルへの応用」はモンテカルロフィルタの課題がまとめられており、更に最適な初期分布の探索手法も提案されていて、非常に面白い資料でした。
これを参考に初期分布の最適化まで含めたモノを作成するつもりだったのですが、それとは程遠いモノが出来上がってしまいました。。

モンテカルロフィルタの問題点

モンテカルロフィルタは、非定常非ガウス型のモデルが構築出来る非常に柔軟な状態推定アルゴリズムです。しかし、以下の事が問題点として挙げられています。

モンテカルロフィルタにおけるパラメタ推定については(1)パラメタに対する尤度関数の推定量がモンテカルロ法による誤差を含むこと、(2)尤度関数の微分が算出困難であること(初期分布探索付き自己組織化状態空間モデルによる金融時系列解析の最前線 P.4)

これを解決するために、Nelder-Mead法(直線探索型の非線形最適化手法)という最適化手法を使用した手法が提案されています。が、この手法では収束条件を大幅に緩和しなければ収束せずパラメタ推定量の分散が大きくなる、という問題点が残ってしまいます。
一方、自己組織化状態空間モデルは、パラメタを状態空間に繰り込んでBayes的推定を行う方法です。
「初期分布探索付き自己組織化状態空間モデルによる金融時系列解析の最前線」ではNM法を用いた、初期分布の構成手法(最尤推定方法)が提案されています。
モンテカルロフィルタについては、簡単な粒子フィルタの実装を参考にして頂きたく。

モデル

以下のモデルのコードを作成してみました。
$x_t=Fx_{t-1}+Gv_t$
$y_t=Hx_{t}+w_t, w_t~C(0,σt^2)$
$v
{t,t_t}~C(0,τt^2)$
$v
{t,{τt^2}}~C(0,ν^2)$
$v
{t,{σt^2}}~C(0,ξ^2)$
$x_t=[t_t,log
{10}τt^2,log{10}σt^2]^T$
$v
{t}=[v_{t,t_t},v_{t,{τt^2}},v{t,{σ_t^2}}]^T$

「知識発見と自己組織型の統計モデル」に記載がありますが、たとえば$σ^2$を推定する場合、正値性を保証する必要があるため$log_{10}σ^2$の推定を行っています。

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

参考コード

# coding: utf-8

from math import log, pow, sqrt
import numpy as np
from scipy.stats import norm, cauchy
from numpy.random import uniform, multivariate_normal
from multiprocessing import Pool
import matplotlib.pyplot as plt


class ParticleFilter:
    nu = 0.01 # システムノイズのスケール超々パラメタ
    xi = 0.03 # 観測ノイズのスケール超々パラメタ
    log_likelihood = 0.0 # 対数尤度
    TIME = 1
    PR=8 # unmber of processing

    def __init__(self, PARTICLES_NUM, k=1, ydim=1, pdim=2):
        self.PARTICLES_NUM = PARTICLES_NUM # 粒子の数
        self.TEETH_OF_COMB = np.arange(0, 1, float(1.0)/self.PARTICLES_NUM)
        self.weights = np.zeros((ydim, self.PARTICLES_NUM))
        self.particles = np.zeros((k*ydim+pdim ,self.PARTICLES_NUM))
        self.predicted_particles = np.zeros((k*ydim+pdim , self.PARTICLES_NUM))
        np.random.seed(555)
        self.predicted_value = []
        self.filtered_value = []
        self.LSM = np.zeros(ydim) # 2乗誤差

        self.F, self.G, self.H = self.FGHset(k, ydim, pdim)
        self.k = k
        self.ydim = ydim
        self.pdim = pdim

    def init_praticles_distribution(self, P, r):
        """initialize particles
        x_0|0
        tau_0|0
        sigma_0|0
        """
        data_particles = multivariate_normal([1]*self.ydim*self.k,
                            np.eye(self.ydim*self.k)*10, self.PARTICLES_NUM).T
        param_particles = np.zeros((self.pdim, self.PARTICLES_NUM))
        for i in xrange(self.pdim):
            param_particles[i,:] = uniform(P-r, P+r, self.PARTICLES_NUM)
        self.particles = np.vstack((data_particles, param_particles))

    def get_system_noise(self):
        """v_t vector"""
        data_noise = cauchy.rvs(loc=[0]*self.PARTICLES_NUM, scale=np.power(10,self.particles[self.ydim])) #size=self.PARTICLES_NUM
        data_noise[data_noise==float("-inf")] = -1e308
        data_noise[data_noise==float("inf")] = 1e308
        parameta_noise_sys = cauchy.rvs(loc=0, scale=self.xi, size=self.PARTICLES_NUM) # noise of hyper parameter in system model
        parameta_noise_ob = cauchy.rvs(loc=0, scale=self.nu, size=self.PARTICLES_NUM) # noise of hyper parameter in observation model

        #print np.vstack((data_noise, parameta_noise_sys, parameta_noise_ob))
        return np.vstack((data_noise, parameta_noise_sys, parameta_noise_ob))

    def calc_pred_particles(self):
        """calculate system function
        x_t|t-1 = F*x_t-1 + Gv_t
        """
        return self.F.dot(self.particles) + self.G.dot(self.get_system_noise()) # linear non-Gaussian  

    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
        scale=np.power(10,locs[-1])
        scale[scale==0] = 1e-308

        self.weights = cauchy.pdf( np.array([y]*self.PARTICLES_NUM) - self.H.dot(locs), loc=[0]*self.PARTICLES_NUM,
                                scale=scale ).flatten()

    def calc_likelihood(self):
        """calculate likelihood at that point
        p(y_t|y_1:t-1)
        """
        res = np.sum(self.weights)/self.PARTICLES_NUM
        #print res
        self.log_likelihood += log(res)

    def normalize_weights(self):
        """wtilda_t"""
        self.weights = self.weights/np.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()     
        """

        self.particles = self.predicted_particles[:,selected_idx]
        self.memorize_filtered_value(selected_idx, y)

    def memorize_predicted_value(self):
        predicted_value = np.sum(self.predicted_particles*self.weights, axis=1)[0]
        self.predicted_value.append(predicted_value)

    def memorize_filtered_value(self, selected_idx, y):
        filtered_value = np.sum(self.particles*self.weights[selected_idx] , axis=1) \
                            /np.sum(self.weights[selected_idx])
        self.filtered_value.append(filtered_value[0])
        self.calculate_LSM(y,filtered_value[:self.ydim])

    def calculate_LSM(self,y,filterd_value):
        self.LSM += pow(y-filterd_value,2)

    def forward(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 FGHset(self, k, vn_y, n_ss_parameters):
        """
        vn_y:input dimension
        n_ss_parameters:number of hyper parameter 
        k:difference
        """
        G_upper_block = np.zeros((k*vn_y, vn_y+n_ss_parameters))
        G_lower_block = np.zeros((n_ss_parameters, vn_y+n_ss_parameters))
        G_lower_block[-n_ss_parameters:, -n_ss_parameters:] = np.eye(n_ss_parameters)
        G_upper_block[:vn_y, :vn_y] = np.eye(vn_y)
        G = np.vstack( (G_upper_block, G_lower_block) )

        H = np.hstack( (np.eye(vn_y), 
                            np.zeros((vn_y, vn_y*(k-1)+n_ss_parameters))
                        ) )

        F_upper_block = np.zeros((k*vn_y, k*vn_y+n_ss_parameters))
        F_lower_block = np.zeros((n_ss_parameters, k*vn_y+n_ss_parameters))
        F_lower_block[-n_ss_parameters:, -n_ss_parameters:] = np.eye(n_ss_parameters)
        if k==1:
            F_upper_block[:vn_y, :vn_y] = np.eye(vn_y)
        elif k==2:
            F_upper_block[:vn_y, :vn_y] = np.eye(vn_y)*2
            F_upper_block[:vn_y, vn_y:k*vn_y] = np.eye(vn_y)*-1
            F_upper_block[vn_y:k*vn_y, :vn_y] = np.eye(vn_y)
        F = np.vstack((F_upper_block, F_lower_block))

        return F, G, H

def get_slected_particles_idx((cum,p)):
    """multiprocessing function"""
    try:
        return np.where(cum>=p)[0][0]

    except Exception, e:
        import sys
        import traceback
        sys.stderr.write(traceback.format_exc())    

if __name__=='__main__':
    n_particle = 10000
    pf = ParticleFilter(n_particle, k=1, ydim=1, pdim=2)
    pf.init_praticles_distribution(0, 8) # P, r

    data = np.hstack((norm.rvs(0,1,size=20),norm.rvs(10,1,size=60),norm.rvs(-30,0.5,size=20)))

    for d in data:
        pf.forward(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()

実験結果

結果の図を載せるまでもない欠陥だらけのコードです...

コメント

今回は中途半端な状態で時間切れとなってしまい残念です。今後まともな結果が出るように修正して、再度上げ直します。
(参考にならない中途半端なモノを公開して申し訳無いのですが、今回は自分用のメモも兼ねてupしました。)

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

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
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