本記事の概要

みなさんご存じAICのことを改めて整理しよう
「パラメータ数×2を引けば（足せば）いいんでしょ」以上のことを言えるようになろう

はじめに

最近統計の本を読んでいて、AIC（赤池情報量規準）が出てきました。もちろん知ってはいますし、大学院の講義で解説もいただいたのですが、残念なことに
「モデルの選択に使う」
「過学習をしづらくなる」
「パラメータ数×2を引く（足す）」
くらいのことしか言えない自分に気づきました。

そこで、本記事では AICについて改めて整理して、もうちょっとかっちょいい説明ができるようになることを目的とします。

本記事の読み方

最尤法、KLダイバージェンス、モデル選択などの基本的な内容から触れています。
これらについてある程度勉強したことのある方は、(確認)とついている章は読み飛ばしていただいても問題ないと思います。

(確認) 最尤推定量

最尤法の確認

まずは最尤法について確認しておきましょう。
データ $y_{i} (i=1,\cdots,N)$ があるパラメトリックモデル $f(y|\theta)$ から得られたと仮定します。
このとき，$y_{i}$の同時分布

f(\boldsymbol{y}|\theta)=f(y_{1}, y_{2},\cdots,y_{N}|\theta)

が尤度です。尤度はパラメータ$\theta$の関数とみることができて、$\theta$を固定したときに手元のデータがモデル $f(y|\theta)$ から得られる確率を表します。

ここで、データをもとにパラメータ$\theta$を定めることを考えます。尤度が大きければ手元にあるデータが得られる確率が大きいわけですから、なるべく尤度が大きくなるように$\theta$を定めればデータにフィットしたモデルが得られそうです。この観点でパラメータを定めるのが最尤法です。

データ $y_{i} (i=1,\cdots,N)$ が独立に得られたと仮定したとき、最尤法による$\theta$の推定値は以下のように定式化されます。

\hat{\theta}  = \underset{\theta} {\operatorname{argmax}} \log{f(\boldsymbol{y}|\theta)}= \underset{\theta} {\operatorname{argmax}} \sum_{i=1}^N \log{f(y_{i}|\theta)}

ここで、$\log{f(\boldsymbol{y}|\theta)}$は尤度の対数をとった値で、対数尤度といいます。対数をとる理由の一つは、同時分布 $f(y|\theta)$ を和の形に変形でき微分の計算で都合がよいからです。

このようにして得られた推定値$\hat{\theta}$が最尤推定量です。

最尤推定量の例

最尤推定量の例を挙げます。
データ $y_{i} (i=1,\cdots,N)$ が 正規分布 $N(\mu, \sigma^2)$ から得られたと仮定します。
このとき，データの対数尤度は

\sum_{i=1}^N \left( \log{\frac{1}{\sqrt{2\pi\sigma^2}}} - \frac{(x_i-\mu)^2}{2\sigma^2} \right)

であり、これを$\mu$, $\sigma^2$で偏微分することによってそれぞれの最尤推定量を求めると以下のようになります。

\hat{\mu}=\frac{1}{N}\sum_{i=1}^N y_{i},\quad \hat{\sigma^2}=\frac{1}{N}\sum_{i=1}^N (y_{i}-\hat{\mu})

よく注意点として言われますが、分散$\sigma^2$の最尤推定量はいわゆる「不偏標本分散」ではなく「標本分散」です。

(確認) KLダイバージェンス

AICを議論する準備として欠かせない、カルバック・ライブラー情報量（KL情報量、KLダイバージェンス）についても確認しておきます。

KL情報量の定義

KL情報量は、2つの分布の距離（のようなもの）を数量的に表すものです。分布$f(y)$の$g(y)$に対するKL情報量は次のように定義されます。

I(g;f)=\int g(y) \log{\frac{g(y)}{f(y)}}dy

かみ砕くと、分布$f(y)$が$g(y)$にどれくらい近いかを表しています。$I(g;f)$は、$0$以上の値をとること、$f(y)$が$g(y)$と完全に一致するときにのみ$I(g;f)=0$となることから、2つの分布の近さを測る値として使用されます。

$I(g;f)$は$g$と$f$を入れ替えた$I(f;g)$とは一致しないので、距離関数の定義は満たさず厳密にいうと距離とはいえません。

(確認) モデル選択

モデル選択の目標

