6/23（月）

信号解析第１０回講義録
日時：２０１４年６月２３日
講義内容：カルマンフィルタの使い方
担当者：情報知能工学科小島史男
1
はじめに
先週はカルマンフィルタの構成法について学習しました。データの AR モデルによるあてはめがレビン
ソンのアルゴリズムを用いて効率よく実現できることを示しました。カルマンフィルタの構成を考える
ときの状態空間モデルは時系列信号である AR および ARMA モデルとどう関係するのでしょうか。今
週は AR モデルにカルマンフィルタを適用する方法について考察します。教科書第９章の内容です。
2
ＡＲモデルと状態空間モデル
状態空間モデルの一般系は
x[n] = F [n]x[n − 1] + G[n]v[n]
y[n] = H[n]x[n] + w[n]
で記述されることは先週学習しました。以下では、ＡＲモデルがこの状態空間モデルで記述できること
を示してみましょう。ＡＲモデル
y[n] =
m
∑
ai y[n − i] + v[n]
i=1
において、



F [n] = F = 


x[n] = [y[n], y[n − 1], · · · , y[n − m + 1]]T
a1 a2 · · ·
1 ··· 0
..
. . ..
. .
.
0 ··· 1
am
0
..
.



 ∈ Rm×m





G[n] = G = 

0

1
0
..
.



 ∈ Rm×1


