多項式フィッティングにおける次数選定：赤池情報量規準についてまとめてみる

赤池情報量規準についてまとめます．

会社の同期と静岡旅行中に赤池情報量規準についてあれこれ話してたので，その辺をまとめた記事です．
私はもともと機械系出身なので情報系の知識は皆無なんですが，
そのぶん初心者向けになったかなーと思います．

※そのせいで結構アバウトな部分があるとは思います．ご了承ください．

これはどんな記事

「多項式フィッティングの次元どうやって決めよう」
「これを最小化すればいいよ　つ（AIC）」
「なにこれ？」

って感じの問題に答えます．

流れは，情報エントロピー，KLダイバージェンス，最尤推定，AICの説明と，
最後にpythonで実際に多項式フィッティングの例題でAICを計算して結果を確認します．

赤池情報量規準とは

多項式フィッティングを行う際などに，どのように次元の数を決めれば良いか，という問題が生じます．次元が少なすぎると関数の表現力が小さすぎて当てはめが上手く行きませんし，多すぎるとオーバーフィッティングして困ります．その問いに答えてくれる一つの方法が赤池情報量規準（AIC: Akaike's Information Criterion）です．

ガウスノイズを含んだあるデータに対して，$M$個のパラメータを持つ関数でフィッティングを行うとします．このとき，AICは以下の式で示されます．

{\rm AIC}=-2\ln L + 2M

ここで$L$は最尤推定下での尤度です．この値が最小になるように$M$の値を決めるのが良いとされています．

確かに，この式を見ると，尤度はなるべく大きく，かつパラメータ数は少なく，という感じにバランスを取っていることは分かりますが，なぜこの形なんでしょう．と，疑問に思ったので調べてみました．

最尤推定って？

ほとんどの人はご存知だと思いますが，一応．

上記の関数フィッティングの問題は，パラメータ$\theta$を持つ関数$y=F(x;\theta)$を用いて，手持ちのデータ関係$D:x\rightarrow y$を近似することでした．最尤推定というのは，このパラメータ$\theta$の決め方の一つの方法です．

まず，ここで言う尤度とは，関数$F(x;\theta)$が存在した時に，データ関係$D:x\rightarrow y$が観測される確率のことを言います．

こちらに分かりやすい文章が載っていたので，お借りします．

表が出る確率が$θ$であるコインを100回投げて70回表が出る確率は，反復試行の確率の公式より

L(\theta)=_{100}C_{70} \theta^{70}(1−\theta)^{30}

これが尤度関数である。これを最大にする$θ$がもっともらしい$θ$である。

このときの最尤パラメータは，対数尤度を取って計算すると$\theta=0.7$が尤度を最大にすることが分かります．

関数フィッティングの場合は，観測ノイズが入っていると仮定しているので，関数$F(x;\theta)$による$y$と$x$の近似関係は，観測ノイズ$w$を用いて

y=F(x;\theta)+w

と表されることになります．$w$が確率的なので，上式に従って$x$を与えた時の$y$の値は確率的に算出されます．例えば，$x_1$を与えたときに$y_1$が得られる確率は$0.5$である，とかです．この確率を全データ$x_i,y_i\in D$についてすべて積算したときの合計が，そのデータ$D$が得られる確率，つまり尤度となります．

式で書くとこんな感じでしょうか．

L(\theta) = \Pi_i p(y_i | x_i; \theta)

ここで$p(y_i|x_i;\theta)$は，パラメータ$\theta$の下で，$x_i$を与えた時に$y_i$が得られる確率です．この尤度が最大になるように$\theta$を定めるのが，最尤推定というわけです．

ちなみに最小二乗法による推定は，観測ノイズが平均$0$のガウスノイズである場合の最尤推定値になることが知られています．

AICの導出の概要

最尤推定は「**与えられた関数$F(x;\theta)$**に対して，尤度が最大になるようにパラメータ$\theta$を決める手法」でした．対して，AICは「どのようにして関数$F(x;\theta)$を決めれば良いのか」に答える方法です．

結果からまとめると，真のモデルの確率分布，と近似モデルの確率分布との差が最小になるようにするためにはAICを最小にすればよい，ということが導けるそうです．つまり，AICが最小のとき，真の確率分布$g(x)$に対して近似モデルの確率分布$f(x;\theta)$が$g(x)$に最も近くなります．

これの説明のためには，確率分布同士の差をどうやって測るか，という問題について説明する必要があります．

※混乱しそうなので，表記をまとめておきます
・$\theta$：関数フィッティングにおけるパラメータ
・$M$：フィッティングパラメータの数．なので$\theta \in R^{M}$．
・$D$：入出力関係$x$と$y$の（ノイズ付き）観測データ
・$F(x;\theta)$：パラメータ$\theta$を持ったデータ$D$の近似モデル関数
・$f(x,\theta)$：近似モデル関数$F$の下での確率分布
・$g(x)$：真の確率分布

