5/19（月）

信号解析第５回講義録
日時：２０１４年５月１９日
講義内容：統計モデルの評価法
担当者：情報知能工学科小島史男
1
はじめに
先週は統計処理の基本となる乱数の生成法から確率分布とその統計モデルについて学習しました。とこ
ろで、不規則信号を取得したとしても、その真の分布は実際にはわかりません。データがあたえられた
とき、どの統計モデルを採用すればよいのでしょうか。今週の講義では与えられたデータから真の統計
モデルを推定する評価法について学習していきます。教科書第４章の 4.2 節以降の内容です。
2
カルバックの情報量とエントロピー最大化原理
それでは、まず真のモデルを g(y)、なんらかの方法で近似した統計モデルを f (y) としましょう。この真の
モデルと統計モデルとの近さを測る評価法として以下のカルバック・ライブラーの情報量 (Kullback-Leibler
information) があります。
[
I(g; f ) = EY
{
g(Y )
log
f (Y )
}]
∫
∞
{
}
g(y)
=
log
g(y)dy
f (y)
−∞
教科書の記法に従ってこれを以降 KL 情報量と呼ぶことにします。すこし天下り的ですが、log x ≤ x − 1
を利用すると
∫
{
∞
}
∫
{
}
∞
g(y)
f (y)
g(y)dy = −
log
g(y)dy
f (y)
g(y)
−∞
−∞
}
∫ ∞
∫ ∞
∫ ∞ {
f (y)
1−
g(y)dy =
g(y)dy −
f (y)dy = 1 − 1 = 0
≥
g(y)
−∞
−∞
−∞
log
となるので、結局 KL 情報量は
I(g; y) ≥ 0
となり、また、f (y) = g(y) のときに最小値 0 を持つことになります。したがって、この情報量が小さ
いほど近似度が良いと判定できることになるのです。また KL 情報量 I(g; y) にマイナス符号を付けた評
価 B(g; f ) = −I(g; f ) を一般化されたエントロピー (entropy) と呼びます。この場合、もっとも良い近
似とはエントロピーを最大化することになります。以上、この評価のもとで、最適な近似モデルを構成
することをエントロピー最大化原理 (entropy mazimaization principle) と呼びます。実際に取得された
有限長の信号から近似モデルを求める方法については次節以降で説明することにして、まずこのエント
ロピーの計算をしてみましょう。真のモデル g(y) が、平均値 µ、分散 σ 2 の正規型確率密度関数
{
1
(y − µ)2
g(y|µ, σ ) := √
exp −
2σ 2
2πσ
}
2
に従うと仮定します。一方近似モデル f (y) が、同じ正規型分布だが異なる母数、つまり平均値 ξ 、分散
τ 2 をもつ、
}
{
2
1
(y
−
ξ)
f (y|ξ, τ 2 ) := √
exp −
2τ 2
2πτ
で構成すると考え、このときの KL 情報量を計算してみましょう。まず
∫
I(g; y) =
∞
−∞
{log g(y) − log f (y)} g(y)dy
ですので、被積分項の対数を計算すると
1
(y − µ)2 1
(y − ξ)2
log g(y) − log f (y) = − log 2πσ 2 −
+ log 2πτ 2 +
2
2
2σ
2
2τ 2
{
}
1
τ2
(y − µ)2 (y − ξ)2
=
log 2 −
+
2
σ
σ2
τ2
となります。項別に積分を実行すると最初の２項は、
∫
∞
τ2
τ2
g(y)dy
=
log
σ2
σ2
−∞
∫ ∞
∫ ∞
2
1
(y − µ)
1
g(y)dy =
(y − µ)2 g(y)dy = 2 × σ 2 = 1
2
2
σ
σ −∞
σ
−∞
log
となり、また３番目の項は
(y − ξ)2 = {(y − µ) + (µ − ξ)}2 = (y − µ)2 + 2(µ − ξ)(y − µ) + (µ − ξ)2
と変形して、それぞれの積分を計算すると
∫
2(µ − ξ)
∞
(y − µ)2 g(y)dy = σ 2
−∞
∞
∫
−∞
(y − µ)g(y)dy = 0
(µ − ξ)2
∫
∞
−∞
g(y)dy = (µ − ξ)2
となるので、結局 KL 情報量は
1
I(g; f ) =
2
{
τ2
σ 2 + (µ − ξ)2
log 2 − 1 +
σ
τ2
}
ということになります。もちろん、以上のことは形式的な式変形ですので、平均値、分散量等具体的に
数値をいれて、いろいろ確かめてみましょう。
3
時系列信号による推定 KL 情報量の計算
前節での KL 情報量は真の分布 g(y) が未知であるので、現実には計算できないことになります。そこで
−1
真の分布の代替えとして、時系列信号 {y[n]}N
n=0 から統計モデル f (y) の KL 情報量を推定する方法に
ついて説明していきます。KL 情報量は
[{
I(g; f ) = EY
log
g(y)
f (y)
}]
= EY {log g(Y )} − EY {log f (Y )}
と分解できます。この第１項は近似モデル f (Y ) に関して不変ですから、無視することができます。す
なわち右辺の第２項を f に関して最大にすれば最良の近似になるはずです。この第２項
EY {log f (Y )} =
∫
∞
−∞
{log f (Y )}g(y)dy
のことを平均対数尤度 (expected log-likelihood) と呼びます。もちろんこの項にも真の分布 g(y) が含ま
れていますので、結局計算はできません。ところで第２回目の講義で説明した標本自己共分散関数の推
−1
定計算のように、時系列信号 {y[n]}N
n=0 によって推定することは可能です。いまこの時系列信号が統計
的に独立であるならば、
−1
1 N∑
log f (y[n]) −→ EY {log f (Y )} as N → ∞
N n=0
が成立することが大数の法則によって保証されています。一般の時系列信号は相互に関連性をもってい
N −1
ますので、この収束性の保証はありませんが、以上の議論の拡張として、時系列信号 YN = {y[n]}n=0
による統計モデルの計算量
{ ∑
N −1
log f (y[n])
log f (y[0], y[1], · · · , y[N − 1])
n=0
l(YN ) =
を対数尤度 (log-likelihood)、またその指数、すなわち
{ ∏
N −1
f (y[n])
f (y[0], y[1], · · · , y[N − 1])
n=0
L(YN ) =
を尤度 (likelihood) と呼び、統計モデルの評価指標に用いられています。上式においては、上部は時系
列信号が独立に観測された場合を、また下部の f は N 個の時系列データによる N 次元の同時確率密度
関数を意味しています。
4
最大尤度法とパラメータ推定
尤度（対数尤度）は時系列信号の背景にある統計量を評価するもっとも重要な指標であることを説明し
ました。尤度関数はよく知られているこれまでのいろいろな計算法に明確な解釈を与えることができま
す。たとえば最小自乗法は誰でも知っている統計処理法ですが、その平均の意味を尤度関数から考えて
みましょう。対数尤度を計算するには、統計モデル f (y) の構造を仮定して、その関数の係数をパラメー
タ θ としてこの値を変化させることにより、いろいろな尤度の値を調べる方法が一般的に用いられます。
すなわち対数尤度を θ の関数
{ ∏
N −1
l(θ) =
f (y[n]|θ)
f (y[0], y[1], · · · , y[N − 1]|θ)
n=0
として表記します。この統計モデル f (y|θ) をパラメトリックモデルと呼び、対数尤度 l(θ) を最大にする
パラメータを求める計算法を最尤法 (maximum likelihood method) と呼びます。以下で統計モデルが
平均値 θ 分散 1 の正規型確率密度関数
{
(y − θ)2
1
f (y|θ) := √ exp −
2
2π
}
で与えられるときに、平均値 θ を最尤推定により計算してみましょう。この場合の対数尤度は、前節で
おこなった計算を施すと
−1
1 N∑
N
l(θ) = − log 2π −
(y[n] − θ)2
2
2 n=0
となります。この対数尤度 l(θ) を最大とする θ は、上式において、θ に関する極値をとり、
−1
1 N∑
θˆ =
y[n]
N n=0
によって計算できます。時系列信号が独立に得られる場合は、データ長を無限大にとれば、その標本平
均の極限の分布が平均値 θ 分散 1 の正規型の統計モデルに一致することになります。このパラメータ最
適化（最大尤度法）は自乗誤差
N
−1
∑
(y[n] − θ)2
n=0
の最小化にほかなりません。これはよく知られている最小自乗法 (least square method) ですが、その最
適値である標本平均は以上の意味で合理的に説明できます。最尤法に関する計算アルゴリズムの詳細は
講義の目的ではありませんので、ここではこれ以上深入りしないことにします。
R を使ったおさらい
5
今回は、２節で計算した KL 情報量を図示してみましょう。簡単のため真の統計モデルは g(y|µ = 0, σ 2 = 1)
とします。このとき近似統計量 f (y|ξ, τ 2 ) をそれぞれ −2 ≤ ξ ≤ 2, 0.2 ≤ τ 2 ≤ 1.8 として、計算したと
きの分布を描きます。図 1 は情報量の分布を立体表示しています。実際、ξ = 0, τ 2 = 1 において、KL
情報量が最小値 0 を示していることが確認できます。
KL 情報量の分布
> xi <- seq(-2, 2, length=50)
> tau2 <- seq(0.2, 0.8, length=50)
> f <- function(xi, tau2) {0.5 * (log(tau2) - 1 + (1 + xi * xi)/tau2)}
> z <- outer(xi, tau2, f)
> persp(xi, tau2, z, theta=30, phi=30, expand=0.5, col=rainbow(50), border=NA)
>dev.copy2eps(file="kullback.eps")
tau
2
z
xi
Figure 1: カルバックの情報量分布
次に、対数尤度関数を計算して標本データの統計モデルの評価を行いましょう。次の計算は、まず
コーシー分布の乱数を標本データとして、対数尤度関数を求めてみます。統計モデル f (y) として正規分
布 dnorm(mean = 0, sd = 1) 、コーシー分布 dcauchy(location = 0, scale = 1) を仮定して対数尤度関
数を次のように定義します。
対数尤度関数の定義
lgauss <- function(n,y) {
va <- 0;
for(i in 1:n){
va <- va + log(dnorm(y[i], mean=0, sd=1));
}
return(va);
}
lcauchy <- function(n,y) {
va <- 0;
for(i in 1:n){
va <- va + log(dcauchy(y[i], location=0, scale=1));
}
return(va);
}
この関数を呼び出した計算結果は下記となりました。
対数尤度の評価値
>
> source("likelihood.R")
> y <- rcauchy(64, location=0, scale=1.0)
> y
> n <- 64
> lgauss(n,y)
[1] -500.7108
> lcauchy(n,y)
[1] -155.7048
>
この例でわかるように対数尤度関数の値はコーシー分布のほうが正規型分布のときより大きな値をとり、
コーシー分布で近似するのが良いという判定ができます。
練習問題：皆さんは正規性白色雑音を時系列として逆に上記の評価をしてみてください。

Download Report