0
H[n] = H = [1, 0, · · · , 0] ∈ R1×m
と設定し、また観測ノイズを考慮しないとすれば、状態空間モデルは
x[n] = F x[n − 1] + Gv[n]
y[n] = Hx[n]
のように書くことができます。システム x[n] の次元 k はここでは AR 次数 m となり、またシステム雑
音の行列 G[n] 観測行列 H[n] はそれぞれ m × 1 および 1 × m の定数行列 (m = 1, l = 1) となります。ま
たシステムノイズ v[n] はスカラーですので分散共分散行列 Q[n] も白色雑音の分散値 σ 2 となります。ま
た AR モデルでは観測雑音 w[n] は存在しません。すなわち R[n] = 0 となります。
N −1
から尤度関数について考えてみ
さて、この状態空間モデルの形で観測データの系列 YN = {˜
y [n]}n=0
ましょう。尤度関数は N 次の同時正規型確率密度関数によって
L(θ) = fN (y1 , y2 , · · · , yN | θ)
によって与えられることは学習しました。ここで θ = {a1 , a2 , · · · , am , σ 2 } です。ところで条件付き期待
値の性質（ベイズの定理）により、
fN (y1 , y2 , · · · , yN | θ) = fN −1 (y1 , y2 , · · · , yN −1 | θ)
×gN (yN | y1 , y2 , · · · , yN −1 , θ)
が成立します。ここで gN は条件付き確率密度関数を意味します。これを繰り返し適用すれば、尤度関
数は
L(θ) = gN (yN | y1 , y2 , · · · , yN −1 , θ)gN −1 (yN −1 | y1 , y2 , · · · , yN −2 , θ)
· · · × g1 (y2 | y1 , θ)
=
N
∏
gn (y[n] | y1 , y2 , · · · , yn−1 , θ)
n=1
のような条件付き期待値のかけ算で与えられることになります。ここで条件付き確率密度関数の条件部
分に観測データ Yn = {˜
y1 , y˜2 , · · · , y˜n } を代入した条件付き確率密度関数 gn (yn | Yn−1 , θ) は観測値 y˜n の
予測分布に相当します。つまり時時刻々の予測分布のかけ算により尤度関数が反復的に計算できること
になります。カルマンフィルタはこの予測計算を簡単に実現できる計算アルゴリズムです。
3
ＡＲモデルを用いたカルマンフィルタ
前節で、AR モデルを状態空間モデルで表すことができましたので、カルマンフィルターを構成するこ
とができます。フィルタの構成は以下のようになります。
予測問題（１ステップ）
：
x
ˆn|n−1 = F x
ˆn−1|n−1
Vn|n−1 = F Vn−1|n−1 F T + σ 2 GGT
濾波問題（フィルタリング）
：
{
K[n] =
1
HVn|n−1 H T
}
Vn|n−1 H T
(
x
ˆn|n = x
ˆn|n−1 + K[n] y˜[n] − H x
ˆn|n−1
)
Vn|n = (I − K[n]H) Vn|n−1
この問題では観測方程式がスカラーであることから、カルマンゲインの更新において逆行列がスカラー
値となるのでの簡単にフィルタを計算することができます。
4
R でカルマンフィルタを実行する
R にはカルマンフィルタのプログラム群がありますが、ここでは教科書にもとづく計算アルゴリズムで
実際に作成してみましょう。つぎのプログラムは教科書１２９ページに記述されているアルゴリズムを
書き下したものです。
✓
カルマンフィルタの計算アルゴリズム
✏
# dim(data) = n; dim(system) = k; dim(noise) = m; dim(observation) = l
kalman <- function(n, k, m, l, kp, x0, xs0, fm, gm, hm, qm, rm, y) {
xh_nm1_nm1 <- numeric(k);
xh_n_nm1 <- numeric(k);
xh_n_n <- numeric(k);
vm_nm1_nm1 <- matrix(0, k, k);
vm_n_nm1 <- matrix(0, k, k);
vm_n_n <- matrix(0, k, k);
xhat <- numeric(n+1);
km <- matrix(0, nrow=k, ncol=l);
xhat[1] <- x0[kp];
xh_nm1_nm1 <- x0;
vm_nm1_nm1 <- xs0 %*% t(xs0) - x0 %*% t(x0);
listt = 1:n;
for(t in listt) {
# One step prediction
xh_n_nm1 <- fm %*% xh_nm1_nm1;
vm_n_nm1 <- fm %*% vm_nm1_nm1 %*% t(fm) + gm %*% qm %*% t(gm);
# Filtering
am <- hm %*% vm_n_nm1 %*% t(hm) + rm;
bm <- diag(l);
hinv <- solve(am,bm);
km <- vm_n_nm1 %*% t(hm) %*% hinv; # Kalman Gain
xh_n_n <- xh_n_nm1 + km %*% (y[t] - hm %*% xh_n_nm1);
vm_n_n <- (diag(k) - km %*% hm) %*% vm_n_nm1;
xhat[t+1] <- xh_n_n[kp];
# Iteration
xh_nm1_nm1 <- xh_n_n;
vm_nm1_nm1 <- vm_n_n;
}
return(xhat);
}
✒
状態空間モデルの次元は (k, l, m) で同じですが、簡単のため k × k 、k × m、l × k の行列は時間 n に依存
しない行列 F, G, H で、またシステム雑音 v[n] および観測雑音 w[n] の分散共分散行列も定数行列 Q, R
で与えています。教科書では記述されていませんが、カルマンフィルタを実行するには、以下のフィル
タの初期値も設定する必要があります。
x
ˆ0|0 = E(x[0])
(
)
V0|0 = E x[0]x[0]T − x
ˆ0|0 x
ˆT0|0
プログラムのなかで % ∗ % は行列の乗算、t(·) は転置行列、また diag(·) は単位行列を表しています。コ
メント#Filtering に続く３行のコマンドは逆行列を求めています。プログラムをみるとこのフィルタを
実行するには下記の情報を事前に与える必要があることがわかります。
(1) 信号のデータ長 n、状態空間表示の次元 (k, m, l)
✑
(2) 状態空間を構成する行列
f m ← F,
gm ← G,
hm ← H
(3) 不規則雑音の分散共分散行列
qm ← Q,
rm ← R
(4) 初期状態の平均値ベクトルと２乗の平均値行列
x0 ← E(x[0]) xs0 ← E(x[0]x[0]T )
(5) 観測信号 {˜
y [i]}N
i=1
なおこのプログラムでは出力情報を簡略化するためカルマンフィルタのすべての情報は取り出さず、状
態ベクトルの kp 番目のコンポーネントのフィルタ情報のみを返しています。
5
R で AR モデルのカルマンフィルタを実現する
カルマンフィルタを実現するには、まず状態空間モデルを構成しなければなりません。ところで AR モデ
ルの場合さて、第８回の講義で説明したレビンソンのアルゴリズムでその構造を決定できます。そこで今
回は、第８回の例にあったアルゴリズム実施例をもとに状態空間モデルを構成することにします。前回の
例では AIC の判定結果から２次の AR モデルとしましたので、状態空間モデルは次数 (k, m.l) = (2, 1, 1)
となります。すなわち各行列は
[
F =
a1 a2
1 0
]
[
G=
1
0
]
H = [1, 0]
となります。また雑音に関する情報は
Q = σ2,
R=0
です。これらに、第８回の講義録で求めた AR 係数および白色雑音の分散の推定値
a
ˆ1 = 1.591856 a
ˆ2 = −0.82323
σ
ˆ 2 = 1.042901
を代入します。またフィルタの初期値設定のため
[
E(x[0]) = [0, 0]
T
T
E(x[0]x[0] ) =
と仮定します。実行画面は次のようになります。
1 0
0 1
]
✓
実行画面
✏
> source("Kalman.R");
> source("signal.dat");
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
n <- 256;
k <- 2;
m <- 1;
l <- 1;
kp <- 1;
x0 <- c(0,0);
xs0 <- diag(2);
a1 <- 1.591856;
a2 <- -0.82323;
sigma2 <- 1.042901;
fm <- matrix(c(a1,a2,1,0), nrow=2,ncol=2);
gm <- matrix(c(1,0), nrow=2, ncol=1);
hm <- matrix(c(1,0), nrow=1, ncol=2);
qm <- matrix(c(sigma2), nrow=1,ncol=1);
rm <- matrix(0, nrow=1,ncol=1);
> xhat <- kalman(n, k, m, l, kp, x0, xs0, fm, gm, hm, qm, rm,y);
>
>
>
>
plot(y, type="l", lty=1);
par(new=T);
plot(xhat, type="l",lty=2);
dev.copy2eps(file="filter.eps", width=10, height=6);
✒
✑
−10
−5
xhat
y
0
5
このカルマンフィルタに第９回の宿題でアップしておいた時系列信号を代入し、フィルタを実行します。
次の図はフィルタリング結果（データ長 n = 256）と元の時系列信号を比較しています。実線が観測デー
タ、点線がフィルタの結果です。
0
50
100
150
150
200
200
250
250
Index
Figure 1: カルマンフィルタの結果
宿題：先週の問題で求めた AR モデルを用いてカルマンフィルタを実行してください。求めたフィルタ
出力の最初の要素 (kp = 1) と観測信号データ yd（先週ダウンロードしたファイルです）とを比較して
みましょう。
6
ARMA モデルをカルマンフィルタに適用するには
ARMA モデルを状態空間表示で記述してみましょう。(l, m) 次の ARMA モデルは
y[n] =
m
∑
l
∑
aj y[n − j] + v[n] −
j=1
bj v[n − j]
j=1
で与えられるとします。v[·] は平均 0 分散 σ 2 の正規性白色雑音とします。ここで
y n+i|n−1 :=
m
∑
aj y[n + i − j] −
l
∑
j=i+1
b[j]v[n + i − j]
j=i
と置きます。ここで、k = max(m, l + 1) とすると ARMA モデルは
y[n] = a1 y[n − 1] + y n|n−2 + v[n]
y n+i|n−1 = ai+1 y[n − 1] + y n+i|n−2 − bi v[n]
y n+k−1|n−1 = ak y[n − 1] − bk−1 v[n]
のような関係式を満足します。そこで、システムの k 次元の状態ベクトルを
(
x[n] = y[n], y n+1|n−1 , · · · , y n+k−1|n−1
)T
と設定すると、AR モデルのときとは少し異なる状態空間表示が可能となります。
x[n] = F x[n − 1] + Gv[n]
y[n] = Hx[n]
上の式の定数行列 F ∈ Rk×k 、G ∈ Rk×1 、H ∈ R1×k はそれぞれ

