機械学習
時系列解析
カルマンフィルター
データサイエンス

カルマンフィルター(Kalman Filter)ってご存知でしょうか?

多分ほとんどの人がカルマンフィルターってなにそれという感じだと思います。

もしかしたら工学部の人は制御工学で触ってたり、もしくはデータ同化を勉強したことがある人なら知ってるかもしれません。

ぼくは普段は機械学習をやってる人なんですけど、カルマンフィルターというのはデータ処理を学ぶにあたっていい題材だと思っていて、なので今回はそのカルマンフィルターの解説をさせてもらおうと思います。


カルマンフィルターについて

まず始めに、カルマンフィルターについての簡単な説明をさせてもらおうと思います。

カルマンフィルターとは、Wikipediaによると


カルマンフィルター (Kalman filter) は、誤差のある観測値を用いて、ある動的システムの状態を推定あるいは制御するための、無限インパルス応答フィルターの一種である。


とあります。

これじゃあんまりイメージつかないと思うので、もう少し具体的に説明してみようと思います。

例えば、なんかしらの現象に対して数理モデルがあるとします。ここでは試験管の中でのミジンコの個体数とか天気予報とか。

これらの各モデルに対して各時刻での値を測るんですけど、実際に観測できるのは正確な値ではなくて、観測された値には必ず誤差が入ってますよね。ミジンコの個体数なんて正確に測れるわけありません。だってめっちゃ小さいし。

カルマンフィルターの基本アイデアは、こうした誤差が入った情報しか手に入らない状況下で、なるべく次の時刻での情報を正確に予想しようというのにあります。

ちなみにこのカルマンフィルターはロケットの軌道推定に有用だということでアポロ計画で利用されたりしています。


機械学習との関わり

機械学習とは既知のデータ群から特徴量を見抜いて予測を立てるアルゴリズムです。

カルマンフィルターとは上述の通り時系列データを次々に予測するアルゴリズムなのですが、このように逐次予測していくアルゴリズムを機械学習ではオンライン機械学習と言ったりします。

まあつまりカルマンフィルターは機械学習の一種であるわけですね。

ただ、よく知られている一般的な機械学習(ニューラルネットとか)とカルマンフィルターが微妙に違うのは、データ予測において内部の数理モデルを仮定するという点です。

このような内部モデルを仮定するのをモデルベース(Model based)と言ったりするらしいのですが、これが正式名称かは知りません(笑)

しかしここでは便宜上モデルベースという呼称が便利なので、ここではカルマンフィルターのような内部モデルを仮定する機械学習アルゴリズムをモデルベース型機械学習と呼ぶことにします。

さて、このように機械学習にはいろんな枠組みがあるので説明が難しいところなのですが、主に内部モデルを仮定するかどうかだったり静的か動的かどうかだったりといろんな区分けがある中で、カルマンフィルターはモデルベース型オンライン機械学習ということとなります。


問題設定

先ほどのカルマンフィルターについての説明からカルマンフィルターは離散時系列データだったら基本的になんでも良いのですが、せっかくですので先ほど述べた試験管におけるミジンコの個体数モデルを問題設定としましょう。

ミジンコの個体数増加モデルは以下のような記述されるとします。

\frac{dx}{dt} = rx(1-\frac{x}{K})

ここで$x$は個体数を表し、$r$が増加率、$K$が環境収容力を表します。

これを実際にグラフとして表すと以下のようになります。

sample.png

最初はすごい勢いで増えるけどだんだんと環境のミジンコの限界収容数に達してきてだんだんと厳しくなってくる感じがちゃんと可視化されているかと思います。

コードは以下の通りとなっています。

import numpy as np

import matplotlib.pyplot as plt

dt = 0.0025
init = 0.0001
max_steps = 100000

def function(x, r = 0.1, K = 20.0):

return x + dt * (r * x * (1.0 - x / K))

def integrate(max_steps=max_steps):
data_steps = np.zeros([max_steps+1])
data_steps[0] = init
x = init

for i in range(max_steps):
x = function(x)
data_steps[i+1] = x

return data_steps

ちなみにこの常微分方程式はベルヌーイ型微分方程式で、初等的に解くことができます。答えはロジスティック曲線と同じ形をしたものとなるので、興味がある人は自分で解いてみてください。

さて、今回は$\Delta t$ごとに状態を観測するとし、観測誤差として各値に平均0、分散0.5を持った正規分布に従うノイズを足します。

sample_noise.png

ちょっとステップ数が細かいのですごいことになってますが、実際にノイズが足し込めてるのがわかると思います。

さて、今回はこのような試験管が10本用意してあるものとします。

問題設定についてはこれでだいたい完了したので、実際にカルマンフィルターを使ってゴリゴリ予測していきたいと思います。


カルマンフィルターの実装

本節ではカルマンフィルターの数式的な解説を行います。

以下で用いる記号として、

$x_t$:時刻tでの個体数

$y_t$:時刻tでの観測値

$M_t$:時刻tでの線形化モデル