カルバック・ライブラー ダイバージェンス

AICではカルバック・ライブラー(KL) ダイバージェンスという値を用いて確率分布同士の差を測っています．確率分布$g(x)$から見た$f(x)$のKLダイバージェンスは以下の式で書けます．

{\rm KL}(g,f) = \int_\Omega g(x) \log \left( \frac{g(x)}{f(x)} \right)dx

$\Omega$は対象とする空間を表しています（とりあえず全空間だと思えばいいでしょう）．この式は$f(x)=g(x)$のときに$KL=0$で，それ以外の時には$KL> 0$となることが分かっているそうです（ギブスの不等式というらしい）．なのでまあ，確率分布同士の近さを図る尺度としては分からなくもないですが，なんでこの形なのか釈然としないですよね．別に近さを測る尺度ならもっと簡単なものでいいじゃん，と．

とりあえずwikiによると

今確率変数X が本当は分布P に従っているのに、誤って分布Q に従っていると判断してしまった場合、 本来の最小値よりも多くのビット数を必要としてしまう。 カルバック・ライブラー情報量は、このような誤りを犯してしまった場合に 余分にかかってしまうビット数の平均値を表す。

らしいですが，私には初見じゃ全く意味が分からなかったです．

いろいろと理由があるそうなのですが，情報量による説明が一番しっくり来ました．
まず，なぜ確率分布の近さをKLダイバージェンスで測るのか，という疑問に答えるためには，情報理論におけるエントロピーを理解する必要があります．

情報におけるエントロピーとは

機械科出身の私にとってエントロピーと言われると熱力学しか出てこないのですが，情報理論においても同様のニュアンス？を持ったエントロピーという概念があるそうです．簡単に言うと，ある事象がどれくらいの情報を含んでいるかを示す指標らしいです．以下の例でそれを説明します．

＜情報量の例題＞
ある乱数生成器が0~1023の整数を生成し，その値を"0","1"で符号化して信号として相手に送るケースを考えます．このとき純粋に0~1023の値をビットで対応させようとしたら，10個のビットが必要になるわけです（2^10=1024）．なので，1つの値を相手に送るために必要な平均ビット数は10になります．

ここで，0~1023の整数値の中で，0が圧倒的に生じる可能性が高い場合はどうすればよいでしょうか．このときは頻発する0に対して，10ビットかけて"0000000000"の信号を送るのは無駄なように思えます．少ないビットで効率的に情報を送るために，"0"が送られてきたときは整数0を表し，"1"が送られてきたときはその後の10ビット分が1~1023の整数を表すように取り決めをしたとしましょう．1~1023の情報を送るためには先頭の"1"を含めて11ビットが必要ですが，もしかしたらこっちの方が効率がいいかもしれません．
ここで，整数0が90%の確率，それ以外の整数が残りの等確率で生じるとしたときの，情報を送るための平均ビット数を計算してみます．$n(i)$を整数$i$を送信するために必要なビット数，$p(i)$がその整数が生じる確率とすれば

\sum_{i=0\sim 1023}n(i)p(i)=1\times1\times \frac{90}{100}+1023\times11\times \frac{10}{100 \times 1023}=2.0

となります．なので，平均2ビットの情報量でデータが送信できてしまうわけです．これは全ての整数を10ビットで送るよりも遥かに効率がいいことが分かります．
＜例題終わり＞

このことから，確率的に観測されるあるデータを情報に変換する場合には，情報量が最小となるような最適な変換方法があると予測できます．この変換における情報量の下限がエントロピー$H$と呼ばれ，以下で表されます．

H=-\sum_i p(x_i)\log p(x_i)

この式は，確率$p(x_i)$で生じる事象には$-\log p(x_i)$の符号長を与えればよい，という式だそうです．上で計算した整数0が90%の確率で生じるケースのエントロピーを計算してみると，$H \simeq 1.47$となるので，$2.0$は結構いい線を行ってたわけですね．

このエントロピーの定義をそのまま連続系に変換すると

H=-\int_\Omega p(x)\log p(x)dx

になります．こいつがKLダイバージェンスに現れるやつですね．

なお，ここではビットを用いて説明したので$log_2$を採用しましたが，別に情報はビットで表す必要はないので，対数に何を取るかについては一般化できるみたいです．

KLダイバージェンスは，なぜあの形なのか

KLダイバージェンスの章で書いた式を変形すると

KL(g,f)=\int_\Omega g(x)\ln g(x)dx - \int_\Omega g(x) \ln f(x)dx

