カルマンフィルタと最尤法
新規作成：2017年04月16日
最終更新：2017年04月16日
カルマンフィルタを実行するには、パラメタを事前に与える必要があります。
そのパラメタを推定する方法が、今回紹介する最尤法です。
この記事では、尤度の説明をしたのちに、ローカルレベルモデルを例とした、状態空間モデルにおける尤度の計算方法を説明します。
またライブラリを使わない自作の尤度計算・最尤法実行メソッドを作って状態空間モデルを推定してみます。
この記事はカルマンフィルタの考え方の続編にあたるものです。
あらかじめこちらの記事を読んでおいた方が理解が深まるかと思います。
ソースコードはまとめてこちらに載せてあります。
スポンサードリンク
目次
- 尤度と最尤法の考え方
- 状態空間モデルにおける尤度計算
- 尤度計算プログラムの実装
- 最尤法の実装
1．尤度と最尤法の考え方
尤度とは、「パラメタを指定したときに、今手持ちのデータを再現することができる確率」を指します。
最尤法とは「尤度が最も大きくなるパラメタを推定する」という考え方を指します。
たとえば、コインを二回だけ投げて結果が「表・裏」だったとします。このコインの「表が出る確率」というパラメタは一体いくつが妥当だと思いますか？
仮説１ パラメタは 1/3 だ！
仮説1を採用した時に「表・裏」という「今手持ちのデータ」が得られる確率は
1/3 × 2/3 ＝ 2/9
仮説２ パラメタは 1/2 だ！
仮説1を採用した時に「表・裏」という「今手持ちのデータ」が得られる確率は
1/2 × 1/2 ＝ 1/4
2/9 ＜ 1/4 より、仮説2を採用した方が「今手持ちのデータを再現することができる確率（＝尤度）」が高い
⇒仮説2「表が出る確率は1/2」を採用する
という流れが最尤法です。
なお、尤度をそのまま使うのではなく、計算の簡単のため対数を取った対数尤度を使うことが多いです。
状態空間モデルにこれを適応します。
まず「テキトーな」パラメタをつかってカルマンフィルターをしてやり「状態」を推定します。
そしてそのときの「今手持ちのデータを再現することができる確率（＝尤度）」を計算してやります。
その「今手持ちのデータを再現することができる確率(＝尤度)」が最も高くなるようにパラメタを変更してやって、最も尤もらしいパラメタを見つけます。
そして、またそのパラメタを使って状態空間モデルを推定しなおしてやれば、最も尤もらしい状態空間モデルが出来上がります。
2．状態空間モデルにおける尤度計算
正規分布と尤度
ローカルレベルモデルでは、確率分布として正規分布を仮定します。
正規分布を仮定した際の尤度の計算方法について、簡単に解説します。
正規分布とは確率分布の一種です。詳しくは「確率密度関数と正規分布」の記事を参照してください。
先ほどのコイン投げの例ですと「コインを投げて表が出る確率」がパラメタとなっていました。
正規分布を使う場合は「期待値（平均値）」と「分散」の2つがパラメタとなります。
尤度の計算に必要な要素
正規分布の確率密度関数を用いて、確率密度を計算する場合は、以下の3つの情報が必要です。
・データ
・期待値
・分散
状態空間モデルでは、以下の値を使います
・データ ＝ 観測値 － 観測値の予測結果 （＝観測値の予測誤差）
・期待値 ＝ ０
・分散 ＝ 観測値の予測誤差の分散
まず、データとして「観測値の予測誤差」を使っていることに注目してください。
予測誤差は上振れすることもあれば下振れすることもあり、その期待値は０だとみなすことができます。
あとは、観測値を予測したときの予測誤差の分散があれば、尤度が計算できますね。
観測値の予測誤差の分散の求め方
去年の「状態」から今年の「観測値」を予測する手順は以下の通りです。
去年 今年
状態 10 →状態方程式→ 10
↓
観測方程式
↓
観測値 10
このとき、2つの方程式が経由されています。
まずは、状態方程式を使って、「去年の状態」から「今年の状態」を予測する際に、予測誤差が生まれます。
また「今年の状態」から「観測値」を予測する際に、また観測方程式を経由するため、観測誤差が付与されます。
よって、「観測値の予測誤差の分散」は以下のようにして計算されます。
今期の状態の予測誤差の分散 ＝ 前期の状態の予測誤差の分散 ＋ 状態方程式のノイズの分散
今期の観測値の予測誤差の分散 ＝ 今期の状態の予測誤差の分散 ＋ 観測方程式のノイズの分散
スポンサードリンク
3．尤度計算プログラムの実装
カルマンフィルタする関数の作成
尤度を計算してみましょう。
カルマンフィルタの考え方で作ったカルマンフィルタ関数を少し改造します。
ソースコードはまとめてこちらに載せてあります。
# --------------------------------------------------------------------------------------
# ライブラリを使わないカルマンフィルタの実装
# 尤度を計算するための情報も出力するようにした
# --------------------------------------------------------------------------------------
localLevelModel_2 <- function(y, xPre, pPre, sigmaW, sigmaV) {
#状態の予測(ローカルレベルモデルなので、予測値は、前期の値と同じ)
xForecast <- xPre
# 状態の予測誤差の分散
pForecast <- pPre + sigmaW
# 観測値の予測誤差
v <- y - xForecast
# 観測値の予測誤差の分散
F <- pForecast + sigmaV
# カルマンゲイン
kGain <- pForecast / (pForecast + sigmaV)
# カルマンゲインを使って補正された状態
xFiltered <- xForecast + kGain * (y - xForecast)
# 補正された状態の誤差の分散
pFiltered <- (1 - kGain) * pForecast
# 結果の格納
result <- data.frame(
xFiltered = xFiltered,
pFiltered = pFiltered,
v = v,
F = F
)
return(result)
}
引数は以下の通りです
y :当期の観測値
xPre :前期の状態
pPre :前期の状態の予測誤差の分散
sigmaW:状態方程式のノイズの分散
sigmaV :観測方程式のノイズの分散
14行目において、「観測値」の予測誤差を計算しています。
17行目において、「観測値の予測誤差の分散」を計算しています。11行目で計算された「状態の予測誤差の分散」に観測方程式のノイズの分散を足しています。
この関数を、ナイル川の流量データに適用します。
詳しいコードの解説はカルマンフィルタの考え方を参照してください。
# ナイル川流量データで計算
Nile
# サンプルサイズ
N <- length(Nile)
# 状態の推定値
x <- numeric(N)
# 状態の予測誤差の分散
P <- numeric(N)
# 「状態」の初期値は0とします
x <- c(0, x)
x
# 「状態の予測誤差の分散」の初期値は10000000にします
P <- c(10000000, P)
P
# 観測値の予測誤差
v <- numeric(N)
# 観測値の予測誤差の分散
F <- numeric(N)
# カルマンフィルタの逐次計算を行う
for(i in 1:N) {
kekka <- localLevelModel_2(Nile[i],x[i],P[i], sigmaW = 1000, sigmaV = 10000)
x[i + 1] <- kekka$xFiltered
P[i + 1] <- kekka$pFiltered
v[i] <- kekka$v
F[i] <- kekka$F
}
観測値の予測誤差を「v」という変数に、「観測値の予測誤差の分散」をFという変数に格納しました。
あとは、これを使って対数尤度を計算するだけです。
なお、vやFには、ナイル川の流量データのサンプルサイズと同じ長さの系列となっていることだけ念のため注意してください。
こんなデータが入っています。
> # 観測値の予測誤差の系列 > v [1] 1120.0000000 41.1187694 -177.4103155 137.7291309 42.8044929 30.0133645 -325.5439019 181.9698000 272.0498772 [10] -32.0492045 -168.3531271 -182.7658464 41.6705418 -85.5944277 -36.4624506 -86.6100528 156.7905194 -266.5696586 [19] -35.5521761 156.0526096 73.8936900 163.9307654 59.6437502 143.5306010 114.7548926 43.7531348 -158.0670485 [28] -45.3642493 -359.1088149 -196.0933351 -109.1175018 -259.6387306 56.5042856 -65.7606982 -179.9950371 83.6317403 [37] -162.9618938 209.0632741 182.5835319 52.2574566 -99.8602199 -177.8823612 -399.8263363 76.1892321 -66.3937622 [46] 369.5429251 249.7086083 -85.7517234 -130.5853626 -38.3069157 -80.9580645 17.9132595 32.0738812 21.4089229 [55] -148.3748306 38.7095516 -72.7480742 -1.0947301 243.2010181 -103.5012477 -53.5397427 44.9243514 12.7877588 [64] 108.3330663 119.0662155 -0.1002622 -75.0731758 133.2083090 -141.7787431 -198.4763350 -171.8567202 71.5714403 [73] 18.2359711 -56.6905897 17.6247252 251.8632962 3.8208622 16.7886325 -13.7469208 31.9668952 -122.6691601 [82] -84.5293245 27.3067976 231.9296966 37.2724484 95.2030649 -119.5166345 38.7715268 80.2971580 -101.3956181 [91] 130.9970380 -18.3926255 -18.4237435 255.5535453 -71.4858325 -218.1734908 13.7674330 -190.9519245 -143.3650759 [100] -78.6341101 > > # 観測値の予測誤差の分散の系列 > F [1] 10011000.00 20990.01 16235.83 14840.78 14261.81 13988.27 13851.15 13780.38 13743.31 13723.73 [11] 13713.35 13707.84 13704.90 13703.34 13702.51 13702.07 13701.83 13701.71 13701.64 13701.60 [21] 13701.58 13701.57 13701.57 13701.57 13701.56 13701.56 13701.56 13701.56 13701.56 13701.56 [31] 13701.56 13701.56 13701.56 13701.56 13701.56 13701.56 13701.56 13701.56 13701.56 13701.56 [41] 13701.56 13701.56 13701.56 13701.56 13701.56 13701.56 13701.56 13701.56 13701.56 13701.56 [51] 13701.56 13701.56 13701.56 13701.56 13701.56 13701.56 13701.56 13701.56 13701.56 13701.56 [61] 13701.56 13701.56 13701.56 13701.56 13701.56 13701.56 13701.56 13701.56 13701.56 13701.56 [71] 13701.56 13701.56 13701.56 13701.56 13701.56 13701.56 13701.56 13701.56 13701.56 13701.56 [81] 13701.56 13701.56 13701.56 13701.56 13701.56 13701.56 13701.56 13701.56 13701.56 13701.56 [91] 13701.56 13701.56 13701.56 13701.56 13701.56 13701.56 13701.56 13701.56 13701.56 13701.56
対数尤度は以下のようにして計算されます。
> # 対数尤度の計算1 > # dnorm関数を使った対数尤度の計算 > sum(log(dnorm(v, mean = 0, sd = sqrt(F)))) [1] -646.3254
まずdnormの解説をします。
これは引数として（データ、期待値、標準偏差）の3つをとり、確率密度を計算してくれる関数です。分散ではなく標準偏差をとるため、F（観測値の予測誤差の分散）の平方根を入れていることに注意してください。
あとは、このdnormの計算結果をすべてのデータでかけ合わせれば、尤度が求まります。
次に、尤度をlog関数を使って対数をとっています。
対数をとると、掛け算は足し算になります。
よって、最後にsum関数を使って、対数をとった確率を足し合わせることにより、対数尤度を求めています。
なお、今回は正規分布を仮定していますので、その確率密度関数を直接使うことで、以下のように計算することもできます。
> # 対数尤度の計算 2 > # 教科書通りの計算式を使った、対数尤度の計算 > -1 * (N/2) * log(2 * pi) - 1/2 * sum(log(F) + v^2 / F) [1] -646.3254
ここで注目してほしいのですが、式の前半部分の『-1 * (N/2) * log(2 * pi)』は、パラメタによらず一定です。
なので、計算量を節約するために、以下のようにします。
dlmパッケージで計算される対数尤度はこの値となります。計算の効率化のための工夫ですね。また、正負が逆になっているので「対数尤度の最大化」は「この指標を最小にする」行為に置き換わります。
> # 対数尤度の計算 3 > # dlmパッケージで計算される対数尤度 > 1/2 * sum(log(F) + v^2 / F) [1] 554.4316
実際にdlmパッケージで計算される対数尤度と比較してみます。
# dlmでの計算結果との比較 # dlmパッケージの読み込み library(dlm) # dlmのパラメタの設定 modelDlm <- dlmModPoly(order=1, m0=0, C0=10000000, dW = 1000, dV = 10000) # カルマンフィルタの実行 Filter <- dlmFilter(Nile, modelDlm)
結果はこちら。同じ値となっていることを確認してください。
> dlmLL(Nile, modelDlm) [1] 554.4316
4．最尤法の実装
では、最尤法してみましょう。
まずは、「パラメタを入れると、対数尤度をすぐに計算してくれる関数」を作ります。
# --------------------------------------------------------------------------------------
# 最尤法の実行
# --------------------------------------------------------------------------------------
calcLogLikLocalLevel <- function(sigma){
# パラメタ（状態・観測方程式のノイズの分散）を入れると、すぐに対数尤度を計算してくれる関数
# データや、「状態」「状態の予測誤差の分散」の初期値の設定
data <- Nile
x0 <- 0
P0 <- 10000000
# 分散は負にならないため、あらかじめEXPをとっておく
sigmaW <- exp(sigma[1])
sigmaV <- exp(sigma[2])
# サンプルサイズ
N <- length(data)
# 状態の予測誤差の分散
P <- numeric(N)
# 「状態の予測誤差の分散」の初期値の設定
P <- c(P0, P)
# 状態の推定値
x <- numeric(N)
# 「状態」の初期値の設定
x <- c(x0, x)
# 観測値の予測誤差の系列
v <- numeric(N)
# 観測値の予測誤差の分散の時系列
F <- numeric(N)
# カルマンフィルタの逐次計算を行う
for(i in 1:N) {
kekka <- localLevelModel_2(Nile[i],x[i],P[i], sigmaW = sigmaW, sigmaV = sigmaV)
x[i + 1] <- kekka$xFiltered
P[i + 1] <- kekka$pFiltered
v[i] <- kekka$v
F[i] <- kekka$F
}
# 対数尤度を返す
return(1/2 * sum(log(F) + v^2 / F))
}
あとで、最適化関数optimを使うために、状態・観測方程式のノイズの分散というパラメタ以外は一切引数に使わないようにしました。
また、分散は負にならないため、あらかじめEXPをとっていることに注意してください。うっかりでsigmaVに-30という値を入れたとしても、エラーにならないようにする配慮です。
正しく動くか確認してみます。
パラメタは、内部でEXPをとっているので、引数として入れる場合はlogを付けて入れておきました。
> # 動作確認 > calcLogLikLocalLevel(sigma=c(log(1000),log(10000))) [1] 554.4316
先ほどの計算結果と一致しているようです。
では、最適化関数optim関数を使って、尤度を最大にするパラメタを推定します。
> # 最尤法の実行
> optim(c(1,1), calcLogLikLocalLevel, method = "L-BFGS-B")
$par
[1] 7.291970 9.622439
$value
[1] 549.6918
$counts
function gradient
29 29
$convergence
[1] 0
$message
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"
optim関数は、汎用最適化関数と呼ばれるもので、「計算に使う関数」と「パラメタの初期値」を与えてやると、「計算結果を最小とするパラメタ」を自動で計算してくれます。なお、今回「calcLogLikLocalLevel」関数で出力している値は、対数尤度と正負が逆になっているので、この値を最小にすれば、対数尤度は最大になります。
最適化のアルゴリズムにはいろいろあるのですが、dlmパッケージの計算結果に合わせるために「method = “L-BFGS-B”」を指定しておきました。
推定されたパラメタは「$par」の欄を見ます。
dlmパッケージを使って計算してみます。
# dlmパッケージを使ったパラメタの推定
buildDlm <- function(theta){
dlmModPoly(order=1, m0=0, C0=10000000, dW=exp(theta[1]), dV=exp(theta[2]))
}
fitDlm <- dlmMLE(Nile, parm=c(1, 1), buildDlm)
結果はこちら
> # dlmパッケージでの計算結果
> fitDlm
$par
[1] 7.291970 9.622439
$value
[1] 549.6918
$counts
function gradient
29 29
$convergence
[1] 0
$message
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"
同じ結果となりました。
なお、内部でEXPをとっているため、実際の分散の値は以下のようになります。
> exp(fitDlm$par) [1] 1468.461 15099.836
状態方程式のノイズの分散が1468.461で、観測方程式が15099.836となりました。
後は、このパラメタを使ってカルマンフィルタを再実行してやればよいです。
これで状態空間モデルが推定できることになります。
参考文献
タイトル・画像はAmazonへのリンクです。
|
カルマンフィルタ ―Rを使った時系列予測と状態空間モデル― 薄い本ですが、内容は濃いです。数式だけでなくR言語を使った解説があるのがポイント。 日本語での解説も丁寧ですのでお勧め。これも良書だと思います。 尤度の計算についても詳しく載っていました。 |
 
|
|
カルマンフィルタの基礎 カルマンフィルタの考え方を理解したければ、まずはこの一冊がおすすめです。 数式は多めですが、数値例も豊富です。検算ができるのが嬉しいですね。 |
 
|
|
状態空間時系列分析入門 このサイトでもよく紹介している状態空間モデルの入門書です。 カルマンフィルタの解説は少ないですが、状態空間モデルというモデルの考え方を解説した、初学者向けの良書です。 |
 
|
スポンサードリンク