$P_t$:時刻tにおける(予測と観測データの)誤差分散行列

$H$:観測演算子

とします。

線形化モデル$F_t$について、先ほど述べた通り、ミジンコの個体数モデルは常微分方程式を用いて表すことができましたが、今回は離散化しており、時系列データのモデルに基づく遷移は線形写像(=行列の掛け算)で表すことになります。

観測演算子$H$について、場合によって観測によって欲しい情報が全て得られるとは限りません。

例えば車の位置座標が欲しいけど実際に得られるのは距離だけで、この場合はユークリッド距離として$L^2$ノルムで与えられ、位置座標をベクトルのような形で得ることができません。このように、予測値を観測情報と同じ形にするために観測演算子というものが用いられます。

さて、準備ができたのでようやくカルマンフィルターに入るわけですが、主にカルマンフィルターは予測ステップ更新ステップに分かれます。

Screen Shot 2019-01-21 at 2.36.08.png

予測ステップは、

x_{t+1} = M_t x_t' \ \ (+ B_t u_t)\\

P_{t+1} = M_{t+1} P_t' M_{t+1}^T

となり、更新ステップは

K_t = P_t H^T (HP_tH^T + R)^{-1} \\

x_t' = x_t + K_t (y_t - Hx_t) \\
P_t' = (I - K_t H) P_t

となります。

ここで更新ステップにおいて求められる更新値と誤差行列をプライム'をつけることで予測ステップと更新ステップでは計算しているものが違うことを強調しました。

さて、カルマンフィルターのイメージとして、予測したあとに観測誤差にある程度近い値に更新する、みたいな感じになります。

Screen Shot 2019-01-21 at 3.10.33.png

ここでは簡単のために予測ステップにおいて線形化モデルによる遷移のみで予測を行いましたが、制御モデルではここで制御信号として$u_t$が入り込むのですが、ここでは簡単のため制御信号を不要とするようなモデルを用いました。


導出

本節では上述したカルマンフィルターの導出を見ていきます。

最初に述べた通り、カルマンフィルターの基本コンセプトはノイズが入った情報を用いてなるべく正確に次の情報を予測するというものです。

さて、カルマンフィルターで目標とするべきは誤差にオーバーフィッティングせずに予測値をなるべく真の値に近づけることです。

つまりカルマンフィルターでは誤差の分散を最小化するように最適化する(最小分散推定)のが目標となるわけです。


予測式

カルマンフィルターの更新式と予測式を導出するにあたって、時刻tのモデルの真の値(true value)を$\hat{x}_t$と表すことにしましょう。

まず、予測値と真の値の差(=ノイズ)$epsilon$は以下の通りとして表せます。

\epsilon = x_t - \hat{x}_t

ここで予測値は真の値にノイズを足し込んだものであると解釈することで、ノイズを以下のように表すことができます。

