AIC(赤池情報量基準)を使った線形回帰のモデル評価 - 機械学習に詳しくなりたいブログ

これまで線形回帰における過学習の抑制方法をいくつか見てきました。

参考

これらの考え方は、パラメータ $\mathbf w$ が大きくなりすぎることに対して制約を与えるというものでした。それに対してAIC(Akaike's Information Criterion：赤池情報量基準)は、モデルの複雑さに対してペナルティを与えて過学習を防止しようという考え方です。そもそもモデルを複雑にしすぎなければ、過学習は起きないでしょうということですね。AICは以下の式で定義されます。

$\mathrm{AIC} = -2\ln L + 2M \tag{1}$

ここで $M$ はモデルのパラメータ数、 $L$ は最大尤度です。そして $\mathrm{AIC}$ が最小になるようなモデルを選択します。原論文を見てみましたが、ちょっと僕には荷が重そうなので公式として使ってみることにします。

$y = -\ln x$ のグラフは下図のような形です。

尤度の増加に伴って $-2\ln L$ は小さくなりますが、その減少速度はどんどん小さくなります。一方で $2M$ は線形に増加していきますので、そのバランスがとれたところをモデル数として採用するっていうイメージでしょうか。

さて、線形回帰において誤差を正規分布として仮定すれば、尤度は

$\displaystyle L = \prod_{i=1}^{N}\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left\{ \frac{ -( t_i - \mathbf{w}^{T} \boldsymbol{\phi}(x_i))^2 }{ 2 \sigma ^2}\right\} \tag{2}$

でした。（参考：線形回帰を最尤推定で解く（尤度とは？））よって、

$\displaystyle \ln L = -\frac{N}{2}\ln(2 \pi \sigma^{2}) + \sum_{i=1}^{N} \left\{ \frac{-(t_i - \mathbf{w}^{T} \boldsymbol {\phi}(x_i))^2}{2 \sigma^{2}} \right\} \tag{3}$

です。ここで残差平方和 $\displaystyle \sum_{i=1}^{N} (t_i - \mathbf{w}^{T} \boldsymbol{\phi}(x_i))^{2}$ を $\mathrm{RSS}$ とおくと、

$\displaystyle \ln L = -\frac{N}{2}\ln(2 \pi \sigma^{2}) - \frac{\mathrm{RSS}}{2 \sigma^{2}}\tag{4}$

と書けます。 $\sigma^{2}$ は誤差の分散ですから、その最尤推定量は

$\begin{eqnarray*} \displaystyle \sigma^{2} &=& \frac{1}{N} \sum_{i=1}^{N} (t_i - \mathbf{w}^{T} \boldsymbol{\phi}(x_i))^{2} \tag{5} \\ &=& \frac{\mathrm{RSS}}{N} \tag{6} \end{eqnarray*}$

です。（誤差の平均は0と仮定しています。分散の定義式に従って計算したものが正規分布の最尤推定量です）　これを式(4)に代入すれば

$\displaystyle \ln L = -\frac{N}{2} \ln \left( \frac{2 \pi \mathrm{RSS}}{N} \right) - \frac{N}{2} \tag{7}$

よって

$\begin{eqnarray*} \mathrm{AIC} &=& -2 \ln L + 2M \tag{8} \\ \displaystyle &=& N \ln \left( \frac{2 \pi \mathrm{RSS}}{N} \right) + N + 2M \tag{9} \\ \displaystyle &=& N \ln 2\pi + N \ln \frac{\mathrm{RSS}}{N} + N + 2M \tag{10} \end{eqnarray*}$

となります。ここで $N \ln 2\pi,N$ は定数ですから、

$\displaystyle N \ln \frac{\mathrm{RSS}}{N} + 2M \tag{11}$

を最小にすればよいことになります。

さて、式(11)を使って最適なモデルパラメータ数を計算してみました。近似モデルはM次多項式で、理想曲線はsin関数です。結果は以下のグラフです。 f:id:opabinia2:20180224122518p:plain

最適なパラメータはM=5となり、これ以上増やしても $\mathrm{AIC}$ は悪化していく一方です。以下の図は左がM=5、右がM=20としたときの近似曲線です。 f:id:opabinia2:20180224122054p:plain

今回のコードです。訓練データ数によって最適パラメータ数がかなりばらつきました。どう見ても過学習だろうと思えるような曲線が最適となったり、あまり安定しません。使い方に問題があるのでしょうか？

# AICで線形回帰モデルの最適パラメータ数を求める

import numpy as np
import matplotlib.pyplot as plt


# y = w0*x^0+....wM*x^M　を、引数xの配列数分求める
def y(w, x, M):
    X = np.empty((M + 1, x.size))
    for i in range(M + 1):
        X[i, :] = x ** i
    return np.dot(w.T, X)


# ランダムシードを固定
np.random.seed(0)

# 訓練データ数
N = 50

# AICを求めるパラメータ数
MAX_N = 20

# 訓練データの列ベクトル
x = np.linspace(0, 1, N).reshape(N, 1)
# 訓練データtの列ベクトル
t = np.sin(2*np.pi*x.T) + np.random.normal(0, 0.2, N)
t = t.reshape(N, 1)

AIC_array = np.empty([0])

for M in range(1, MAX_N+1):

    # 行列Phiを作成
    Phi = np.empty((N, M+1))
    for i in range(M+1):
        Phi[:, i] = x.reshape(1, N) ** i

    # 係数wの列ベクトルを解析的に求める
    w = np.linalg.solve(np.dot(Phi.T, Phi), np.dot(Phi.T, t))

    # 求めた係数wを元に、新たな入力x2に対する予測値yを求める
    x2 = np.linspace(0, 1, 100)
    y1 = y(w, x2, M)

    # 残差平方和を求める
    RSS = np.sum((y(w, x.T, M).T - t)**2)

    # AICを求める
    AIC = N*np.log(RSS/N) + 2*M
    print(M, RSS, AIC)
    AIC_array = np.append(AIC_array, AIC)

    # 結果の表示
    # plt.xlim(0.0, 1.0)
    # plt.ylim(-1.5, 1.5)
    # plt.scatter(x, t)
    # plt.plot(x2, y1.T)
    # plt.show()

plt.xlabel("M")
plt.ylabel("AIC+const")
plt.plot(range(1, MAX_N+1), AIC_array)
plt.show()