AICのモチベーションとなる、モデル選択についてです。
データ $y_{i} (i=1,\cdots,N)$ について、以下の2つのモデルを考えます。

$g(y)$：真のモデル
$f(y|\theta)$：統計モデル（我々が仮定したモデル）

真のモデル$g(y)$がどのようなものかは神のみぞ知るものですが、データからより良い統計モデル$f(y|\theta)$を求めることが我々の目標です。
この「良いモデル」の意味は、今手元にあるデータをよく予測（説明）するモデルではなく、「今後得られるデータを良く予測（説明）してくれるであろうモデル」です。未知のデータに対して予測する力が高いモデルをどのように選ぶのかがモデル選択の本質的な目標になります。

過学習

統計モデルはパラメータ$\theta$によって決まるので、よりデータにあった統計モデルになるような推定パラメータ$\hat{\theta}$ を決める方法の一つが、最尤法でした。

一方で統計モデルはパラメータの値だけでなく、モデルの形そのものでも決まります。
例えば下の図は、10個のデータに対してそれぞれ1次，2次，4次，7次の多項式回帰を行ったものです。

各統計モデルのパラメータは最尤法により推定されたものですが、実際にはどの統計モデルがよいかを決めるのがモデル選択です。
単に予測値（緑線）とデータの値との差が小さいものを良いモデルと考えた場合、選択されるのは右下の7次関数になります。しかし、このモデルは手元のデータに過度にフィットしている、いわゆる過学習とよばれる状態と考えられ、「今後得られるデータを良く予測する」モデルであるといえるかどうかは大きく疑問です。実際、モデルのパラメータを増やして表現力を高くすればするほど手元のデータとの予測差は小さくできてしまいます。

モデル評価とKL情報量

統計モデルの未知のデータに対しての振る舞いをどう評価するかという観点で有用なのが、前出のKL情報量です。
真のモデル $g(y)$ と統計モデル $f(y|\theta)$ のKL情報量

I(g;f)=\int g(y) \log{\frac{g(y)}{f(y|\theta)}}dy

は、$g(y)$ と統計モデル $f(y|\theta)$ の近さを数量化します。
つまり、統計モデルの中で $I(g;f)$ を最も小さくするであろうモデルを選択することは、モデル選択の視点でよい選択方法といえそうです。
このKL情報量を用いたモデル選択の延長上にAICがあります。

AICの成り立ち

ようやくここからAICが登場します。KL情報量からAICにつながる内容です。

平均対数尤度

簡略化のため、2つの統計モデル$f_{1}(y|\hat{\theta_{1}})$、$f_{2}(y|\hat{\theta_{2}})$のどちらかを選択する状況を考えましょう。

$\hat{\theta_{1}}、\hat{\theta_{2}}$はどちらも、データ$y_{i}(i=1,\cdots,N)$から求めたモデルパラメータの最尤推定値とします。

選択の基準は真のモデル $g(y)$ とのKL情報量をより小さくできるのはどちらであるかです。

まず、KL情報量$I_{m}(g;f_{m}) (m=1,2)$は次のように変形されます。

\begin{align}
I_{1}(g;f_{1})=\int g(y) \log{g(y)}dy-\int g(y) \log{f_{1}(y|\hat{\theta_{1}})}dy \\
I_{2}(g;f_{2})=\int g(y) \log{g(y)}dy-\int g(y) \log{f_{2}(y|\hat{\theta_{2}})}dy
\end{align}

ここで、$g(y)$を知ることはできないため、$I_{m}(g;f_{m})$を厳密に計算することはできません。

ただし、$I_{1}(g;f_{1})$と$I_{2}(g;f_{2})$を比較するのであれば、まだ対処できることがあります。

第1項
- $\int g(y) \log{g(y)}dy$ はどちらの式にも共通している、統計モデルに因らない項です。よって、$f_{1}$と$f_{2}$を比較するうえではこの項は無視できます。
第2項
- 負の符号がついているので、$I_{m}(g;f_{m})$ を小さくしたいなら $\int g(y) \log{f_{m}(y|\hat{\theta_{m}})}dy$ が大きくなれば良いことがわかります。
- 第2項は、その式の形から、対数尤度 $\log{f_{m}(y|\hat{\theta_{m}})}$ の真の分布 $g(y)$ でとった平均（期待値）と捉えることができます。