となります．関数フィッティングの話を思い出すと，$g(x)$はデータ$D$を生み出した真の確率分布，$f(x)$は$g(x)$を近似しようとして作った近似モデルの確率分布でした．ここで，以下のエントロピーを定義します．

H(g,g) = -\int_\Omega g(x)\ln g(x)dx,\\
H(g,f) = -\int_\Omega g(x)\ln f(x)dx.

情報量の概念に沿って考えると，$H(g,g)$は確率分布$g$から生成されるデータを最適にエンコードしたときの情報量，$H(g,f)$は確率分布$g$から生成されるデータを$f$から生成されたと近似してエンコードしたときの情報量，となります．なので，$H(g,g)<H(g,f)$ですね．

この定義をKLダイバージェンスを使って書き直すと

KL(g,f)=-H(g,g) +H(g,f)

こうなります．この式を見ると，これは**「$g(x)$から生成されたデータに対して，$g(x)$を使って最適にエンコードしたときのエントロピー（情報量）と，$g(x)$の近似モデル$f(x)$を使ってエンコードしたときのエントロピー（情報量）の差」**を表していることが分かります．つまり，近似によって抜け落ちてしまった情報量ってことですね．

これで何となくKLダイバージェンスのイメージが掴めたと思います．AICに戻りましょう．

AICはKLダイバージェンスの第2項の近似

KLダイバージェンスで確率分布の近さを図れることが分かりました．それならば，真の分布と近似分布のKLダイバージェンスを最小にしようとするのは必然でしょう．AICがやっているのは，その最小化です．

KL(g,f)=-H(g,g) +H(g,f)

先ほど説明したこの式を最小化したいのですが，第1項は真の分布$g$しか出てこないので値は確定しています．なので問題は，どうやって$f$を選べば第2項の$H(g,f)$が最小化されるのか，になります．

KLダイバージェンスの第2項は，引数とパラメータもちゃんと書くと

H(g(x),f(x|\theta)) =  -\int_\Omega g(x)\ln f(x|\theta)dx

でした．この値を最小化するために式を解析的に書き下したいのですが，そのときにいくつかの条件（仮定？）が入ります．それが以下の２つ．
１．サンプル数は無限大とする
２．モデルのパラメータ$\theta$は最尤推定によって決める．このとき，尤度を$\theta$でテーラー2次近似したときの極小点として$\theta$を求める．

１に関して，$H(g,f)$には未知の確率分布$g$が入ってきているので，直接計算はできません．AICではこの$g$を，サンプル数が無限大の経験分布で近似するということを行っているそうです．経験分布というのはデルタ関数からなる確率分布で，n個のサンプルデータのすべてをn個のデルタ関数で表したやつだと思えばよさそうです．

２に関して，$H(g,f)$の計算のためには，まず$f(x|\theta)$のパラメータである$\theta$を決めないといけません．AICではこれを，最尤推定下での$\theta$を使って$H(g,f)$を計算するとしているそうです．この最尤推定は一般的に解くことはできないのですが，ここでは尤度を$\theta$で2次までテーラー展開した関数に対して，最尤推定を行っています．式で書くと

\frac{\partial H(g,f(\theta))}{\partial \theta}=0

の式の左辺をテーラー近似して，$\theta$について解いているみたいです．この解を最尤推定の解$\theta^*$としておきましょう．

これらの近似と仮定を使うと，KLダイバージェンスの第2項は

-\log f(x|\theta^*)+M

と書けるそうです．$f(x|\theta^*)$は最尤推定下での近似モデルの尤度です．AICがやりたかったことは，KLダイバージェンスの最小化でしたが，結果的にそれはこの式を最小化すればよいということが分かります．

歴史的経緯で，この式に2を掛けて

{\rm AIC}=-2\log f(x|\theta^*)+2M

としているそうです．

それにしても，ここで$\theta$のパラメータ数である$M$が急に現れるのは違和感ですね...　数式的には尤度関数をテーラー2次近似をした時の2次の項

(\theta - \theta^*)^TS(\theta-\theta^*)

の（$S$はHessian行列）$S$のトレース$\rm{tr}(S)$から$M$が現れてくるそうです．にしても，ここまできれいな式になるもんなんですねぇ．うーん，感動．

いや，でもトレースから次数を出してるってことは，2次近似したらヘッシアンが出てくるって事実から，AICの$M$の項が生じてるってことですよね．それはいいのか？近似した時の高次項が影響してるってことですからね．でも別に二次近似は最尤推定値を求めるために使ってるだけだからいいのかな...
この辺はちゃんと式を追わないとイメージ掴めないですね．

実際に多項式フィッティングでAICを試してみる

AICが本当に役に立つのか，pythonで確かめてみます．

オリジナルデータを次の3次多項式

y = 4 -5x -10x^2 + x^3 + w

