1. No category

信号解析 第６回講義録
日時：２０１４年５月２６日
講義内容：最小自乗法とデータ曲線
担当者：情報知能工学科 小島史男
1
はじめに
前回はカルバックの情報量について説明し、データから真の統計モデルを推定する評価法について学習
しました。今週はまず前回の復習を兼ねて、統計モデルの構造を与えて、それに含まれるパラメータを
データから推定する手法について説明します。前回の講義で、最小自乗法が統計モデルとして正規型を
仮定したときの最尤推定と一致することを示しました。今週の講義では、データから最小自乗法により
回帰曲線を求める計算法について説明します。赤池情報量基準 (AIC) は回帰曲線を求めるメタ最適化手
法として有名です。これで信号解析前半が終了します。説明する事項がたくさんありますが、基本的事
項はしっかりおさえていきましょう。教科書第５章の内容です。
2
回帰モデルによるデータのあてはめ
説明に入る前に、まず例題をみてみましょう。次の図は阪神タイガースが優勝した年（今年もひょっと
すると？）の阪神電鉄の株価１年間の毎日の最終終値を時系列的にまとめたものです。青色の不規則な
変動が実際のデータ、ピンク色の曲線が最小自乗法により求めた回帰モデルによる曲線です。このよう
Figure 1: 阪神タイガースが優勝した年の株価変動
なことはどうしてできるのかが今回の講義の目的です。ではまず回帰モデルとは何かについて説明して
−1
いきます。いまデータ長 N の時系列信号 {y[n]}N
n=0 について、各時刻 n におけるデータ y[n] が m 個の
変数 {xn1 , xn2 , · · · , xnm の線形結合と不確定要素 [n] との和で与えられるとします。すなわち
y[n] =
m
∑
ai xni + [n]
(n = 0, 1, 2, · · · , N − 1)
i=1
のことを回帰モデル (regression model) と呼びます。回帰モデルでは y[n] を目的変数、{xni }m
i=1 をこの
データが生成する説明変数と呼びます。説明変数の個数 m は次数 (order) と呼び、説明変数の前につい
ている係数 {ai }m
i=1 を回帰係数 (regression coefficients) と呼びます。第２項の [n] は説明変数では記述
できない不規則変動をモデル化したもので、残差 (residual) と呼びます。ここで残差 [n] の統計モデル
は平均 0 分散 σ 2 の独立な正規分布で与えたとします。以上の設定のもとで、適当な次元 m と説明変数
{xni }m
i=1 を具体的に与え、パラメータベクトル
θ = (a1 , a2 , · · · , am , σ 2 )
を適当に与えれば回帰モデルが計算できます。残差 [n] が互いに独立な分布であることを考慮すれば、
尤度および対数尤度はそれぞれ
L(θ) =
l(θ) =
N∏
−1
n=0
N
−1
∑
g(y[n]|θ, xn1 , · · · , xnm )
log g(y[n]|θ, xn1 , · · · , xnm )
n=0
と表記できます。またこのとき統計モデルは平均 0 分散 σ 2 の正規分布としましたので、尤度関数、対数
尤度関数の確率密度は


1
1
√
exp − 2
 2σ
2πσ
g(y[n]|θ, xn1 , · · · , xnm ) =
1
1
log g(y[n]|θ, xn1 , · · · , xnm ) = − log 2πσ 2 − 2
2
2σ
(
m
∑
y[n] −
)2 

ai xni
i=1
(
y[n] −
m
∑

)2
ai xni
i=1
となります。したがって、対数尤度関数は次のようになります。
(
−1
m
∑
N
1 N∑
l(θ) = − log 2πσ 2 − 2
y[n] −
ai xni
2
2σ n=0
i=1
)2
この尤度関数を最大にするパラメータベクトル θ を求めることにより統計モデルのもとで最良な近似値
（最尤推定量）を求めることができます。分散 σ 2 に関して、上記の対数尤度を偏微分してみましょう。
(
−
m
∑
∂l(θ)
N
1 N∑1
=
−
+
y[n]
−
ai xni
∂σ 2
2σ 2 2(σ 2 )2 n=0
i=1
この極値を求めると、
(
−1
m
∑
1 N∑
y[n] −
ai xni
σ
ˆ =
N n=0
i=1
)2
=0
)2
2
−1
m
となりますから、各時刻の時系列データ {y[n]}N
n=0 を上式に代入し、任意の回帰係数 {ai }i=1 を与えれ
2
2
ば分散の推定値が求まります。一方 σ = σ
ˆ を対数尤度関数に代入すると、対数尤度は
l(a1 , a2 , · · · , am ) = −
N
N
log 2πˆ
σ2 −
2
2
で与えられますので、これを {ai }m
ˆ 2 を最小にする、す
i=1 について最大にするには、分散の最尤推定量 σ
なわち以下の最小自乗問題を解くことになります。
(
−1
m
∑
1 N∑
y[n] −
ai xni
N n=0
i=1
3
)2
ハウスホルダー変換による計算アルゴリズム
前節では、回帰係数の最尤推定量は残差の最小自乗法解であることを示しました。ここでは、その解法
アルゴリズムについて簡単に説明しましょう。まず回帰モデル全体を行列でまとめていきましょう。N
次元ベクトル y, 、m 次元ベクトル a, N × m 行列 Z を



y=


y[0]
y[1]
..
.
y[N − 1]







=




[0]
[1]
..
.
[N − 1]







a=




a1
a2
..
.
am







Z=




x01
x11
..
.
···
···
x0m
x1m
..
.
···
xN −1,1 · · · xN −1,m






と表記すると、回帰モデルは簡単に
y = Za + と記述できます。つまり最尤推定量はベクトル a について以下の N 次元ユークリッドノルムの最小解を
求めることと同等です。
min kk2N = min ky − Zak2N
a
a
これを直接求めてもよいのですが、計算を効率よく実行するために直交行列 U を用いて以下のように変
換します。
ky − Zak2N = kU (y − Za)k2N = kUy − UZak2N
このような変換をしてもノルムの値は変化しませんから、変換後の最小化問題を解いてもよいことにな
ります。この最小解は a に関する線形代数方程式
(UZ)T UZa = (UZ)T Uy
を解けば良いことになります。UZ を解きやすい形に変形できればアルゴリズムは簡単になります。そ
こで考えられたのが、ハウスホルダー変換 (Householder transformation) です。説明変数の行列 Z の右
側に目的変数ベクトルを加えた拡大行列を
X = [Z | y]
のように定義します。この行列 X に適当なハウスホルダー変換を施して (N × (m + 1)) の上三角行列

s11 · · ·
 ..
..
 .
.

 0


S = UX = [UZ | Uy] =  0

 0

 ..
 .
0
s1m
..
.

s1,m+1
..
.












· · · smm sm,m+1
···
0
sm+1,m+1
···
0
0
..
..
..
.
.
.
···
0
0
に行変換します。これを利用すると
kUy − UZak2N
=
=

s1,m+1

..

.

 s
 m,m+1


0

..


.
0


s11

  .
  ..
 
 
− 0
  .
  .
  .

0

 
s1,m+1
s11



.
.
..

 −  ..
s
0
m,m+1
· · · s1m
..
..
.
.
··· 0
..
..
.
.
··· 0



a1

  ..
 .

 a
m

2



N

a1
· · · s1m
..   ..
..
.
.  .
am
· · · smm
2

2
 + sm+1,m+1
m
と変形できて、この式から最小解は下記の線形代数方程式を解けば良いことになります。





s1,m+1
a1
s11 · · · s1m

 ..
..
..   ..  = 
.
.

 .
.
.
.  .  
sm,m+1
am
0 · · · smm
後退代入をもちいることにより、回帰係数は後ろから順次求めていくことができます。以上の計算アル
ゴリズムをまとめてみましょう。
sm,m+1
a
ˆm =
smm
si,m+1 − si,i+1 a
ˆi+1 − · · · − sim a
ˆm
a
ˆi =
(i = m − 1, m − 2, · · · , 2, 1)
sii
また残差 の分散 σ 2 の推定量は
2
σ
ˆm
=
sm+1,m+1
N
によって計算できることになります。
4
赤池情報量基準 (AIC)
前節までの議論で最初の例題を解くほとんどの準備が整いました。最後にひとつだけ問題があります。
回帰モデルの次元 m はどう決めたらいいのでしょうか。データが多ければ大きくとればよいのでしょう
か。もしこの考えでいくと、前節の行列のサイズがどんどん大きくなり、そのうちに実現不可能な大き
さになってしまいます。できれば小さいほうが望ましいですね。このパラメータの次元決定に関しては、
有名な赤池の情報量基準（Akaike Information Criteria, 以降 AIC 基準と呼ぶ）があります。この方法
は実際には計算できない平均対数尤度と時系列データにより計算した対数尤度とのバイアス値（偏り）
を最適化するために求められた指標です。具体的な計算は省略しますが、AIC 基準とは
ˆ + 2k
AIC = −2l(θ)
をモデル選択の基準として用いる方法です。ここで右辺最初の項は最大対数尤度ですが、２項目の k は
パラメータ θ の次元です。AIC を最小にするモデルを選択することにより、できるだけ少ないパラメー
タ数で不規則信号の背景にある統計量を合理的に設計できる手法であり、簡便な手法であることから広
く用いられています。ではここでの問題に AIC を適用してみましょう。前節の計算アルゴリズムにより
最大対数尤度は
N
2
ˆ = − N log 2πˆ
σm
−
l(θ)
2
2
であること、およびパラメータ θ の次元が m + 1 であることに注意すれば、この回帰モデルの赤池情報
量基準は
(
)
2
AIC(m) = N log 2πˆ
σm
+ 1 + 2 (m + 1)
となります。回帰モデルの次元は小さいほどよいのですから、実際は適当な最大次元 M （通常は 20 前
後に設定）をさだめ、m ≤ M として m = 1 から順番に次元をあげていって AIC(m) が最小値となった
ところで計算を止め、そのときの次元を最適な m として、そのときのみ線形代数方程式を解き回帰モデ
ルをつくります。
5
R を使ったおさらい
では最後に、今回の講義で最初に示した回帰曲線のあてはめを R を使って実行しましょう。ここでは教
科書の５章で実際に解いた例を参考にします。回帰モデルを
y[n] = a +
ls
∑
bi sin(iωn) +
i=1
lc
∑
ci cos(iωn) + [n]
i=1
で与えます。教科書の添え字では回帰曲線の次元 m と重なりますのでここでは ls および lc を使うこと
にします。回帰曲線の次数 m に対して
[ ]
m
ls =
2
{
lc =
ls
(m : odd)
ls − 1 (m : even)
とします。説明変数およびパラメータは m が奇数のときは
{1, sin(ωn), cos(ωn), · · · , sin(ls ωn), cos(lc ωn) }
{a, b1 , c1 , · · · , bls , clc }
のように、また偶数の時は
{1, sin(ωn), cos(ωn), · · · , sin(ls ωn), cos(lc ωn) }
{a, b1 , c1 , · · · , bls , clc }
で与えられることになります。また説明変数に含まれる ω ですが、およそ 1 年の周期変動を考慮して
ω = 2π/365 と設定します。それではもとの時系列信号のシミュレーションデータを作成します。実際の
データの次元はわからないのですが、シミュレーションデータの作成に当たっては、周期関数に正規性白
色雑音が加法的に加わったとして作成します。そのコマンドラインとデータ作成のための関数 signal.R
を下記に示します。図 2 はその標本時系列信号の例です。
シミュレーションデータの生成
>
>
>
>
>
>
>
source("signal.R")
N <- 256
m <- 7
theta <- c(2, 1, 3, 2, 4, 7, 2)
y <- signal(N, m, theta)
plot(y, type="l", lty=1)
dev.copy2eps(file="simData.eps", width=10, height=6)
関数 signal.R のプログラム
signal <- function(n, m, theta) {
omega <- 2 * pi / 365;
ls <- m %/% 2;
if(m %% 2 == 0) lc <- ls - 1 else lc <- ls;
y <- numeric(n);
e <- rnorm(n, sd=1);
for (i in 1:n) {
y[i] <- theta[1] + e[i];
for (j in 1:ls) y[i] <- y[i] + theta[2*j] * sin(j * omega * i);
for (j in 1:lc) y[i] <- y[i] + theta[2*j+1] * cos(j * omega * i);
}
return(y);
}
5
0
−5
−10
y
10
15
20
このデータの回帰曲線を教科書 5.3 節に従って求めてみます。考えられる最高次数を mmax = 21 としま
0
50
100
150
200
Index
Figure 2: 試験データの標本時系列信号
250
す。これから教科書 66 ページにある説明変数行列 Z に目的変数ベクトル y を付加した (N ×(mmax+1) =
256 × 22) の行列 X を求め、上三角行列 S を求め、赤池情報量基準 AICj (j = 0, · · · , m) の評価を実施
します。ここでは計算の過程で R に実装されている QR 変換を求めました。これには 2 つの関数を用い
ました。また AIC の計算においては、変換された行列の (mmax + 1) 列目の情報のみ使いますので、こ
れをベクトル sm として抽出しておきます。その結果は図 3 のようになります。
AIC の計算
>
>
>
>
>
>
>
>
>
source("setx.R")
source("aic.R")
mmax <- 21
x <- setx(N, mmax, y)
s <- qr.R(qr(x, tol=1e-07))
sm <- s[,m+1]
aicv <- aic(N, mmax, sm)
barplot(aicv)
dev.copy2eps(file="aic.eps", width=10, height=6)
行列 X の作成 (setx.R)
setx <- function(N, m, y) {
ls <- m %/% 2;
if(m %% 2 == 0) lc <- ls - 1 else lc <- ls;
omega <- 2 * pi / 365.0;
x <- matrix(0, N, m+1);
for (i in 1:N) {
x[i,1] <- 1.0;
x[i, m+1] <- y[i];
for (j in 1:ls) x[i,j*2] <- sin(omega * i * j);
for (j in 1:lc) x[i,j*2+1] <- cos(omega * i * j);
}
return(x);
}
AIC の評価 (aic.R)
aic <- function(N, m, sm) {
sigma2 <- numeric(m+1);
aicv <- numeric(m+1);
sigma2[m+1] <- sm[m+1]^2 / N;
aicv[m+1] <- N * log(2*pi*sigma2[m+1] + 1) + 2 * (m + 1);
for (j in m:1) {
sigma2[j] <- sm[j]^2 / N + sigma2[j+1];
aicv[j] <- N * log(2*pi*sigma2[j] + 1) + 2 * j;
}
return(aicv);
}
0
50
100
150
200
250
300
350
AIC の値が最小となるの次元は図 3 から mopt = 9 であることがわかります。またコマンドラインで
Figure 3: AIC の評価
which.min(aicv) または order(aicv) でも確認することができます。この次元における回帰係数を求めて
みます。この回帰係数から回帰曲線を求め、元の観測データと比較したのが図 4 です。十分な近似が得
られていることがわかります。
回帰曲線
> source("regression.R")
> which.min(aicv)
[1] 9
> order(aicv)
[1] 9 10 11 12 13 14 15 16 17 18 19 20 21 22 8 7 6 1
> mopt <- 9
> thetaHat <- regression(mmax, mopt, s)
> thetaHat
[1] 0.5864402 3.1175653 1.4338817 3.9043239 4.5895112
[8] -0.2813638 0.4618989
>
>
+
>
>
>
2
4
5
7.3056253
3
3.0952581
source("regCurve.R")
curve(regCurve(mopt, thetaHat, x), 0, 256,
xlim = c(0, 256), ylim = c(-10, 20), xlab="", ylab="")
par(new=T)
plot(y, xlim = c(0, 256), ylim = c(-10, 20))
dev.copy2eps(file="comparison.eps", width=10, height=6)
回帰係数を求める (regression.R)
regression <- function(m, j, s) {
ahat <- numeric(j);
ahat[j] <- s[j, m+1] / s[j, j];
if (j > 1) {
for (i in (j-1):1) {
ahat[i] <- s[i, m+1] / s[i, i];
for (k in (i+1):j) {
ahat[i] <- ahat[i] - s[i, k] * ahat[k] / s[i, i];
}
}
}
return(ahat);
}
回帰曲線 (regCurve.R)
regCurve <- function(m, theta, x) {
omega <- 2 * pi / 365;
y <- theta[1];
ls <- m %/% 2;
if(m %% 2 == 0) lc <- ls - 1 else lc <- ls;
if (m > 0) {
for (j in 1:ls) y <- y + theta[2*j] * sin(j * omega * x);
for (j in 1:lc) y <- y + theta[2*j+1] * cos(j * omega * x);
}
return(y);
}
y
−10
−5
0
5
10
15
20
今週学んだ重要なことは、観測信号に対する回帰曲線の近似精度が赤池情報量基準 AIC という数値で評
0
50
100
150
200
250
Index
Figure 4: 回帰曲線と観測信号の比較
価できることです。この結果から、回帰曲線の次元を高くとればよいというわけではないことがわかり
ます。
下記は宿題ではありませんが、興味のある受講生はチャレンジしてみてください。
課題： 数年前の阪神 HD の株価データを挙げておきますので、その回帰曲線を求めて見てください