\int g(y) \log{f_{m}(y|\hat{\theta_{m}})}dy =E_{g}[\log{f_{m}(y|\hat{\theta_{m}})}]

以上のことから、KL情報量によるモデル選択は以下のような問題として定式化できます。

統計モデルの平均対数尤度 $E_{g}[\log{f_{m}(y|\hat{\theta_{m}})}]$を求める。
推定値が最も大きいモデル$f_{m}$を選択する。

平均対数尤度の推定

平均対数尤度$E_{g}[\log{f_{m}(y|\hat{\theta_{m}})}]$を求めて、それをもとにモデル選択をしたいということですが、ここで問題が一つあります。

$E_{g}[\log{f_{m}(y|\hat{\theta_{m}})}]$は真の分布$g(y)$での期待値ですが、どのようにして求めればよいのでしょうか？

その通り、$g(y)$を知ることができないから悩んでいたわけで、期待値を計算したくてもそれは不可能です。

よって、何とか $E_{g}[\log{f_{m}(y|\hat{\theta_{m}})}]$ を推定したい、そのために良い推定量をみつけたい、というモチベーションになります。

最大対数尤度

$E_{g}[\log{f_{m}(y|\hat{\theta_{m}})}]$ はどのように推定すればよいでしょうか？
真の分布について手元にある情報は、得られたデータ$y_{i}(i=1,\cdots,N)$です。よって、真の分布$g(y)$を次のような分布$\hat{g}(y)$で推定します。

\hat{g}(y)=\frac{1}{N}\quad (y=y_{1},\cdots,y_{N})

つまり、観測されたデータ$y_{i}(i=1,\cdots,N)$の各点に均等に確率$\frac{1}{N}$を割り振ります。
このとき、この推定された分布$\hat{g}(y)$でとった対数尤度 $\log{f_{m}(y|\theta_{m})}$ の期待値は次のようになります。

E_{\hat{g}}[\log{f_{m}(y|\hat{\theta_{m}})}]=\sum_{i=1}^N{\hat{g}(y_{i})\log{f_{m}(y_{i}|\hat{\theta_{m}})}}=\frac{1}{N}\sum_{i=1}^N{\log{f_{m}(y_{i}|\hat{\theta_{m}})}}

期待値の記号は、$E_{\hat{g}}[*]$ （$\hat{g}$による期待値）であることに注意です。

右辺の $\frac{1}{N}\sum_{i=1}^N{\log{f_{m}(y_{i}|\hat{\theta_{m}})}}$ を詳しくみてみましょう。

和の部分 $\sum_{i=1}^N{\log{f_{m}(y_{i}|\hat{\theta_{m}})}}$ は、統計モデル $f_m$ の対数尤度のパラメータに最尤推定量 $\hat{\theta}$ を代入した値になっています。
最尤推定量は対数尤度を最大にするパラメータの値なので、そのパラメータを代入した$\sum_{i=1}^N{\log{f_{m}(y_{i}|\hat{\theta_{m}})}}$ は対数尤度の最大値、つまり最大対数尤度と呼ぶことができます。

最大対数尤度による推定とバイアス

前項の議論から、次のことが見えてきました。

求めたいもの：平均対数尤度 $E_{g}[\log{f_{m}(y|\hat{\theta_{m}})}]$
問題点：真の分布 $g(y)$ での期待値は求めることができない
解決法：最大対数尤度 $\sum_{i=1}^N{\log{f_{m}(y_{i}|\hat{\theta_{m}})}}$ を使用して推定する

もし、最大対数尤度が平均対数尤度の推定量としてバイアスがなければ、次の式が成り立ちます。

E_{g} \left[\frac{1}{N} \sum_{i=1}^N \log{f_{m}(y_{i}|\hat{\theta}_{m})} - E_{g} \left[\log{f_{m}(y|\hat{\theta}_{m})} \right] \right]=0

便宜上、両辺を$N$倍して次のように書きかえます。

E_{g} \left[\sum_{i=1}^N \log{f_{m}(y_{i}|\hat{\theta}_{m})} - N \left( E_{g} \left[\log{f_{m}(y|\hat{\theta}_{m})} \right] \right) \right]=0

言葉で書くと、次のようになります。

E_{g} \left[ 最大対数尤度 - N\times 平均対数尤度 \right]=0

バイアスがないなら、平均対数尤度（の$N$倍）を推定するのに最大対数尤度をそのまま使用することは妥当だ、ということになります。