\epsilon = M(\hat{x}_{t-1} + \epsilon') - M(\hat{x}_{t-1})

ここで$M(\cdot)$をモデルを作用させたものとします。

このノイズの計算式において、第一項の$M(\hat{x}_{t-1} + \epsilon')$をテイラー展開することによって以下の式を得ることができます。

\epsilon = M(\hat{x}_{t-1} + \epsilon') - M(\hat{x}_{t-1}) \\

= \frac{\partial M}{\partial x}(\hat{x}_{t-1})\epsilon' + O(\epsilon'^2)

ここでモデル$M(\cdot)$の接線形モデル(ヤコビアン)を以下の通り定義しておきます。

M_t := \frac{\partial M}{\partial x}(\hat{x}_{t})

また、誤差$\epsilon$をランダウの記法$O(\cdot)$を用いて表しましたが、これを近似値として以下の通りとなります。

\epsilon \approx \frac{\partial M}{\partial x}(\hat{x}_{t-1})\epsilon' \\

= M_{t-1} \epsilon'

さて、これらを元に予測誤差分散行列$P_t$の導出式は以下のようにして得られます。

P_t = \langle \epsilon \ \epsilon^T \rangle \approx \langle M_{t-1}\epsilon' \ \epsilon'^T M_{t-1}^T \rangle = M_{t-1} \langle \epsilon' \ \epsilon'^T \rangle M_{t-1}^T 

これより以下を得ます。

P_t = M_{t-1} P_{t-1}' M_{t-1}^T


更新式

カルマンフィルターにおいて肝となるのが予測の後の更新式で、この更新ステップを経ることで観測結果をそのまま使う場合よりも良い予測を行うことができるというわけです。

さて、実際の更新式を以下のように定義します。

x_t' = x_t + K(y_t - H(x_t))

ここで観測誤差$y_t - H(\hat{x}_t)$を$\epsilon^o$として表すことにすると上式は以下のように表せます。

x_t' - \hat{x}_t = (x_t - \hat{x}_t) + K\{[y_t - H(\hat{x}_t)] - [H(x_t) - H(\hat{x}_t)]\} \\

\Leftrightarrow \epsilon' = \epsilon + K\{\epsilon^o - [H(x_t) - H(\hat{x}_t)]\}

ここで観測演算子について以下のように線形近似します。

H(x_t) = H(\hat{x}_t + \epsilon) \approx H(\hat{x}_t) + H\epsilon

これより以下の式を得られます。

\epsilon' = \epsilon + K(\epsilon^o - H\epsilon) \\

= (I - KH)\epsilon + K\epsilon^o

よって更新ステップにおける予測誤差分散行列を以下のように記述することができます。

P_t' = \langle \epsilon' \epsilon'^T \rangle \\

\approx \langle [(I - KH)\epsilon + K\epsilon^o)((I - KH)\epsilon + K\epsilon^o]^T \rangle \\
= (I - KH) \langle \epsilon \epsilon^T \rangle (I - KH)^T + K \langle \epsilon^o (\epsilon^o)^T \rangle K^T \\

ここで2行目から3行目においてにおいて$\epsilon$と$\epsilon^o$には相関がないと仮定しました。

そうして以下を得ます。

P_t' = (I - KH) P_t (I - KH)^T + K R K^T

さて、カルマンフィルターは誤差の分散を最小化するように最適化する(最小分散推定)のが肝でしたが、これはつまり誤差分散行列$P_t'$のトレースを最小化するようなKを探すことと一致します。

よって

\frac{\partial \ tr(P_t')}{\partial K} = 0

を満たす$K$を求めることとなります。

すでに$P_t'$を$K$に関する式で表すことができているので、行列の微分に関する公式

\frac{\partial \ tr(XYX^T)}{\partial X} = X(Y + Y^T) \\

\frac{\partial \ tr(XY)}{\partial X} = Y^T

を使いながらひたすらゴリゴリ計算していくと、最終的に$K$を得ることができます。

K = P_t H^T (H P_t H^T + R)^{-1}

また、これを$P_t'$に関する式に代入することで以下を得られます。

P_t' = (I - KH) P_t


実装

長いこと数式をいじりましたので、今度はプログラミングの出番です。

改めてカルマンフィルターの数式をまとめると、予測ステップは、

x_{t+1} = M_t x_t' \ \ (+ B_t u_t)\\

P_{t+1} = M_{t+1} P_t' M_{t+1}^T

となり、更新ステップは

K_t = P_t H^T (HP_tH^T + R)^{-1} \\

x_t' = x_t + K_t (y_t - Hx_t) \\
P_t' = (I - K_t H) P_t

なります。

ではこれを素直に実装していきましょう。

ちなみにここではPython + Numpyでいきます。

では、まず最初にカルマンフィルターを作る前にモデルの接線形モデルについて考えなければいけませんが、接線形モデルは一階微分の形をしていて、今回は一回常微分方程式の形をしたモデルを仮定してるのでヤコビアンは数理モデルの式に値を対角成分に持つ行列を作るだけで完了します。

def model(x, r = 0.1, K = 20.0)):

M = np.zeros([len(x)][len(x)])
for i in range(len(x)):
M[i][i] = r * x * (1 - x[i] / K)

return M

ということでまず予測ステップ。

def forecast(x, P):

M = model(x)
x = M.dot(x)
P = M.dot(P).dot(M.T)

return x, P

次に更新ステップ。

def update(x, y, P):

S = np.linalg.inv(H.dot(P).dot(H^T) + R)
K = P.dot(H).dot(S)
x = x + K.dot(y - H.dot(x))
P = (I - K.dot(H)).dot(P)

return x, P

ここで実装に際して、Iは単位行列で、観測演算子Hも今回は全点観測可能なので単位行列です。

また、観測誤差の分散は0.5であると仮定しているのでRは対角成分が0.5で他の成分は0である(他の観測に影響を及ぼさないと仮定)ので0.5*Iとして表すことができます。


結果

さて、実装した結果は以下の通りです。

result.png

相変わらずステップの刻みが細かすぎてアレなのですが、一応足し込んだノイズよりもある程度小さい幅でモデルの予測ができてるような気がします。


まとめ

今回はカルマンフィルターを紹介しました。

機械学習をやるにあたってこういった最適化アルゴリズムを自分でスクラッチで書けると機械学習erとしてワンランク上がるんじゃねって思ったのですが、でもカルマンフィルターってあんまり有名じゃないなぁって思ったので今回は実装と理論を共に載せて紹介することとしました。

今回はちょっと力尽きてデータの可視化などにかなり手を抜いてしまったのですが、今回これを読んで興味を持った人はぜひとも色々試してみてください。

このカルマンフィルターは結構他の手法と合わせることができて、例えばアンサンブル学習の手法を合わせてみたり、もしくはベイズ最適化を合わせてみたり、まあ工夫の方法は色々あります。

また、これらの手法を少し工夫すればデータ同化の最新の研究に行き着きます。

ぜひ皆さんもカルマンフィルターに挑戦してみてください。

ではでは、お疲れ様でした!