から生成します．$w$は正規分布からなる乱数です．こいつから得られたデータに対して，多項式フィッティングを行って多項式係数を推定します．
このとき基本的にオリジナルデータが何次の多項式から生成されているかは分からないので，とりあえず適当に次数を決めてフィッティングします．次のコードでは次数を１～７に対してフィッティングを行い，対数尤度，AICをそれぞれ計算しています．
実際のノイズは分散10としてガウスノイズを与えたんですが，推定時はノイズは未知として，観測データと最小二乗の推定結果から観測ノイズの分散を計算して，正規分布を仮定して尤度を計算してます．

aic_polyfit.py

import numpy as np
import math
from matplotlib import pyplot
from numpy.linalg import pinv


# 最小二乗法で多項式フィッティング

c = np.array([4, -5, -10, 1])
x = np.linspace(0,10,100)
y_measured = c[0] + c[1] * x + c[2] * x**2 + c[3] * x**3 + np.random.normal(0,10,x.size)
pyplot.plot(x,y_measured,label='measured')

for dim in range(1,8):    
    x_mat = x * 0 + 1
    for i in range(1, dim+1):
        x_mat = np.c_[x_mat, x**i] 
    c_est =  np.linalg.pinv(x_mat).dot(y_measured)
    
    y_est = x * 0
    for i in range(0, dim+1):
        y_est = y_est + c_est[i] * x**i
        
    pyplot.plot(x,y_est,label='fit(dim = ' + str(dim) + ')')
   
    # 分散の計算
    sigma = np.sum((y_est - y_measured)**2) / x.size

    # 対数尤度の計算
    p = lambda x,mu,sigma: (math.exp(-(x-mu)**2 / (2*sigma)) / math.sqrt(2*math.pi*sigma))    
    L = 0.0
    for i in range(0, x.size-1):
        L = L + np.log(p(y_measured[i], y_est[i], sigma))
        
        
    AIC = L - (dim+1)
    print('dim = ' + str(dim) +', L = ' + str(L) + ', AIC = ' + str(AIC))

pyplot.legend()
pyplot.grid()

オリジナルデータとフィッティングの結果はこんな感じ．

さすがに１次と２次じゃうまく近似できていません．
各次数に対して対数尤度とAICを計算したのが下の表です．

対数尤度が最大になっているのは次数６の時ですが，AICが最小になっているのは次数４のときですね．
まぁまぁな結果って感じでしょうか．

ほんとはビシッと次数３を示してほしかったんですが，多項式近似だと高次と低次へのノイズの影響とか変わりそうですし，そのまま適用はできないのかもしれないですね．

計算してみて思ったのですが，
サンプル数が増えていくと対数尤度はそれに従って値が減少していきます．
でも次元の値はサンプル数に対して不変（あたりまえ）ですから，
AICの$-2\ln L + 2M$って$n\rightarrow \infty$のときはほとんど$+2M$の項の影響なんて消えちゃうのでは？

なにか間違えてるのかなぁ．

このあたりの知見がある方，コメントいただけると嬉しいです．

※追記
今回のプログラムではサンプル数が$n=100$としましたが，ここを$n=10000$に変えて計算してみると，ほぼ100%の確率でAIC的に次数３がベストという結果が返ってきました．サンプル数$100$は少なすぎたんですね．

AICのまとめと使用上の注意

まとめです！

• AICは真の確率分布と，最尤推定されたパラメータ$\theta^*$を持つ近似確率分布のKLダイバージェンスを最小化するためのもの
• 具体的にはKLダイバージェンスの第2項の近似（第1項はモデルに依存しない）
• 近似にはテーラー2次近似を使ってる
• サンプル数は無限大と仮定してる

あとこれも

• KLダイバージェンスは，ある確率分布に沿って生成されるデータを真の分布と近似分布でエンコードした時のエントロピーの差
• エントロピーとは情報を最適にエンコードしたときに必要となるビット数（の一般化）

このことから，AICを使うときは，最尤推定されたパラメータを使わないと意味がない，尤度関数が複雑すぎるとテーラー近似から離れすぎてダメ，といった使用上の注意点が浮かび上がります．

この問題を解決するためのAICの派生形もたくさんあるらしいです．いつか勉強します！

その他

参考にしたページ
エントロピーについて：https://www.eidos.ic.i.u-tokyo.ac.jp/~tau/lecture/komaba_joho/gen2/slides/3-entropy.pdf
AICについて１：http://www.is.titech.ac.jp/~shimo/pub/koza2000/ohp-20000703b.pdf
AICについて２：http://www4.ncsu.edu/~shu3/Presentation/AIC.pdf
AICについて３：http://lef-t.blogspot.jp/2014/01/aicabc.html
AICについて４：http://takashiyoshino.random-walk.org/memo/keikaku2/node5.html