しかし残念ながら、この関係は成り立ちません。 最大対数尤度は平均対数尤度（の$N$倍）の推定量としてバイアスがある ことがわかっています。

その理由の直観的な理解は以下です。

最大対数尤度のパラメータ（最尤推定値）は手元のデータ $y_{i}(i=1,\cdots,N)$ を用いて求められる
最大対数尤度そのものも、手元のデータ $y_{i}(i=1,\cdots,N)$ を用いて求められている
パラメータの推定と対数尤度の計算の両方に、手元のデータ $y_{i}(i=1,\cdots,N)$ を使っている
データの二度使いによって、バイアスが生じる

機械学習の文脈でいうと、学習データを使ってモデルを評価しているのと同様の状況といえそうです。

バイアス補正（AICの意味）

残念ながら、最大対数尤度は平均対数尤度（の$N$倍）の推定量としてバイアスがあることがわかりました。
であれば、そのバイアスがどれくらいか分かれば、それを使って補正してあげればよいことになります。

言葉で書くと、次のようになります。

E_{g} \left[ (最大対数尤度 -バイアス補正項) - N\times 平均対数尤度 \right]=0

このバイアス補正項を近似的に求めて、それによってモデル選択の基準を与えたのがAICです。

\begin{align}
\mathrm{AIC} & =-2(最大対数尤度)+2(バイアス補正項) \\
& = -2(最大対数尤度)+2(統計モデルの自由度パラメータ数)
\end{align}

$-2$倍しているのは見た目を整える意味だと思います（特に正規分布の場合を意識？）
また最大対数尤度にマイナスがかかっているので、AICが小さいほうが良いモデルという意味になっています、

バイアス補正項が、なぜ近似的に「2×(統計モデルの自由度パラメータ数)」となるかの証明は、こちらの論文がわかりやすかったです。
統計学におけるモデル―情報量基準の観点から(山口健太郎) https://repository.kulib.kyoto-u.ac.jp/dspace/bitstream/2433/56989/1/yamaguchi.pdf

いったんまとめ

以上のとおり、AICの式で出てくる「モデルの自由度パラメータ数」は、
最大対数尤度を平均対数尤度（の$N$倍）の推定量として使用する場合のバイアス補正項
であることが整理できました。
最後に簡単なシミュレーションによって、最大対数尤度による推定のバイアスについて確認してみたいと思います。

シミュレーションデータでの検証

検証する点

以下の2点を検証します。

最大対数尤度は平均対数尤度（の$N$倍）の推定量としてバイアスがあるのか
データの二度づかいをしなければ平均対数尤度（の$N$倍）をバイアスがなく推定できるのか

シミュレーションコード

コードの全体を示します。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from scipy.stats import norm

# シミュレーションの設定
np.random.seed(42)  # 再現性のためのシード設定
n_simulations = 1000  # シミュレーション回数
n_samples = 100  # サンプルサイズ
true_coefficients = [2, -3, 5, 4]  # 真の線形回帰係数（切片, x1, x2, x3）
true_sigma_squared = 1  # 真の標準偏差の2乗

# 結果の記録
log_likelihoods = []
true_log_likelihoods = []
bias_values = []
bias_new_data_values = []

# 説明変数（x1, x2, x3）を生成
x1 = np.random.normal(0, 1, n_samples)
x2 = np.random.normal(0, 1, n_samples)
x3 = np.random.normal(0, 1, n_samples)

# 説明変数の結合 (切片項を含む)
X = np.column_stack((x1, x2, x3))
X = np.hstack((np.ones((n_samples, 1)), X))  # 切片項（1の列）を追加