a1
 .
 ..

F =
 ak−1
ak
1 ···
.. . .
.
.

0
.. 
. 

,
0 ··· 1 
0 ··· 0



G=


1
−b1
..
.
−bk−1



,


H = [1, 0, · · · , 0]
と表示できます。ただし、AR 係数、MA 係数ともそれぞれの次元 (m, l) を超える分については、以下
のように 0 と設定します。
ai = 0 for i > m bi = 0 for i > l
AR モデルと同様に ARMA モデルにおいても、観測雑音 w[n] は存在しません。すなわち分散共分散行
列 R[n] = R = 0 となります。カルマンフィルタを実現するには、システムの構成要素行列 F および G
に含まれる AR 係数および MA 係数を推定する必要がありますが、AR モデルにおけるレビンソンのア
ルゴリズムは使えません。ARMA モデルでは蓄積された観測データから、尤度関数を求め、最大尤度を
とるパラメータを推定値として、係数決定が可能となりますが、これを実行するには AIC の次元決定と
も絡んで AR モデルのときほど簡単には求まりません。詳細は教科書の第１０章および付録 A を参照し
てみてください。次回はいよいよカルマンフィルタを用いた予測の方法について学習していきます。
7
休講のお知らせおよび試験の実施日について
来週６月３０日（月）は休講とします。

Download Report