統計的システム論スペクトル解析と状態空間モデル

統計的システム論
スペクトル解析と状態空間モデル
林　和則
京都大学大学院情報学研究科
E-mail: [email protected]
講義資料：　http://www.msys.sys.i.kyoto-u.ac.jp/~kazunori/sst.html
1
2
AGENDA
v スペクトル密度
v Wiener-Khinchineの公式
v 線形時不変システム
v 線形確率システム
v 線形定常確率システム
v 一般状態空間モデル
v 参考文献
3
フーリエ解析
連続時間
(
周フーリエ級数：
期
1
離
cn =
T
散
離散時間
離散フーリエ変換：
T
2
x(t)e
j 2 Tn t
N 1
dt
T
2
1
x(n) =
N
cn e
)
n=
（
非フーリエ変換：
周
期
連 X(⇥) =
続
1
x(t) =
2
x(n)e
j 2N nk
n=0
j 2 Tn t
x(t) =
Xk =
N 1
2
Xk ej N nk
k=0
離散時間フーリエ変換：
x(t)e
j t
dt
X(⇥) =
x(n)e
j⇥n
n=
X(⇥)ej t d⇥
1
x(n) =
2
X(⇥)ej⇥n d⇥
）
4
Wiener-Khinchineの公式
-スペクトル密度：
{x(t), t = 0, ±1, · · · } 平均0の２次定常過程とし，
Rxx (l) が l=
|Rxx (l)| <
相関関数　とする
Sxx ( ) =
Rxx (l)e
j l
,
⇥<
l=
Sxx (z) =
Rxx (l)z
l
l=
-Winer-Khinchineの公式：
1
Rxx (l) =
2⇥
⇥
⇥
Sxx ( )ej l d
<⇥
5
スペクトル密度の性質
(ⅰ) Sxx ( ) = Sxx ( ) 対称性
(ⅱ)　非負性
Sxx ( ) 0
(ⅰ) Rxx (l) = Rxx ( l) より明らか
(ⅱ)　lim
N
1
E
2N + 1
1
= lim
N
2N + 1
x(l)e
N
N
= Sxx ( )
N
E{x(l)x(k)}e
j (l k)
l= N k= N
1
l= 2N
j l
l= N
2N
= lim
2
N
|l|
2N + 1
Rxx (l)e
j l
(∵チェザロ和に関する補題(3))
6
線形時不変システム
u(t)
y(t)
G(z)
-インパルス応答： {g(i), i = 0, 1, · · · }
g(i) = 0, i < 0
(因果的　）
-伝達関数： G(z) =
g(i)z
i=0
i
G( ) =
g(i)e j i
(c.f. 周波数応答：　）
i=0
y(t)
u(t)
-安定：任意の有界入力　に対して出力　が有界
i=0
|g(i)| < ⇥
すべての極が単位円内部に存在
7
線形時不変システムと確率過程
Suu (z)
u(t)
G(z)
は安定とし，入力　は平均値0，スペクトル密度　y(t)
の２次定常過程とする．このとき出力　は２次定常過程で
あり，そのスペクトル密度は
Syy (z) = G(z)G( z)Suu (z)
(Syy ( ) = |G( )|2 Suu ( ))
y(t)
の分散は
1
2
⇤y =
2⇥
⇥
⇥
|G(ej )|2 Suu ( )d
8
線形時不変システムと確率過程(続)
u(t) が印加されて十分時間が経過した後の出力：
-入力　y(t) =
g(i)u(t
i)
i=0
y(t)
-　の平均：
y¯ = E{y(t)} = 0
Ryy (l) = E{y(t + l)y(t)}
y(t)
-　の相関関数：
=
g(i)g(k)E{u(t + l
i)u(t
g(i)g(k)Ruu (l + k
i)
i=0 k=0
=
i=0 k=0
k)}
9
線形時不変システムと確率過程(続々)
2
Ryy (0)
i=0 k=0
|g(i)||g(k)|
2
u
2
u
i=0
|g(i)|
<⇥
y(t) は２次過程
y(t)
-　のスペクトル密度：
Syy (z) =
Ryy (l)z
l
l=
=
z
g(i)g(k)Ruu (l + k
l
i=0 k=0
l=
=
i)
g(i)z
i
i=0
= G(z)G(z
g(k)z k
k=0
1
)Suu (z)
z
l=
(l+k i)
Ruu (l + k
i)
10
ベクトル定常過程のスペクトル解析
-平均0のn次元ベクトル定常過程： {x(t), t = 0, ±1, · · · }
x(t)
Rxx (l) = E{x(t + l)xT (t)}
-　の共分散行列：
x(t)
-　のスペクトル密度行列：
Sxx (z) =
Rxx (l)z
l
l=
Sxx ( ) =
Rxx (l)e
j l
,
⇥<
l=
-Wiener-Khinchineの公式：
1
Rxx (l) =
2⇥
⇥
⇥
Sxx ( )ej l d
<⇥
11
スペクトル密度行列の性質
(ⅰ) Sxx ( ) = STxx ( )
エルミート性
(ⅱ)　非負定値性
Sxx ( ) 0
(ⅰ) Sxx ( ) =
Rxx (l)e
j l
l=
=
RTxx ( l)e
j l
l=
=
RTxx (l)ej
l
= STxx (
)
l=
x1 (t)
..
x(t) =
,
.
(ⅱ)　xn (t)
n
n
=
ai xi (t)
とすると
i=1
n
R (l) =
n
n
ai ak E{xi (t + l)xk (t)} =
i=1 k=1
ai ak Rik (l)
i=1 k=1
12
スペクトル密度行列の性質(続)
両辺をフーリエ変換して
n
n
S⇥⇥ ( ) =
ai ak Rik (l) e
l=
n
j l
i=1 k=1
n
=
ai ak Sik ( )
0
Sik ( ) : Sxx ( の(i,k)成分
)
i=1 k=1
Sxx ( )
任意の a1 , · · · , an について成立するので，　は
非負定値
13
例：１次の自己回帰(AR)モデル
白色ガウス雑音
wt
N (0,
2
)
1
z
G(z) =
入出力関係： yt+1 = ayt +
1
a2
a
yt
a2 wt
2
入力のスペクトル密度： Sww ( ) = ⇥
出力の相関関数： Ryy (l) =
=
g(i)g(k)Rww (l + k
i)
i=0 k=0
2 |l|
a
⇥ 2 (1 a2 )
出力のスペクトル密度： Syy ( ) = |G(e )| Sww ( ) =
1 + a2 2a cos
j
2
14
線形離散確率システム
wt
Gt
xt+1
z
1
xt
vt
Ht
yt
Ft
-状態空間モデル： xt+1 = Ft xt + Gt wt
yt = Ht xt + vt , t = 0, 1, · · ·
Gt
¯0
x
x0
：状態ベクトル (　は平均　，共分散行列　のガウス確率変数)
0
Rp
：観測ベクトル
Rr ：プラント雑音 (平均0のガウス白色雑音ベクトル)
Rp
：観測雑音 (平均0のガウス白色雑音ベクトル)
Rn n ：状態遷移行列
Qt St
wt
T
T
w
v
E
=
s
s
vt
STt Rt ts
Rn r ：駆動行列
Ht
Rp
xt
yt
wt
vt
Ft
Rn
n ：観測行列
x0 は　と独立
wt , vt , t = 0, 1, · · ·
15
例：太陽のウォルフ黒点数のモデル
-YuleによるARモデル( yt :黒点数)：
yt = 1.34254yt 1 0.65504yt
2
-状態空間モデルによる表現１：
0
0.65504
xt+1 =
xt +
1 1.34254
yt = 0
+ wt + 13.854
0.65504
wt
1.34254
1 xt + wt
-状態空間モデルによる表現２：
t+1
=
1.34254
1
yt = 1
0
t
0.65504
0
t+
1
w
0 t+1
16
状態ベクトルの性質
xt , t = 0, 1, · · · はガウス-マルコフ過程
xt = Ft
1 xt 1
+ Gt
= Ft
1 (Ft 2 xt 2
= Ft
..
.
1 Ft 2 xt 2
1 wt 1
+ Gt
+ Ft
2 wt 2 )
+ Gt
1 Gt 2 wt 2
1 wt 1
+ Gt
1 wt 1
t 1
=
(t, s)xs +
(t, k + 1)Gk wk
k=s
推移行列
t 1
=
(t, 0)x0 +
(t, s) =
Ft 1 Ft 2 · · · Fs , t > s
I, t = s
(t, k + 1)Gk wk
k=0
-ガウス確率変数 x0 , w0 , · · · , wt 1 の線形結合はガウス確率変数
xk と　ws , · · · , wt 1 は独立
-任意の k s < t について　p(xt |xs , xs 1 , · · · , x0 ) = p(xt |xs )
17
状態ベクトルの性質(続)
ガウス過程は平均と共分散行列で完全に特徴付けられる
¯ t = E{xt }, Cxx (t, s) = E{[xt
x
¯ t ][xs
x
¯ s ]T }
x
xt の平均および共分散行列はそれぞれ次で与えられる
¯t =
x
Cxx (t, s) =
(t, 0)¯
x0
⇥(t, s) s , t s
T
t ⇥ (s, t), t < s
ただし，
t
E
= Cxx (t, t)
wt
vt
wsT
vsT
=
Qt
STt
St
Rt
ts
t 1
= ⇥(t, 0)
T
⇥
(t, 0) +
0
⇥(t, k + 1)Gk Qk GTk ⇥T (t, k + 1)
k=0
18
状態ベクトルの性質(続々)
平均：
t 1
E{xt } = E
(t, 0)x0 +
(t, k + 1)Gk wk
=
k=0
˜ 0 = x0
x
s
共分散行列（ t のとき）：
¯0
x
t 1
Cxx (t, s) = E
(t, 0)¯
x0
T
s 1
⇥(t, 0)˜
x0 +
⇥(t, k + 1)Gk wk
⇥(s, 0)˜
x0 +
k=0
⇥(s, l + 1)Gl wl
l=0
s 1
= ⇥(t, 0)
T
0 ⇥ (s, 0)
+
⇥(t, k + 1)Gk Qk GTk ⇥T (s, k + 1)
k=0
s = t とすると　の式が得られる．
t
t
t+1
= Ft
T
t Ft
+ Gt Qt GTt
s のとき　より
(t, 0) = (t, s) (s, 0)
s 1
Cxx (t, s) = ⇥(t, s) ⇥(s, 0)
0⇥
T
(s, 0) +
⇥(s, k + 1)Gk Qk GTk ⇥T (s, k + 1)
k=0
= ⇥(t, s)
s
19
観測ベクトルの性質
E
wt
vt
wsT
vsT
=
Qt
STt
St
Rt
ts
yt
はガウス過程であり，平均と共分散行列は次で与えられる
¯ t = Ht x
¯ t = Ht (t, 0)¯
y
x0
Cyy (t, s) =
Ht ⇥(t, s) s HTs + Ht ⇥(t, s + 1)Gs Ss , t > s
Ht t HTt + Rt , t = s
Ht t ⇥T (s, t)HTs + STt GTt ⇥T (s, t + 1)HTs , t < s
共分散行列（ t >　s のとき）：
yt
¯ t = Ht (xt
y
¯ t ) + vt
x
t 1
= Ht (t, s)˜
xs + Ht
(t, k + 1)Gk wk + vt
˜ s = xs
x
¯s
x
k=s
Cyy (t, s) = E{[yt
=E
¯ t ][ys
y
¯ s ]T }
y
t 1
T
˜ s + vs ]
(t, k + 1)Gk wk + vt [Hs x
Ht (t, s)˜
xs + Ht
k=s
20
線形定常確率システム
xt+1 = Ft xt + Gt wt
wt
vt
E
yt = Ht xt + vt
wsT
vsT
=
Qt
STt
St
Rt
ts
Ft , Gt , Ht , Qt , Rt , St が　に依存しない）
t
時間不変（　xt+1 = Fxt + Gwt
xt0
yt = Hxt + vt , t = t0 , t0 + 1, · · ·
推移行列：
(t, s) =
¯t =
平均： x
Ft 1 Ft 2 · · · Fs , t > s
I, t = s
(t, t0 )¯
xt0 = Ft
t0
N (¯
xt0 ,
t0 )
(t, s) = Ft
s
, t
¯ t0
x
t t0 1
共分散行列：
t
= Ft
t0
T t t0
t0 (F )
+
Fk GQGT (FT )k
k=0
t+1
=F
T
tF
+ GQGT
s
21
状態ベクトルの性質
(| i (F)| < 1, i = 1, . . . , n) とする．　と
t0
F は漸近安定　xt
すると，　は平均値0，共分散行列
Cxx (t
s) =
Ft s , t s
(FT )s t , t < s
のガウス-マルコフ過程になる．ただし，　は次の一意的な解
= F FT + GQGT
F は漸近安定なので
t0
lim
t
=
t0
lim Ft
t0
Fk GQGT (FT )k
= 0 よって
t0
¯t = 0
lim x
tに依存しない，解になっている
k=0
一意性の確認：　も解であるとすると
= F(
)FT = Fk (
)(FT )k
0 (k
)
22
倍化(doubling)アルゴリズム
T
T
=
F
F
+
GQG
=
0
初期値　に対する
の解を
t+1
t
0
とする．このとき，アルゴリズム
t
Nk := Mk 1 Nk 1 MTk 1 + Nk 1
Mk := M2k 1
を用いると Nk = 2k ．ただし， M0 = F, N0 = GQGT
1
= GQGT = N0
2
= FN0 FT + GQGT = M0 N0 MT0 + N0 = N1
3
= FN1 FT + GQGT
4
= F2 N1 (FT )2 + FGQGT FT + GQGT
= M1 N1 MT1 + N1 = N2
23
観測ベクトルの性質
(| i (F)| < 1, = 1, . . . , n) とする．　と
t0
F は漸近安定　yt
すると，　は平均値0，共分散行列
Cyy (t
s) =
HFt s HT + HFt s 1 GS, t > s
H HT + R, t = s
H (FT )s t HT + ST GT (FT )s t 1 HT , t < s
の定常ガウス過程になる．
24
同時密度関数の分解
結合密度関数の分解:
25
同時密度関数の分解（続）
モデルを仮定（マルコフ性）
：システムモデル
：観測モデル
26
状態空間モデル
一般状態空間モデル：
：システムモデル
：観測モデル
非線形・非ガウス型状態空間モデル：
：システムモデル
：非線形関数
：観測モデル
：白色雑音
線形・ガウス型状態空間モデル：
：システムモデル
：観測モデル
：行列
：白色ガウス雑音
27
講義に出てくる話題の関係
結合密度関数
グラフィカルモデル
分配法則
確率伝搬法
マルコフ性
状態空間モデル
非線形・非ガウス
線形・
ガウス
カルマンフィルタ
サンプリング法，モンテカルロ法
粒子フィルタ
マルコフ連鎖モンテカルロ法
28
参考文献
1. 
片山徹，新版応用カルマンフィルタ，朝倉書店，2000.
2. 
片山徹，非線形カルマンフィルタ，朝倉書店，2011.
3. 
A. Papoulis and S. U. Pillai, Probability, Random Variables and Stochastic Process,
McGraw-Hill, 2002.
4. 
C. M. Bishop, Pattern Recognition and Machine Leaning, Springer, 2006.

Download Report