# シミュレーション開始
for _ in range(n_simulations):

    # 応答変数 y の生成（真のモデル）
    noise = np.random.normal(0, np.sqrt(true_sigma_squared), n_samples)
    y = (true_coefficients[0] + true_coefficients[1] * x1 +
         true_coefficients[2] * x2 + true_coefficients[3] * x3 + noise)

    # 線形回帰モデルの学習
    model = LinearRegression(fit_intercept=False)
    model.fit(X, y)

    # 予測値と残差の計算
    y_pred = model.predict(X)
    residuals = y - y_pred

    # 推定された最大対数尤度
    sigma_squared = np.var(residuals)
    log_likelihood = -0.5 * n_samples * (np.log(2 * np.pi * sigma_squared) + 1)
    log_likelihoods.append(log_likelihood)

    # 真の対数尤度（真のパラメータを使用）
    true_log_likelihood = -0.5 * n_samples * (np.log(2 * np.pi * true_sigma_squared) + 1)
    true_log_likelihoods.append(true_log_likelihood)

    # バイアスの計算（同じデータを使用）
    bias = log_likelihood - true_log_likelihood
    bias_values.append(bias)

    # 新たにサンプリングしたデータを使用してバイアスを計算
    x1_new = np.random.normal(0, 1, n_samples)
    x2_new = np.random.normal(0, 1, n_samples)
    x3_new = np.random.normal(0, 1, n_samples)
    
    X_new = np.column_stack((x1_new, x2_new, x3_new))
    X_new = np.hstack((np.ones((n_samples, 1)), X_new))

    noise_new = np.random.normal(0, np.sqrt(true_sigma_squared), n_samples)
    y_new = (true_coefficients[0] + true_coefficients[1] * x1_new +
             true_coefficients[2] * x2_new + true_coefficients[3] * x3_new + noise_new)

    y_pred_new = model.predict(X_new)
    residuals_new = y_new - y_pred_new

    sigma_squared_new = np.var(residuals_new)
    log_likelihood_new = -0.5 * n_samples * (np.log(2 * np.pi * sigma_squared_new) + 1)
    bias_new_data = log_likelihood_new - true_log_likelihood
    bias_new_data_values.append(bias_new_data)

# 箱ひげ図のプロット
plt.figure(figsize=(10, 6))
plt.boxplot([bias_values, bias_new_data_values], labels=['same data', 'new data'])
plt.title('box plot-bias')
plt.ylabel('bias')
plt.grid(True)
plt.show()

シミュレーションの設定を簡単に記します。

真のモデル：$y=2-3x_{1}+5x_{2}+4x_{3}+\epsilon_{i},\quad \epsilon_{i}～N(0, 1)$
- 説明変数が3つのガウス型線形回帰モデルです
サンプルサイズ：100
シミュレーション回数：1000回
シミュレーションの流れ：
1. 真のモデル式から $y_{i}(i=1,\cdots,100)$ を生成
2. $y_{i}(i=1,\cdots,100)$ をもとに最尤推定
3. 最尤推定されたモデルをもとに $y_{i}(i=1,\cdots,100)$ の対数尤度（最大対数尤度）を計算
4. 真のモデルの対数尤度（平均対数尤度のN倍）を計算
5. 新たに $y_{i}^{new}(i=1,\cdots,100)$ を生成して対数尤度を計算（データを二度づかいせずに新たなデータで計算した対数尤度）
6. 推定値と真値との差をもとにバイアスを確認

シミュレーション結果

シミュレーション結果を箱ひげ図で表示します。

左側が最大対数尤度のバイアスの箱ひげ図です。
右側は最尤推定に使用したデータとは別のデータで、対数尤度を計算したときのバイアスの箱ひげ図です。

中央値は左側の箱ひげ図が正の方にずれていますし、箱の位置・箱ひげ図全体の位置を見てもやはり左側の箱ひげ図が正の方にずれているのがよくわかります。
一方右側の箱ひげ図は0を中心におおむね上下対称になっています。

シミュレーションでも、最大対数尤度がバイアスをもっていそうなことが確認できました。
また、新たなデータを使用すればバイアスなく対数尤度を推定できそうであることも示唆されます。

まとめ

本記事ではAICの理解の確認を目的として、最大対数尤度のバイアス補正というAICの成り立ちを改めて確認しました。また簡単なシミュレーションで、実際のバイアスについても検証しました。

シミュレーションをみると、クロスバリデーションの有用性も改めて実感できます。やはりデータの二度づかいは分析する上で常に気を配らないといけない点なのですね。（そしてそれをある程度補正してくるAICの有用性も）

この記事を読んでくださった方々の、AICの説明がかっちょよくなることに少しでも貢献できていれば幸いです！

参考文献

小西貞則, 2010, 多変量解析入門, 岩波書店 https://www.iwanami.co.jp/book/b265441.html
山口健太郎, 2008, 統計学におけるモデル：情報量基準の観点から, 京都大学文学部科学哲学化学史研究室 https://repository.kulib.kyoto-u.ac.jp/dspace/bitstream/2433/56989/1/yamaguchi.pdf