AICとかBICとかを自分で計算する人はあまりいないと思いますが、訳あってやってみたので記事にしました。たぶん、今回も知っている人は知る必要のない内容です。

対数尤度

モデルを評価する場合、よく使われるのが対数尤度です。なんか、確率のような気がしますが、合計（？）が常に１になるとは限らないので、確率ではないそうです。ベイズなんとか、とかでよく出てきますが、そういう計算をするときは、1になるよう正規化して、掛け算したりなんだりという計算をするようです（分布の更新とか、そういう計算に使う）。と、色々難しそうなのですが、ここではモデルの評価について考えているので、とりあえず、「当てはまりがよさそうな時に大きい値が出るやつ」の対数をとったもの、ということだけ考えて、先に進みます。

さて、ガウス型モデル（誤差がガウス分布に従うもの）の対数尤度は、実は残差平方和や残差平方和平均があれば、簡単に計算できるそうです。

ちゃんとした式としては、こうなるのですが、

見ての通り、第二項は残差の平方和になっていますので、シグマをそれで近似できるとして下記のように変形します。

分散をRSS/nとやって（サンプルの分散で）置き換えているようなのですが、これn-1じゃなくていいのかな？なんか勘違いしているかも。まー、nが大きい時は似たようなものなので、いいんでないかな！

# gaussian likelihood
def GL(rss, n):
    return -n*0.5*(np.log(2.0*np.pi*(rss/n))+1)

赤池情報基準

あの、有名な、あれですな。あの、情報量に基づいたモデル選択基準というやつです。前にGoogleの検索画面に、提唱者の赤池氏が出てましたね。というわけで、有名なやつです。

ln(L)の部分が先ほどの対数尤度で、kはモデルの自由パラメータです。要するに、あれこれいじるあのパラメータですな。カーネルの幅とか、なんとかの係数とかです。回帰分析なら、回帰係数と切片がパラメータみたいな感じ。

def AIC(rss,n,p):
    return n*np.log(rss/n) + 2*p

私はかなりのボンクラなので、初めてこの式を見たとき、nが1000とかのときに、こんなの計算して意味あるのかなぁ？とつい思ってしまいました。浅はかですね。

実はAICというのは、相対的な基準なので値の差を見てモデルの評価基準とするものです。たとえば、モデルAとモデルBの尤度が同じく-1998で、Aがパラメータ二個、Bが一個なら、AICはそれぞれ2002と2000となるので、Bの方が（AIC的な基準で）良いモデルです。パラメータ一個分差が出るので差は2です。こう考えると、データがどんだけ増えようが、AICの「差」には影響なさそうです。もちろん、実際のモデルでは尤度の方も変化するため、そこがトレードオフとなっていい感じの選択基準となってくれます。よって、データが同じじゃないと、評価基準になりません。要するに、「2000と2001なんだから、たいして差はない」みたいな話ではないということです（よく知らん頃はそう思ってた）。

というわけで、上記の「AICを計算するかんすう」も、出てくる値は結構適当。ただ、一応（AIC的な基準において）良いモデルであれば小さい値、悪いモデルなら大きい値になるようになっているはずです。ちなみに、値の大小が尤度と逆なので注意（AICの第一項にマイナスがついてるため）。これでいいの～？と思う人は、英語版のWikipediaを見たら、ちゃんと大丈夫そうだという計算例が載っていますので、そちらを見てみてください。

ベイズ情報基準

計算式は下記のとおりです。AICに似ていますが、第二項が尤度で、第一項がペナルティですね、この式だと。log(8) ≒ 2.079...程度なので、見ての通りn>7のときには、AICより厳しい感じになるのが特徴だそうです。

ベイズとついているので、ベイジアンな指標なのかと思い込んでいたところ、どうやら、導出の際にベイズ統計的な議論というか、仮定というかがなされているだけであり、この値そのものは別にベイズ的な（？）方法で計算するわけではないようです。そのため、計算は下記の通りかんたん。

def BIC(rss,n,p):
    return n*np.log(rss/n) + np.log(n)*p

やれ、事前分布だの事後分布だのということはしなくても、BICの計算はできる模様です。これも、AICと同じで、単純化した計算式です。まぁ、最初のところで尤度の計算を（たぶん）ちゃんとしているので、それを使えばいいような気もするのですが、どういう式が一般的なのかいまいちわからなかったので、こちらもWikipediaに倣っています。

計算例

さて、理論的な議論は偉い統計学者がすでに済ませているのでいいとして、我々一般人からすると、どんな感じで差が出るのか気になるところです。そんなわけで、実際に計算してみました。

想定としては、データ数がNのとき、モデルの分散が1/Nになるという状況を考えます（ガウス分布なので、そんなもん？）。また、分散はパラメータ数kに反比例します（希望的観測）そういった状況で、パラメータ数を1..10、データ数Nを100...1000まで変化させてAIC/BICをグラフにプロットしてみました。

BICです。

AICです。

パラメータを増やせば増やすほど、それに反比例（パラメータを二個にすると分散が半分という夢のような状況）して誤差が小さくなるので、AIC的にはいくらパラメータを増やしても文句を言わない感じになっています。見ての通り、データ数が多い所ほど、パラメータ増やせば増やすほどAICは良くなっています。

コストと精度

さて、もうちょっと真面目にやってみます。AICの計算式は下記の通りですので、パラメータを一個増やすとAICは2増えます。

def AIC(rss,n,p):
    return n*np.log(rss/n) + 2*p

逆に言えば、第一項が、2程度減少しなければ、パラメータは増したことによりAICは全体として増加するということです（この2って何なんでしょう？ちゃんとした根拠があるのだとは思うのですが…）。

というわけで、sympyで計算してみると（だって、面倒なんだもの。こんなの慣れてる人なら暗算するレベルだよね）、パラメータをひとつ増やしたとき、rssがexp(-2/n)程度減少すれば、帳尻が合うことが分かります。データ数nが1000の場合、これは数値で言うと0.998程度です。n=100のときは約0.98、n=10の時は約0.819です。というわけで、データ数が1000の場合は、パラメータ一つにつき0.2%程度、100の時は2%程度、10の時に至っては18%程度の向上が見込めなければ、AICは（パラメータを増やしたモデルの方が）優れているよ、とは言ってくれません。というわけで、AICはデータ数がもたらす影響も、こんな感じで織り込んでいるようです。

おわり

今回はガウス型を仮定しているため、データを増やせば誤差のガウス分布の裾が指数的な感じで減少していくことが想定されます。なので、n=100のときに精度が1%向上するのと、n=1000の時に制度が1%向上するのとは雲泥の差です（n=100のときは、ラッキーパンチで1%程度向上したのかもしれないですが、n=1000ならラッキーでそんなことが起きる可能性は低いため）。それを図にしてみてみようと思ったのですが、面倒臭くなってしまったので、これで終わりです。もし、AICやBIC以外の評価指標についても勉強したら、追記するかもです。

以上です。