狠狠撸

北川源四郎「時系列解析入門」第 9 章
柏野雄太
バクフー株式会社
July 8, 2015

注意
この本は名著ではありますが，役割的に「時系列解析のチートシート」
的な存在で，全体的に記述があっさりしすぎています．
特に 9 章，10 章の状態空間モデルの部分は記述が少なすぎますので，
以下の書籍等で足りない部分を補う必要があります．
カルマンフィルタは通常工学部の大学 4 年生から修士 1 年のときに，
半期の講義＋演習でガッツリやるようなレベルです．
状態空間モデル:
樋口知之「予測にいかす統計モデリングの基礎」
樋口知之+「データ同化入門」
カルマンフィルタ:
足立修一?丸田一郎「カルマンフィルタの基礎」
R のパッケージ dlm:
Petris+「R によるベイジアン動的線型モデル」
状態空間モデルの古典的入門書:
コマンダー＆クープマン「状態空間時系列分析入門」
2 / 1

状態空間モデル
xn = Fnxn?1 + Gnvn (9.1)
yn = Hnxn + wn (9.2)
yn: 観測される時系列次元 l
xn: 状態 (state) 次元 k
vn: システムノイズ（状態ノイズ) 平均ベクトル 0, 分散共分散行列 Qn
に従う m 次元の正規白色ノイズ
wn: 観測ノイズ平均ベクトル 0, 分散共分散行列 Rn に従う l 次元の正
規白色ノイズ
Fn: k × k 次元
Gn: k × m 次元
Hn: l × k 次元
3 / 1

状態空間モデルの解釈
xn = Fnxn?1 + Gnvn (9.1)
yn = Hnxn + wn (9.2)
回帰モデルとしての解釈:
観測モデルは観測 yn を表現する回帰モデルで，xn が回帰係数，シス
テムモデルはその時間発展
信号理論としての解釈:
システムモデルは信号発生メカニズム，観測モデルはノイズが付加さ
れる様子を記述．
4 / 1

状態空間モデルの例: AR モデル
AR モデル
yn =
m∑
i=1
aiyn?i + vn (9.3)
xn = Fxn?1 + Gvn (9.4)
xn = (yn, yn?1, . . . , yn?m+1)T
F =
?
?
?
?
?
a1 a2 · · · am
1
...
1 0
?
?
?
?
?
, G =
?
?
?
?
?
1
0
...
0
?
?
?
?
?
(9.5)
H =
(
1 0 · · · 0
)
5 / 1

カルマンフィルタによる状態の推定
観測値 Yj =? xn の推定
Yj が与えられたときの xn の条件付き分布 p(xn|Yj) を求める問題
j < n: 予測 (prediction)
j = n: フィルタ (?ltering)
j > n: 平滑 (smooting)
6 / 1

カルマンフィルタの逐次計算
x1|0 ?→ x2|0 ?→ x3|0 ?→ x4|0 ?→ x5|0 ?→
?
x1|1 =? x2|1 ?→ x3|1 ?→ x4|1 ?→ x5|1 ?→
?
x1|2 ←? x2|2 =? x3|2 ?→ x4|2 ?→ x5|2 ?→
?
x1|3 ←? x2|3 ←? x3|3 =? x4|3 ?→ x5|3 ?→
?
x1|4 ←? x2|4 ←? x3|4 ←? x4|4 =? x5|4 ?→
?
[ノーテーション] ?: フィルタ, =?: 予測, ←?: 平滑化, ?→: 長期予測
8 / 1

9.3 平滑化のアルゴリズム
YN = y1, · · · , yN が与えられたとき，途中状態の xn を推定する．
[固定区間平滑化]
An = Vn|nFT
n+1V ?1
n+1|n
xn|N = xn|n + An(xn+1|N ? xn+1|n)
Vn|N = Vn|n + An(Vn+1|N ? Vn+1|n)AT
n (9.14)
9 / 1

