【学習動機】
分析に用いているモデルの評価にAICが使えるのか不安になったので, AICについて学び直すことにした. 小西・北川『情報量規準』でAICの導出を追い直しつつ, Rで数値実験しているものを探していたら, ノンパラでよく見かける竹澤先生の書いた『AICをぶっ壊す』というものに出会った.
http://cse.naro.affrc.go.jp/takezawa/AIC%E3%82%92%E3%81%B6%E3%81%A3%E5%A3%8A%E3%81%99.pdf
このドキュメントを読んで, コードがもっと速くできると思ったので, コードを書き換えつつ, AICについて学び直す.
理論面はこちらで学ぶ. 記法や呼称もこちらに準拠する.
小西・北川『情報量規準 (シリーズ・予測と発見の科学)』
統計学の入門を一通り終えた人には, この本をオススメする. とくに, 3章まで読むことをオススメ. というか他の章はまだ読めていないので, どんなものなのか分からない.
【学習内容】
AICの復習
最初に, カルバック・ライブラー情報量(以下, KL情報量, 相対エントロピーともいう. 符号を逆にするとボルツマンのエントロピーという.)からAIC(赤池情報量規準)の導出への流れを高速で確認する. ところどころ飛躍が激しいところもあるが, 詳細はテキストで確認されたい.
以下, 確率変数\(X\)の真の分布を\(G(X)\), あてはめるモデルの分布を\(F(X)\)とする. 統計学では, データというものは真の分布から生成されてくるものだという大前提を置く(追記(2023/02/24) : この大前提は, 統計学の主義によるそうなので, 必ずしもこの限りではないようだ. このあたりの話がいつか整理できたら書くことにする. 今はこの大前提を置く主義をとるものとして話を進める.). この大前提は頻度論の枠組みだろうとベイズの枠組みだろうと変わらない. そして, その真の分布がどのようなものかは分からないが, 真の分布に近似したものをモデリングという営みによって作成していく. モデルをつくったとして, そのモデルがどれくらい良いものなのかを知りたいというのは自然な気持である. そのモデルがどれくらい真の分布に近いのかを測る尺度の1つがKL情報量である.
確率分布間の近さのようなもの(数学でよく出てくる距離の公理は満たしていない)として定義されるKL情報量は\[I(g;f)=E_G[log(\frac{g(X)}{f(X)})]=E_G[log(g(x))]-E_G[log(f(X))]\]と分解でき, 第1項は真の分布に依存し, モデルに依存しない. 真の分布\(G(X)\)の負のエントロピーともいう. 第2項の\(E_G[log(f(X))]\)は平均対数尤度と呼ばれる. これはモデルの分布\(F(Z)\)で定義される汎化損失の符号が逆になったものでもある.
KL情報量の符号を逆転させたボルツマンのエントロピーは, モデルから得られたサンプルの相対度数が真の分布と一致する確率の対数に比例する量として解釈できる.
また, KL情報量は\(G(X)\)を\(F(X)\)でシミュレートしたとき, \(G(X)\)からの乱数のように見える確率の減少速度として捉えられる話やKL情報量が自然に出てくる話は黒木玄先生のアップしている資料や過去のツイートから学ぶことができる. まだしっかり理解できていないので, ここでは説明しない.
『Kullback-Leibler情報量とSanovの定理』
http://www.math.tohoku.ac.jp/~kuroki/LaTeX/20160616KullbackLeibler.pdf
(ちなみに, KL情報量をテイラー展開すると二次の項にフィッシャー情報量が出てくる.)
したがって, KL情報量を最小化すれば, 真の分布に近い分布がモデルから得られると考えられる. 上に書いたように第1項は真の分布にのみ依存するので, モデルの変化に対して不変である. このため, 第2項の平均対数尤度を最大化すれば(汎化損失を最小化すれば), 良いモデルが得られると考えることができる.
さて, 平均対数尤度を最大化したいわけだが, 平均対数尤度は普通未知である(つまり, 人工データでも作らない限り分からない)真の分布に依存しているので, この量を上手く推定したいというモチベーションが生まれる. 自然に思いつくものとして, 対数尤度を計算してサンプルサイズで割ったものは平均対数尤度の一致推定量となる. (最尤推定法は自然な推定法なのは, この事実があるためでもある.). つまり対数尤度は, \(nE_G[log(f(Z|\hat\theta)]\)の推定量となる.
しかし, この対数尤度による推定量はバイアスを持ち, \[b(G)=E_{G(x_n)}[log(f(X_n|\hat\theta(X_n)))-nE_{G(z)}[log(f(Z|\hat\theta(X_n))]]\] と表され, この大小はパラメータの次元に依存している. このバイアスは3つの項に分解され, 頑張って計算すると漸近的に\[b(G)=tr(I(\theta_0)J(\theta_0)^{-1})\]となる. ただし, \[I(\theta)=\int{g(x)\frac{\partial log(f(x|\theta))}{\partial \theta}\frac{\partial log(f(x|\theta))}{\partial \theta'}}dx\], \[J(\theta)=-\int{g(x)\frac{\partial^2 log(f(x|\theta))}{\partial \theta \partial \theta'}}dx\]である.
想定しているパラメトリックなモデルの族の中に\(g(x)\)が含まれている場合は, \(I(\theta_0)=J(\theta_0)\)となる.
このIとJをデータから推定する一致推定量を\(\hat I\), \(\hat J\)として, バイアスの推定量\(\hat b = tr(\hat I \hat J^{-1})\)が得られる.
モデルの対数尤度をバイアスで補正した量が平均対数尤度を推定したものとなっている. これに-2を乗じたものが情報量規準TIC(竹内情報量規準)と呼ばれる.
上に書いたように, 真のモデルがモデル族に含まれる場合, IとJが等しくなり, バイアスが\(p\), つまりモデルの自由パラメータ数となる. このバイアスで補正して, 平均対数尤度を推定したものを-2倍したものが, 情報量規準AIC(赤池情報量規準)である.
-2倍しているのはカイ二乗統計量との関係でスケールを合わせているだけである. KL情報量に戻ってみると, 第二項の平均対数尤度が大きければよく, つまり平均対数尤度の符号を反転させたものを小さくすればよい. ゆえにAIC, またはTICが小さければ小さいほど, モデルは真の分布に近いことと解釈できる. TICに比べてAICは簡単に計算できて便利だ. ただ, 導出に使われる仮定があってその限界を持っていることや, AICもTICも確率変数なので不確実性な揺れを持っていることを忘れてはならない.
数値実験
AICの復習を終えたところで本題に入る.
『AICをぶっ壊す』というドキュメントでは, 真の分布が確率\((cr1,1-cr1)\)で\(N(m_1,sd_1^2)\)と\(N(m_2,sd_2^2)\)の混合正規分布であるが, 想定するモデルが正規分布となっているケースで, 平均対数尤度と対数尤度のバイアスを数値計算し, それとAICのバイアス部分, つまり正規分布の自由パラメータの個数である2との違いを見ている. ドキュメントでは混合正規分布の他に一様分布と対数正規分布で計算しているが, ここでは混合正規分布のみ扱う.
以下のように関数を書きかえて定義する. 引数はサンプルサイズと真の分布のパラメータ5つをとる.
rewrite <- function(nd, cr1, m1,sd1,m2,sd2){
set.seed(20190824)
# generate data (nd : sample size, cr1 : miture proportions)
# 確率(cr1,1-cr1)でN(m1,sd1^2)とN(m2,sd2^2)から生成する混合正規乱数生成
gdata1 <- function(nd, cr1){
n1 <- sum(runif(nd) < cr1)
c(rnorm(n1,m1,sd1),rnorm(nd-n1,m2,sd2))
#c(dqrng::dqrnorm(n1,m1,sd1), dqrng::dqrnorm(nd-n1,m2,sd2))
#{dqrng}は乱数生成高速化パッケージだが、この程度だと関数呼び出しあたりが悪いのかむしろ遅かった。
}
lf1h <- numeric(3000)
lf1f <- numeric(3000)
for(iii in 1:10){
for(kkk in 1:3000){
xxt <- gdata1(nd, cr1)
mean1 <- mean(xxt)
sdxxt <- sqrt(sum((xxt-mean1)^2)/nd)
lf1h[kkk] <- nd * log(1/(sqrt(2*pi)*sd1)) - 0.5 * nd
lf1 <- numeric(20)
for(kk in 1:20){
xxf <- gdata1(nd, cr1)
lf1[kk] <- nd * log(1/(sqrt(2*pi)*sd1)) - 0.5 * sum((xxf-mean1)^2)/sdxxt^2
}
lf1f[kkk] <- mean(lf1)
}
print(mean(lf1h - lf1f))
}
}混合正規乱数の最適化についてはこのページを参考にした.
http://www.okadajp.org/RWiki/?R%E3%82%B3%E3%83%BC%E3%83%89%E3%81%AE%E6%9C%80%E9%81%A9%E5%8C%96%E4%BE%8B%EF%BC%9A%E6%B7%B7%E5%90%88%E6%AD%A3%E8%A6%8F%E4%B9%B1%E6%95%B0%E3%81%AE%E7%99%BA%E7%94%9F%E3%82%B3%E3%83%BC%E3%83%89
これにより, ドキュメントと同様の結果が出るコードが走る. 混合正規分布の割合パラメータを0, 0.25, 0.5, 0.75, 1のケースで確認する.
rewrite(nd=50, cr1=0.0, m1=2,sd1=3,m2=2,sd2=7)## [1] 2.181348
## [1] 2.236288
## [1] 2.09223
## [1] 2.095036
## [1] 2.510326
## [1] 2.165027
## [1] 2.051087
## [1] 2.084752
## [1] 2.171209
## [1] 2.248626rewrite(nd=50, cr1=0.25, m1=2,sd1=3,m2=2,sd2=7)## [1] 2.530909
## [1] 2.58092
## [1] 2.373811
## [1] 2.473825
## [1] 2.734282
## [1] 2.486
## [1] 2.389523
## [1] 2.442735
## [1] 2.366403
## [1] 2.529645rewrite(nd=50, cr1=0.5, m1=2,sd1=3,m2=2,sd2=7)## [1] 3.094842
## [1] 3.027269
## [1] 2.839038
## [1] 2.985219
## [1] 2.986768
## [1] 2.98658
## [1] 2.938266
## [1] 2.913371
## [1] 2.812462
## [1] 3.133202rewrite(nd=50, cr1=0.75, m1=2,sd1=3,m2=2,sd2=7)## [1] 3.630774
## [1] 3.527391
## [1] 3.205531
## [1] 3.396918
## [1] 3.753438
## [1] 3.520168
## [1] 3.45365
## [1] 3.335673
## [1] 3.162852
## [1] 3.696705rewrite(nd=50, cr1=1.0, m1=2,sd1=3,m2=2,sd2=7)## [1] 2.181348
## [1] 2.236288
## [1] 2.09223
## [1] 2.095036
## [1] 2.510326
## [1] 2.165027
## [1] 2.051087
## [1] 2.084752
## [1] 2.171209
## [1] 2.248626混合正規分布の割合パラメータが0と1のケースは普通の正規分布になっていて, 真の分布が想定しているモデル族に含まれるのでAICのバイアス部分の2で上手くいってそうなことが分かる. 一方でそれ以外のケースではAICのバイアス部分の2よりも大きくなっていて, AICのバイアス部分では過少推定しているであろうことが見て取れる.
せっかく書き換えをしたので速度比較をするため, ドキュメントのコードからコメントを抜いてそのまま書き写した関数も用意した.
allcopy <- function(){
gdata1 <- function(iseed, nd, cr1){
cr1 <- 0.5
xxa <- NULL
set.seed(iseed)
for(ii in 1:nd){
xx1 <- rnorm(n=1, mean=2, sd=3)
xx2 <- rnorm(n=1, mean=2, sd=7)
uni1 <- runif(n=1, min=0, max=1)
xxa[ii] <- xx1
if(uni1 > cr1) xxa[ii] <- xx2
}
return(xxa)
}
nd <- 50
lf1h <- NULL
lf1f <- NULL
for(iii in 1:10){
for(kkk in 1:3000){
iseed <- kkk + iii*912
xxt <- gdata1(iseed, nd, cr1)
mean1 <- mean(xxt)
sd1 <- sqrt(sum((xxt-mean1)^2)/nd)
lf1h[kkk] <- nd * log(1/(sqrt(2*pi)*sd1)) - 0.5 * nd
lf1 <- NULL
for(kk in 1:20){
iseed <- kkk + kk * 471 + iii * 412
xxf <- gdata1(iseed, nd, cr1)
lf1[kk] <- nd * log(1/(sqrt(2*pi)*sd1)) - 0.5 * sum((xxf-mean1)^2)/sd1/sd1
}
lf1f[kkk] <- mean(lf1)
}
print(mean(lf1h - lf1f))
}
}速度比較をしてみると
system.time(allcopy())## [1] 2.865864
## [1] 2.775913
## [1] 2.744208
## [1] 2.960808
## [1] 2.986057
## [1] 3.040051
## [1] 2.98391
## [1] 2.904098
## [1] 2.917919
## [1] 2.878093## user system elapsed
## 218.67 0.01 218.95system.time(rewrite(nd=50, cr1=0.5, m1=2,sd1=3,m2=2,sd2=7))## [1] 3.094842
## [1] 3.027269
## [1] 2.839038
## [1] 2.985219
## [1] 2.986768
## [1] 2.98658
## [1] 2.938266
## [1] 2.913371
## [1] 2.812462
## [1] 3.133202## user system elapsed
## 11.53 0.02 11.54だいぶ速くなったように思える. 並列化や関数呼び出しの省略を上手く使えば, もっと速くなりそうな気もするが, コードの分かりやすさ(と面倒くささ)から今回はこれでとどめておく.
さて, 速度は今回の話の本質ではない. AICがいい感じに平均対数尤度の-2倍を推定できていないことが数値実験から分かったというのが今回の本質である.
AICの定義するバイアス部分では平均対数尤度を対数尤度で推定したときのバイアスの補正が上手くいかない例が今回の数値実験例であった. 導出のところで確認すべきところだったのだが, 上手くいかないケースは, 小標本の場合, p/nが大きい場合, 真の分布が含まれていない場合, 大標本の場合とがある. 加えて, 非正則の場合もある.
庄野『モデル選択手法の水産資源解析への応用－情報量規準とステップワイズ検定の取り扱い－』
https://www.jstage.jst.go.jp/article/jjb/27/1/27_1_55/_pdf
は, 小標本の場合, 大標本の場合,真の分布が含まれていない場合の数値実験をして, 各情報量基準でモデル選択をして正解しているか否かを見ている. これを見てみると, AICはそこそこ頑張っているように見える.
バイアスの推定が上手くいっていない場合にどの程度上手くいかないかが分からないが, 多少仮定から外れていてもモデル選択は上手く行えるのかもしれない. ただ, 渡辺澄夫先生曰く, 複雑なモデルつまり,構造をもつモデル(混合正規分布、神経回路網、ベイズネットワーク、隠れマルコフモデル、 縮小ランク回帰、深層学習などの隠れた変数・階層構造・モジュール構造を持つモデル)ではAICなどではモデル評価が上手くできないらしい. 渡辺澄夫先生のHP『広く使える情報量規準(WAIC)』
http://watanabe-www.math.dis.titech.ac.jp/users/swatanab/waic2011.html
今回の数値実験混合正規分布もそれにあたり, このような場合はWAICを使えば良いらしい. 今回の実験は, 平均が同じ2つの正規分布の混合分布であったが, これダメなタイプなのかいまいち分かっていない. 今度, 単峰と複峰で乱数を発生させてAICとWAICの挙動を実験してみたい. AICがダメなパターンは尤度関数または事後分布が正規分布で近似できないケースなので、尤度関数の形状を眺めてみたい. いつかそんな記事を書く.
[追記：2019.10.12 書けたのでリンク貼る.]
ドキュメントのタイトルは『AICをぶっ壊す』という過激なものであったが, ぶっ壊すというよりはAICの数理的な構造を理解して, その妥当性と限界を知ることが重要なように思われる. これらを理解した先にWAICを理解して使えるようにしていけばいいような気がしてならない.
真の分布が正則で学習モデルに含まれいるなら, WAICはAICに, 正則で学習モデルに含まれていないならTICに, 確率変数として等価だそうだ。
http://watanabe-www.math.dis.titech.ac.jp/users/swatanab/cross-val.html
【学習予定】
WAICは存在は知っているがまだ全然理解できていない. 数値実験したり, 数式を読み解きながら, ゆっくりでもいいから理解していきたい.
渡辺『ベイズ統計の理論と方法』
https://www.amazon.co.jp/dp/4339024627/ref=cm_sw_r_tw_dp_U_x_dPuzDbACQ08KB
