Feel＆Think

前回は、最も単純な状態空間モデルを用いてカルマンフィルタの方法論を説明しました。今回は、ExcelのVBAを使ってカルマンフィルタの動作を体験してみます。まずは、前回の復習から始めましょう。カルマンフィルタでは、次のような状態空間モデルを前提とします。

観測方程式
「観測値」 ＝ 「状態」 ＋ ノイズ

状態方程式
「状態」 ＝ 「前回の状態」 ＋ ノイズ

例えば、「状態」を「変形している状態」とすれば、「観測値」には「変形量」が対応します。また、カルマンフィルタでは以下の補正の式（フィルタ方程式）を使います。

補正後の状態　＝　補正前の状態 ＋ カルマンゲイン ×（本物の観測値－補正前の状態）

カルマンゲイン ＝ 状態の予測誤差の分散 ÷（状態の予測誤差の分散 ＋ 観測方程式のノイズの分散）

状態の予測誤差の分散 ＝ 前回の状態の予測誤差の分散 ＋ 状態方程式のノイズの分散

補正後の状態の予測誤差の分散 ＝ 補正前の状態の予測誤差の分散 ×（１－カルマンゲイン）

これらによる「状態」の補正の手順をVBAでプログラミングし、カルマンフィルタの効果を見てみます。なお、観測は繰り返して行いますが、観測される値はいつも同じである場合を想定します。興味のある方は、ExcelのVisualBasicのエディタに下記のとおり入力し実行してみてください。
なお、見ているシートのE列の4行目から8行目には、初期値などとして、1，１０，１０，１０，１００を入力しておいてください。

Public Sub kalman()

	y = Cells(4, 5)	'観測値,一定＝1
	xPre = Cells(5, 5) 	'前回の状態,初期値＝10
	pPre = Cells(6, 5)	'前回の状態の予測誤差の分散,初期値＝10
	sigmaW = Cells(7, 5)	'状態方程式のノイズの分散＝10
	sigmaV = Cells(8, 5)	'観測方程式のノイズの分散＝100

	For i = 1 To 20
		'状態の予測（予測値は,前回の値と同じ）
		xForecast = xPre
		'状態の予測誤差の分散
		pForecast = pPre + sigmaW
		'カルマンゲイン
		kGain = pForecast / (pForecast + sigmaV)
		'カルマンゲインを使って補正された状態
		xFiltered = xForecast + kGain * (y - xForecast)
		'補正された状態の予測誤差の分散
		pFiltered = (1 - kGain) * pForecast
		'値の格納
		Cells(10 + i, 6) = xFiltered
		Cells(10 + i, 7) = pFiltered
				'繰り返しのため値を更新
		xPre = xFiltered
		pPre = pFiltered
	Next
End Sub

結果は下図のようになり、観測値をもとにカルマンフィルタで補正を繰り返すと、「状態」の推定値が観測値に収束していきます。収束の速さはノイズの分散で左右されるので、分散の値を変えて試してみてください。

次回は、より一般的な状態空間モデルを用いて説明します。