9.4 状態の長期予測
Yn = {y1, · · · , yn} に基づいて，1 期先予想をクリア消し，j 期先の状
態 xn+j(j > 1) を推定する．
カルマンフィルタにより一期先の xn+1|n と xn+1|n を求める．
yn+1 は得られないが，Yn+1 = Yn と仮定する =?
xn+1|n+1 = xn+1|n, Vn+1|n+1 = Vn+1|n
カルマンフィルタの n + 1 期に対する 1 期先アルゴリズムから 2 期先
予測が得られる．
xn+2|n = Fn+2xn+1|n
Vn2|n = Fn+2Vn+1|nFT
n+2 + Gn+2Qn+2GT
n+2 (9.15)
[長期予測]
xn+i|n = Fn+ixn+i?1|n
Vni|n = Fn+iVn+i?1|nFT
n+i + Gn+iQn+iGT
n+i (9.16)
10 / 1

9.5 時系列の予測? BLSALLFOOD の例
Figure: AR モデルの長期予測分布
12 / 1

9.6 時系列モデルの尤度計算とパラメータ推定
時系列モデルの尤度
l(θ) = ?
1
2
{
lNlog2π +
N∑
n=1
log|dn|n?1|
+
N∑
n=1
(yn ? yn|n?1)T
d?1
n|n?1(yn ? yn|n?1)
}
(9.23)
ただし，このままでは計算が大変．
データに欠損がない AR の場合はユール?ウォーカー，最小二乗法,
PARCOR 法などのほうがいい．
仕方なく (9.23) を計算する場合でも，次元圧縮を考える．
13 / 1

9.6 時系列モデルの尤度計算とパラメータ推定
次元圧縮の方法: 時系列の次元が l = 1 で wn の分散が Rn = σ2 で一定の
場合
1 R = 1 としてカルマンフィルタ (カルマンゲインを計算すると R = σ2
でも R = 1 でも同じ)
2 以下の式により ?σ2 を求める
?σ2
=
1
N
N∑
n=1
(yn ? yn|n?1)2
?dn|n?1
(9.28)
3 以下の式により対数尤度 l(θ?) を求める
l(θ?
) = ?
1
2
{
Nlog2π?σ2
+
N∑
n=1
log ?dn|n?1+N
}
(9.29)
4 1, 2, 3 のステップを繰り返して，対数尤度 l(θ?) を最大化して，最尤
推定値 θ? を求める
14 / 1

9.7 欠損値の補間
時系列の状態空間モデルを用いると欠測値があっても厳密に尤度計算
ができ，パラメータの最尤推定値を求めることができる．
l(θ) = ?
1
2
∑
n∈I(N)
{
llog2π + log|dn|n?1+N |
+ (yn ? yn|n?1)T
d?1
n|n?1+N (yn ? yn|n?1)
}
(9.31)
パラメータが決まりモデルが与えられたら，カルマンフィルタの予測
分布 {xn|n?1, Vn|n?1}，フィルタ分布 {xn|n, Vn|n} が求まり，その後で
平滑値 xn|N から欠測値 yn|N = Hnxn|n を求める．
推定誤差分散は dn|N = HnVn|N HT
n + Rn
15 / 1

9.7 欠損値の補間: BLSALLFOOD の例
Figure: 欠測値の補間
16 / 1

Python 版カルマンフィルター
Python 版カルマンフィルターの更新式
from?numpy?import?dot
x?=?dot(F,?x)?+?dot(B,?u)
P?=?dot(F,?P).dot(F.T)?+?Q
y?=?z?-?dot(H,?x)
S?=?dot(H,?P).dot(H.T)?+?R
K?=?dot(P,?H.T).dot(np.linalg.inv(S))
x?=?x?+?dot(K,y)
P?=?(I?-?dot(K,?H)).dot(P)
17 / 1

MATLAB 版カルマンフィルター
MATLAB 版カルマンフィルターの更新関数
function?[xhat_new,P_new,?G]?=?kf(A,?B,?Bu,?C,?Q,?R,?u,?y,?xhat,?P)
??xhat?=?xhat(:);
??u?=?u(:);
??y?=?y(:);
??xhatm?=?A*xhat?+?Bu*u;
??Pm?=?A*P*A’?+?B*Q*B’;
??G?=?Pm*C/(C’*Pm*C?+?R);
??xhat_new?=?xhatm?+?G*(y?-?C’*xhatm);
??P_new?=?(eye(size(A))?-?G*C’)*Pm;
end
18 / 1

狠狠撸

Gunosy2015 07-07

More Related Content

Gunosy2015 07-07