カルマンスムージング(KS) - WHAT' CHA GONNA DO FOR ME?

前回は、カルマンフィルターをやったので、今回はスムージング（平滑化）をやろうと思う。

<a href="http://lofas.hatenablog.com/entry/2015/02/21/144227">カルマンフィルタ - WHAT' CHA GONNA DO FOR ME?</a>

モデルは前回と同様、下記のシンプルなモデルとする。パラメータも同じ値をセットした。

${ \begin{eqnarray} y_t &=& x_t+\sigma_{\epsilon}\epsilon_t, \ \epsilon_t \sim\ \mathcal{N}(0,1), \nonumber \\ x_t &=& \phi x_{t-1}+\sigma_{\eta}\eta_t, \ \eta_t \sim\ \mathcal{N}(0,1). \nonumber \end{eqnarray} }$

前回は明示的に書かなかったが、フィルタリングにおいては、 ${t}$ 期の観測値 ${y_t}$ が得られた下での ${t}$ 期の潜在変数 ${x_t}$ の分布 ${p(x_t|y_t)}$ を求めにいった。

これに対して、スムージングにおいては、 ${T}$ 期までの観測値 ${y_{1:T}}$ がすべて得られたもとでの各 ${t}$ 期の潜在変数 ${x_t}$ の分布 ${p(x_t|y_{1:T})}$ が求められる。

フィルタリングに比べて、より多くの観測値の情報を用いて潜在変数を推定することになるため、一般的に精度はより高くなる。

アルゴリズムとしては、カルマンフィルターで ${1}$ 期から ${T}$ 期までの潜在変数の平均と分散を計算した後、それらの値を使いながら ${T}$ 期から ${1}$ 期まで戻ってくるようなイメージとなる。

for t in np.arange(T,0,-1)-1:        
    if t==T-1 :
        x_smtd[t] = x_fltd[t]
        P_smtd[t] = P_fltd[t]
    else:
        x_smtd[t] = x_fltd[t]+P_fltd[t]*phi/P_fcst[t+1]*(x_smtd[t+1]-x_fcst[t+1])
        P_smtd[t] = P_fltd[t]+(P_fltd[t]*phi/P_fcst[t+1])**2*(P_smtd[t+1]-P_fcst[t+1])

なぜこのようなアルゴリズムになるのか興味がある方は、前回挙げたHamiltonやTsayに加えて、以下の本も参考になると思います。
HamiltonやTsayはファイナンス色が強いので、あまり馴染みのない方はこちらの方がいいかもしれません。

データ同化入門 (予測と発見の科学)

作者: 樋口知之
出版社/メーカー: 朝倉書店
発売日: 2011/09/15
メディア: 単行本
購入: 2人クリック: 8回
この商品を含むブログ (1件) を見る

フィルタリングとスムージングそれぞれについて、平均と上下2標準偏差の信頼区間をプロットしておく。

まずは、フィルタリング。

f:id:lofas:20150225132521p:plain

次に、スムージング。

f:id:lofas:20150225132530p:plain

スムージングの方が、線がより滑らかになっていることが見て取れる。

コード全体はこちら。

# -*- coding: utf-8 -*-
from __future__ import print_function, division
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd

T=200
phi=0.7
sig_eps=1
sig_eta=1

np.random.seed(0)
x=np.zeros((T,1))
y=np.zeros((T,1))

#潜在変数の作成
for t in range(T-1):
    x[t+1]=phi*x[t]+np.random.normal(0,sig_eta,1)
#観測値の作成
for t in range(T):
    y[t]=x[t]+np.random.normal(0,sig_eps,1)
    
#フィルタリングに必要な変数の箱
x_fcst=np.zeros((T,1))
P_fcst=np.zeros((T,1))
y_fcst=np.zeros((T,1))
MSE_y=np.zeros((T,1))
x_fltd=np.zeros((T,1))
P_fltd=np.zeros((T,1))
x_smtd=np.zeros((T,1))
P_smtd=np.zeros((T,1))

#潜在変数の平均と分散の初期値
x_fcst[0] = 0
P_fcst[0] = 1000

#カルマンフィルター
for t in range(T):
    
    y_fcst[t] = x_fcst[t]
    MSE_y[t] = P_fcst[t]+sig_eps**2
    
    x_fltd[t] = x_fcst[t]+P_fcst[t]/MSE_y[t]*(y[t]-y_fcst[t])
    P_fltd[t] = P_fcst[t]-P_fcst[t]**2/MSE_y[t]
    
    if t != T-1 :
        x_fcst[t+1] = phi*x_fltd[t]
        P_fcst[t+1] = phi**2*P_fltd[t]+sig_eta**2

#カルマンスムージング        
for t in np.arange(T,0,-1)-1:        
    if t==T-1 :
        x_smtd[t] = x_fltd[t]
        P_smtd[t] = P_fltd[t]
    else:
        x_smtd[t] = x_fltd[t]+P_fltd[t]*phi/P_fcst[t+1]*(x_smtd[t+1]-x_fcst[t+1])
        P_smtd[t] = P_fltd[t]+(P_fltd[t]*phi/P_fcst[t+1])**2*(P_smtd[t+1]-P_fcst[t+1])

#Kalman Filter
plt.plot(x_fltd-2*np.sqrt(P_fltd),lw=1,color="b")
plt.plot(x_fltd,lw=1,color="k")
plt.plot(x_fltd+2*np.sqrt(P_fltd),lw=1,color="b")
plt.title("Kalman Filter",fontsize=16)

#Kalman Smoothing
plt.plot(x_smtd-2*np.sqrt(P_smtd),lw=1,color="b")
plt.plot(x_smtd,lw=1,color="k")
plt.plot(x_smtd+2*np.sqrt(P_smtd),lw=1,color="b")
plt.title("Kalman Smoothing",fontsize=16)