7/7（月）

信号解析第１１回講義録
日時：２０１４年７月７日
講義内容：状態推定問題
担当者：情報知能工学科小島史男
1
はじめに
これまで状態空間モデルによって時系列信号を記述し、カルマンフィルタを構成する方法について学習し
てきました。カルマンフィルタはシステムの状態量の観測信号に関する平均過程および誤差分散過程を
逐次求めていく計算アルゴリズムです。これらの結果をどのように利用していくのかが本日の学習テー
マです。受講生のみなさんはカルマンフィルタの実行例のレポートを課しましたが、理解できましたで
しょうか。それを使いこなすことができましたでしょうか。今週の講義では、まず状態推定問題の枠組
みについて述べ、使用する観測データ量に応じて、予測 (prediction)、濾波 (filtering)、平滑 (smooting)
の３種類の推定問題について考えてみましょう。さらに状態空間モデルに関する尤度関数について考え、
それが ARMA 過程のモデリングに適用できることをしまします。最後に予測について R での実行例を
通じて考えていくことにします。一度に多くのことを学習しますので、よく復習しておいてください。
教科書では主に第９章の内容です。
2
状態推定問題とは
カルマンフィルタでは、現在時刻まで時々刻々入手される観測データ Yn = {˜
yi }ni=1 から、現在時刻 n にお
ける条件付き平均値ベクトル x
ˆn|n := E(x[n]|Yn ) および分散共分散行列 Vn|n := Cov(x[n] − x
ˆn|n , x[n] −
x
ˆn|n |Yn ) を求めました。観測値 Yj = {˜
yi }ji=1 が与えられたとき、状態 x[n] の推定値 x
ˆn|j := E(x[n]|Yj )
を求める問題を状態推定問題と呼びます。特に j < n の場合、観測区間より先の状態を推定することに
なり、これを予測問題 (prediction) といいます。また逆に j > n の場合は過去の状態を推定するという
ことで、平滑化問題 (smoothing) と呼びます。カルマンフィルタを解くことは j = n に相当しますが、
この問題はろ波問題 (filtering) ともいいます。つまり状態推定問題は観測データの条件によってろ波問
題、平滑化問題、予測問題と３つの問題があります。本日の講義では、先週学習したカルマンフィルタ
（つまりろ波問題）の結果を用いて平滑化および予測の問題を解く方法について学習します。
3
平滑化問題
カルマンフィルタでは、観測データが入ってくるごとに、１時刻先の予測とフィルタの計算を繰り返し
ました。しかしながら観測データを一時的にバッファに格納しておき、そのデータを用いて過去のデー
タを推定することも考えられます。このような問題を平滑化問題と呼びます。いま固定区間の観測デー
タ YN = {˜
yi } N
i=1 が取得されたとします。このとき事前にカルマンフィルタの結果
{
}
x
ˆn|n−1 , x
ˆn|n , Vn|n−1 , Vn|n
(n = 1, 2, · · · , N )
は計算されているものとします。このとき、各時間ステップのカルマンフィルタの１期先予測およびフィ
ルタリング結果を使って平滑化の計算を以下のように実行することができます。
−1
A[n] = Vn|n F [n + 1]T Vn+1|n
(
)
(
)
x
ˆn|N
= x
ˆn|n + A[n] x
ˆn+1|N − x
ˆn+1|n
Vn|N
= Vn|n + A[n] Vn+1|N − Vn+1|n A[n]T
この計算アルゴリズムではまずカルマンフィルタを観測データの上限 N まで解き、そこから逆向きに上
の式を計算していけばよいことがわかります。したがって、この蹴り返し演算の初期値はカルマンフィ
ルタの最後の状態
{
}
x
ˆN |N , VN |N
で与えます。
4
状態の長期予測
カルマンフィルタでは常にワンステップの予測を行います。これを繰り返し適用することで、現在時刻
n の観測データ Yn = {˜
yi }ni=1 から j ステップ先の状態 x[n + j] (j > 1) の予測を行うことができます。観
測データの集合を
Yn = {˜
y1 , y˜2 , · · · , y˜n }
Yn+1 = {˜
y1 , y˜2 , · · · , y˜n , 0}
Yn+2 = {˜
y1 , y˜2 , · · · , y˜n , 0, 0}
···
···
と考えれば、形式的に Yn+j = Yn (j > 1) が成立します。このとき明らかに、
E (x[n + j]|Yn+j ) = E (x[n + j]|Yn )
が成立するので、カルマンフィルタのワンステップ予測のみを j 回反復すれば予測をおこなうことがで
きます。すなわち状態の長期予測の計算アルゴリムは、
x
ˆn+i|n = F [n + i]ˆ
xn+i−1|n
Vn+i|n = F [n + i]Vn+i−1|n F [n + i]T + G[n + i]Q[n + i]G[ n + i]T
(i = 1, 2, · · · , j)
のように構成することができます。この繰り返し演算の初期値も平滑化のときと同じで、カルマンフィ
ルタの最後の状態 {ˆ
xN |N , VN |N } で与えます。
5
時系列信号の予測
前節の結果にもとづいて、時系列信号の予測の方法についてさらに考えていきましょう。状態空間モデ
ルでは時系列信号 y[n] は次のようなモデルで与えられます。
y[n] = H[n]x[n] + w[n] (l × 1)
観測信号の現在時刻 n までの情報 Yn = {˜
yi }ni=1 に関する j ステップ先の予測平均値を yˆn+j|n := E(y[n +
j]|Yn ) また予測分散共分散行列を dˆn+j|n := Cov(y[n + j], y[n + j]|Yn ) で表すことにすると、次の方程
式が得られます。
yˆn+j|n = E(H[n + j]x[n + j] + w[n + j]|Yn )
= H[n + j]ˆ
xn+j|n
dˆn+j|n = Cov(H[n + j]x[n + j] + w[n + j], H[n + j]x[n + j] + w[n + j]|Yn )
= H[n + j]Cov(x[n + j], x[n + j]|Yn )H[n + j]T + H[n + j]Cov(x[n + j], w[n + j]|Yn )
+Cov(w[n + j], x[n + j]|Yn )H[n + j]T + Cov(w[n + j], w[n + j]|Yn )
= H[n + j]Vn+j|n H[n + j]T + R[n + j]
y[n + j] の観測値の集合 Yn に関する予測分布は平均値 yˆn+j|n および分散共分散行列 dˆn+j|n を持つ正規
分布で与えられます。以上のように、予測を行うには
(1) 一定のデータ取得を通じて時系列モデルを構成し、状態空間モデルを記述する
(2) カルマンフィルタを解いて、現在時刻における平均値と分散共分散を求める
(3) 長期予測モデルを解く。
(4) これをもとに時系列予測の平均値と標準偏差を各ステップで求める。
6
状態空間モデルと尤度計算
状態空間モデルの尤度関数についてはここまで詳しく説明していません。そこで改めて以下にこのこと
について考察を進めることにします。状態空間モデルの尤度関数は条件付き確率密度関数の積によって
L(θ) = gN (yN |˜
y1 , y˜2 , · · · , y˜N −1 , θ) gN −1 (yN −1 |˜
y1 , y˜2 , · · · , y˜N −2 , θ) · · · g1 (y2 |˜
y1 , θ)
=
=
N
∏
n=1
N
∏
gn (yn |˜
y1 , y˜2 , · · · , y˜n−1 , θ)
gn (yn |Yn−1 , θ)
n=1
によって与えられることは先週学習しました。ここでパラメータ θ は推定すべき母数（パラメータ）の
ベクトルです。たとえば AR(m) モデルを構成するときは
[
θ = a1 , a 2 , · · · , a m , σ 2
]
がこれに相当します。ところで観測値の集合 Yn−1 = {˜
yi }n−1
i=1 による y[n] の予測分布 g(y[n]|Yn−1 , θ) は
どうなるのでしょうか。これはカルマンフィルタの結果から直接求めることができます。まず予測誤差
の平均 yˆn|n−1 := E(y[n]|Yn−1 ) と分散共分散行列 dˆn|n−1 := Cov(y[n], y[n]|Yn−1 ) は観測モデルとワンス
テップ予測式より、
yˆn|n−1 = E(y[n]|Yn−1 ) = E(H[n]x[n] + w[n]|Yn−1 ) = H[n]ˆ
xn|n−1
dˆn|n−1 = Cov(H[n]x[n] + w[n], H[n]x[n] + w[n]|Yn−1 )
= H[n]Cov(x[n], x[n]|Yn−1 )H[n]T + H[n]Cov(x[n], w[n]|Yn−1 )
+Cov(w[n], x[n])H[n]T + Cov(w[n], w[n]|Yn−1 )
= H[n]Vn|n−1 H[n]T + R[n]
が成立します。予測分布 gn (y[n]|Yn−1 , θ) は平均 yˆn|n−1 および分散共分散行列 dˆn|n−1 に従う正規型確率
密度関数で与えられるので、予測分布は
(
gn (y[n]|Yn−1 , θ) =
1
√
2π
)l
dˆn|n−1
− 2l
{
exp −
)T (
)−1 (
)}
1(
y[n] − yˆn|n−1
dˆn|n−1
y[n] − yˆn|n−1
2
となります。ここで l は観測モデルの次元です。この式を尤度関数に順次代入すれば状態空間モデルの
対数尤度は結局
[
]
N
N {
}T {
}−1 {
}
∑
∑
1
log dˆn|n−1 +
y˜n − yˆn|n−1 (θ)
dˆn|n−1 (θ)
y˜n − yˆn|n−1 (θ)
l(θ) = − lN log(2π) +
2
n=1
n=1
によって求めることができます。一般的な状態空間モデルにおいても上記の対数尤度関数を最大にする
未知パラメータ θ をもとめることは可能となります。したがって、ARMA モデルの係数同定にはこの計
算手法を適用することが可能です。しかしながら、この最適化問題は θ に関して非線形となるので、推
定計算を更新するためにカルマンフィルタを何度も解かねばならず、計算面で実用的とはいえません。
すでに先週の宿題で行ったように、時系列信号がＡＲモデルで与えられる場合には、事前に適当なＡＲ
次元でユールウオーカーの方法で標本自己共分散関数を用いてＡＲ係数の決定をおこない、ＡＩＣ法で
最適な次元を定めた上で、カルマンフィルタの計算が確実に実行できます。したがって、この結果を用
いれば、今回の平滑化や長期予測を実行することができます。カルマンフィルタ、および平滑化、予測
法のアルゴリズムを理解するためには、適当なＡＲモデルでその精度を確認することだとおもいます。
7
R で長期予測を実現する
簡単のため状態空間モデルの行列はすべて定数行列として話をすすめていきます。このとき長期予測の
計算アルゴリズムは
x
ˆn+i|n = F x
ˆn+i−1|n
Vn+i|n = F Vn+i−1|n F T + GQGT
yˆn+i|n = H x
ˆn+i|n
dˆn+i|n = HVn+i|n H T + R
(i = 1, 2, · · · , j)
となります。これを起動するにはカルマンフィルタの現在時刻 n の情報 {ˆ
xn|n , Vn|n } が必要となります
ので、前回のカルマンフィルタのプログラムを少し書き加えておきます。なお今回はファイルに書き込
みを行うためのプログラムが長くなりましたので、それぞれ２つに分けて記述しています。
✓
カルマンフィルタの書き換え（その１）
✏
# dim(data) = n; dim(system) = k; dim(noise) = m; dim(observation) = l
kalmanrev <- function(n, k, m, l, x0, xs0, fm, gm, hm, qm, rm, y) {
# Main part
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);
km <- matrix(0, nrow=k, ncol=l);
xh_nm1_nm1 <- x0;
vm_nm1_nm1 <- xs0 %*% t(xs0) - x0 %*% t(x0);
ym <- matrix(0, nrow=n, ncol=l);
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;
# Iteration
xh_nm1_nm1 <- xh_n_n;
vm_nm1_nm1 <- vm_n_n;
ym[t,] <- hm %*% xh_n_n;
}
# Save computational results
}
✒
✑
✓
カルマンフィルタの書き換え（その２）
✏
kalmanrev <- function(n, k, m, l, x0, xs0, fm, gm, hm, qm, rm, y) {
# Main part
# Save computtional results
if (k < 2) {
cat("k becomes less than 2\n");
return();
}
out <- file("/Users/mean.dat", "w");
writeLines(paste("xh_nm1_nm1 <- c"), out, sep="\n");
for (i in 1:(k-1)) {
writeLines(paste(xh_n_n[i]), out, sep="");
writeLines(paste(","), out, sep="\n");
}
writeLines(paste(xh_n_n[k]), out, sep="");
writeLines(paste(")"), out, sep="\n");
close(out);
out <- file("/Users/covariance.dat", "w");
writeLines(paste("vm_nm1_nm1 <- c"), out, sep="\n");
for (i in 1:(k-1)) {
for (j in 1:k) {
writeLines(paste(vm_n_n[i,j]), out, sep="");
writeLines(paste(","), out, sep="\n");
}
}
for (j in 1:(k-1)) {
writeLines(paste(vm_n_n[k,j]), out, sep="");
writeLines(paste(","), out, sep="\n");
}
writeLines(paste(vm_n_n[k,k]), out, sep="");
writeLines(paste(")"), out, sep="\n");
close(out);
return(ym);
}
✒
✑
✓
長期予測の計算アルゴリズム（その１）
✏
# dim(data) = n; dim(system) = k; dim(noise) = m; dim(observation) = l
prediction <- function(j, k, m, l, kp, xh_n_n, vm_n_n, fm, gm, hm, qm, rm) {
# Main part
xp_npim1_n <- numeric(k);
xp_npi_n <- numeric(k);
vp_npim1_n <- matrix(0, k, k);
vp_npi_n <- matrix(0, k, k);
yp <- numeric(l);
dh <- matrix(0, nrow=l, ncol=l);
xhp <- numeric(j);
dhp <- numeric(j);
xp_npim1_n <- xh_n_n;
vp_npim1_n <- matrix(vm_n_n, nrow=k, ncol=k);
listt = 1:j;
for(t in listt) {
# prediction
xp_npi_n <- fm %*% xp_npim1_n;
vp_npi_n <- fm %*% vp_npim1_n %*% t(fm) + gm %*% qm %*% t(gm);
# output
yp <- hm %*% xp_npi_n;
dh <- hm %*% vp_npi_n %*% t(hm) + rm;
xhp[t] <- yp[kp];
dhp[t] <- dh[kp,kp];
# Iteration
xp_npim1_n <- xp_npi_n;
vp_npim1_n <- vp_npi_n;
}
# Save computational results
}
✒
✑
✓
長期予測の計算アルゴリズム（その２）
✏
prediction <- function(j, k, m, l, kp, xh_n_n, vm_n_n, fm, gm, hm, qm, rm) {
# Main part
# Save computational results
if (j < 2) {
cat("j becomes less than 2\n");
return();
}
out <- file("/Users/yhatp.dat", "w");
writeLines(paste("yhp <- c"), out, sep="\n");
for (i in 1:(j-1)) {
writeLines(paste(xhp[i]), out, sep="");
writeLines(paste(","), out, sep="\n");
}
writeLines(paste(xhp[j]), out, sep="");
writeLines(paste(")"), out, sep="\n");
close(out);
out <- file("/Users/yphigh.dat", "w");
writeLines(paste("yphigh <- c"), out, sep="\n");
for (i in 1:(j-1)) {
writeLines(paste(xhp[i]+sqrt(dhp[i])), out, sep="");
writeLines(paste(","), out, sep="\n");
}
writeLines(paste(xhp[i]+sqrt(dhp[j])), out, sep="");
writeLines(paste(")"), out, sep="\n");
close(out);
out <- file("/Users/yplow.dat", "w");
writeLines(paste("yplow <- c"), out, sep="\n");
for (i in 1:(j-1)) {
writeLines(paste(xhp[i]-sqrt(dhp[i])), out, sep="");
writeLines(paste(","), out, sep="\n");
}
writeLines(paste(xhp[i]-sqrt(dhp[j])), out, sep="");
writeLines(paste(")"), out, sep="\n");
close(out);
return(ym);
}
✒
このプログラムを実行することにより、予測を実現します。それぞれのプログラムは（その１）および
（その２）で一つのプログラムとなっていますので、注意してください。また（その２）ではファイルに
結果を格納していますがその格納先はそれぞれ適当に変更してください。W indows 系のプログラムで
は、ファイルの格納先は C : を頭につける必要があります。
まず実験に用いるデータを準備しておきます。今回はデータ長 256 の時系列信号のうち最初の 128 個
のデータを現在時刻までのデータ (n = 128) として、次の 36 ステップ先までの信号を予測することにし
ます。信号としては以前講義録で用いたデータ signal.dat を利用します。与えられたデータからレビン
✑
ソンのアルゴリズムで AR モデルを決定します。R のコンソールに一行ずつコマンドを書き込んでいく
のは結構面倒です。下記の図に示すように、ソースファイルの読み込み (> の記号のついたところ）以外
はあらかじめテキストファイルに書き込んでおく一括してカット・エンド・ペーストすれば一気にプロ
グラムが実行できます。ただしそのときには命令文の最後に必ず ; の記号を加えておきます。このプロ
グラムを実行すると mmax を最高次元とした AR モデルの係数情報が theta に出力されます。
✓
✏
AR モデリング
#AR modeling (Part I)
#
>
>
>
Read files
source("signal.dat")
source("acfv.R")
source("levinson.R")
# Use cut and paste
n <- 128;
mmax <- 64;
yt <- c(y[1:n]);
chat <- acfv(n,mmax,yt);
theta <- levinson(n,mmax,chat);
✒
✑
次にカルマンフィルタ（修正版）を用いて n = 196 までのフィルタを実行します。ここではコンソー
ルへ指令するコマンドが多いので、カット・エンド・ペーストを実行するとあっという間に実行できま
す。このプログラムを実行するとあとの予測で必要な初期推定値 mean.dat と covariance.dat という２
つのファイルがデフォルトのデレクトリに格納されます。
✓
✏
カルマンフィルタの実行
# Kalman filtering (Part II)
# Read file
> source("KalmanRev.R")
# Use cut and paste
k <- 32;
m <- 1;
l <- 1;
x0 <- numeric(k);
xs0 <- diag(k);
a <- theta[k+1,1:k];
ez <- cbind(diag(k-1), numeric(k-1));
fm <- rbind(a,ez);
sigma2 <- theta[k+1,mmax+1];
gm <- matrix(c(1,numeric(k-1)), nrow=k, ncol=1);
hm <- matrix(c(1,numeric(k-1)), nrow=1, ncol=k);
qm <- matrix(c(sigma2), nrow=1,ncol=1);
rm <- matrix(0, nrow=1,ncol=1);
yh <- kalmanrev(n, k, m, l, x0, xs0, fm, gm, hm, qm, rm, y);
✒
最後に予測を行います。実行ファイル prediction.R とカルマンフィルタで生成された２つのデータ
ファイル mean.dat と covariance.dat を読み込み、長期予測を行います。このプログラムを実行すると
時系列の予測情報 yphat.dat および標準予測誤差区間の上限 yphigh.dat、下限 yplow.dat がファイルに
格納されます。
✑
✓
長期予測の実現
✏
# Prediction (Part III)
#
>
>
>
Read files
source("prediction.R")
source("mean.dat")
source("covariance.dat")
# Use cut and paste
j <- 36;
kp <- 1;
xhp <- prediction(j, k, m, l, kp, xh_n_n, vm_n_n, fm, gm, hm, qm, rm);
✒
✑
以上の結果をグラフ表示し、その結果をファイルに格納します
✓
結果の表示
✏
#Results
#
>
>
>
Read files
source("yhatp.dat")
source("yphigh.dat")
source("yplow.dat")
# Cut and Paste
plot(c(y[n:(n+j)]), ylim=c(-5, 5), type="l", lty=1, ann=F, axes=F);
par(new=T);
plot(c(xh_n_n[1], yhp), ylim=c(-5, 5), type="l", lty=2, ann=F, axes=F);
par(new=T);
plot(c(xh_n_n[1], yphigh), ylim=c(-5, 5), type="l", lty=3, ann=F, axes=F);
par(new=T);
plot(c(xh_n_n[1], yplow), ylim=c(-5, 5), type="l", lty=4, ann=F);
dev.copy2eps(file="prediction.eps", width=10, height=6);
✒
さてこれらの結果を AR モデルの次元ごとに表示してみると以下のようになりました。図にあると
おり AR モデルの次元が低い m = 16 の場合はほとんど追従しません。m = 32 ぐらいまでくるとすこし
近づいてきますがまだ十分とはいえません。しかし m = 64 になるとプロファイルは追跡できているこ
とがわかります。このようなことはなぜ起こるのか、これは長期予測のモデルが単に１期先予測を繰り
返し行っているだけでしかないからです。つまり予測式は各時刻にシステムの入力はありません。そこ
で１期先予測を繰り返し実行すれば減衰していきます。この状態はカルマンフィルタの初期状態で決定
されます。したがって次元が高くないとシステムへの追従性はなくなっていきます。このように予測に
はモデリングの精度が問われます。
✑
4
2
0
−2
−4
0
5
10
15
20
25
30
35
30
35
−4
−2
0
2
4
Figure 1: 長期予測 (m = 16)
0
5
10
15
20
25
Figure 2: 長期予測 (m = 32)
4
2
0
−2
−4
0
5
10
15
20
25
Figure 3: 長期予測 (m = 64)
30
35

Download Report