カルマンフィルタ(KF) - WHAT' CHA GONNA DO FOR ME?

ある時系列データが観測値として得られているとして、これらの時系列データがどのようにして生成されているか、その仕組みに興味があるとする。
ただ、いま手元にあるのは観測値として得られた時系列データだけなので、その背後にある仕組みについては何かしらの方法で推定する必要がある。
状態空間モデルというフレームワークを用いることで、この推定が可能となる。

今回は、次のようなシンプルな状態空間モデルを考えてみたい。

${ \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} }$

ここで、 ${y_t}$ は ${t}$ 期の観測値、 ${x_t}$ は ${t}$ 期の観測値を駆動する潜在変数とする。
${x_t}$ がAR(1)モデルに従いながら背後で動いていて、この挙動が興味の対象であるものの、我々が観測できるのは ${x_t}$ にノイズが加わった ${y_t}$ だけ、ということをこのモデルは言っている。

上記の式のうち、上の式は観測方程式、下の式は状態方程式と呼ばれており、これらのセットが状態空間モデルである。
観測方程式、状態方程式が共に線形で、かつノイズが正規分布に従う場合に限り、カルマンフィルタというアルゴリズムにより潜在変数 ${x_t}$ の推定は効率的に行われることが知られている。

というわけで、前置きが長くなったが、今回はシミュレーションデータを使ってカルマンフィルタをPythonで書いてみたよ、というお話。

まず、観測値 ${y_t}$ 、潜在変数 ${x_t}$ のそれぞれを以下のように作る。

In [19]: 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)
    ...: 
    ...: 
    ...: #観測値と潜在変数のプロット
    ...: plt.plot(np.hstack((y,x)))
    ...: #plt.title("Observed Value vs Latent Variable")
    ...: plt.legend(['Observed Value','Latent Variable'])

Out[19]: <matplotlib.legend.Legend at 0x11d0f8f0>

f:id:lofas:20150221142701p:plain

潜在変数に対してノイズが加わっている観測値は、その分変動が大きいことが見て取れる。

なお、モデルの未知パラメータである ${\phi}$ 、 ${\sigma_{\epsilon}}$ 、 ${\sigma_{\eta}}$ は、それぞれ0.7、1、1をセットした。

次はメインのカルマンフィルタ部分。

In [29]: 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

なぜこのようなアルゴリズムになるのか興味がある方は、例えば以下の本を参照してください。

Time Series Analysis

作者: James D. Hamilton
出版社/メーカー: Princeton Univ Pr
発売日: 1994/01/11
メディア: ハードカバー
購入: 1人クリック: 5回
この商品を含むブログ (8件) を見る

Analysis of Financial Time Series (Wiley Series in Probability and Statistics)

作者: Ruey S. Tsay
出版社/メーカー: Wiley
発売日: 2010/08/30
メディア: ハードカバー
クリック: 5回
この商品を含むブログを見る

推定された潜在変数と真の潜在変数とを比べてみる。

f:id:lofas:20150221142723p:plain

まあ大体トラックできているような感じはする。
この差分も念のためプロットしておく。

f:id:lofas:20150221142739p:plain

最後に、観測値とも一緒にプロットしてみる。

f:id:lofas:20150221142755p:plain

観測値に含まれるノイズに騙されることなく、よりマイルドな動きをする真の潜在変数に近い値が推定できていることが、なんとなく見て取れる。

今回はデータをシミュレーションから作ったため、パラメータはあらかじめこちらで与えた。
しかし、実際の応用では観測データしか手元に無いため、潜在変数とあわせてパラメータも推定の対象となることが多い。
これについては、また別の機会に書くかもしれないし書かないかもしれない。

コード全体はこちら。

# -*- coding: utf-8 -*-
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)
    
#観測値と潜在変数のプロット
plt.plot(np.hstack((y,x)))
#plt.title("Observed Value vs Latent Variable")
plt.legend(['Observed Value','Latent Variable'])

#フィルタリングに必要な変数の箱
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_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

#真の潜在変数と推定された潜在変数のプロット
plt.plot(np.hstack((x,x_fltd)))
plt.legend(['True Latent','Estimated Latent'])

#推定された潜在変数の真の潜在変数との差分
plt.plot(x-x_fltd)

#真の潜在変数と推定された潜在変数とプロット
plt.plot(np.hstack((x,x_fltd,y)))
plt.legend(['True Latent','Estimated Latent','Observed Value'])