赤池情報量基準(AIC)を多項式の回帰で実装

こんにちは、さくまっちょです。

今回は本業のデータサイエンス関連で気になっていたこととして赤池情報量基準(AIC)を実際に実験用にデータを作って計算しその意味や威力を確認してみました。

機械学習・データサイエンスやっているとモデルの汎化性能を上げるために複雑なモデルを使わない方が良いという説明が PRMLやはじパタなど各種機械学習系のテキストには必ずと言っていいほど載っています。
数式もちゃんと載っているのですが、いざ自分で計算してみようと思っても実装をしている記事などがなかなか出てこず、仕方がないので数式を調べながら実装してみました。

同じ疑問をお持ちの方のお役に立てば幸いです。

もしよければコメントいただきたいのですが普段みなさんってAICって計算しないのでしょうか？
もしくはいくつかのライブラリや計算クラスではAICやそれに準ずるモデル選択のメソッドが用意されているから自分では実装しないのでしょうか？

すべての実装についてはこちらの僕のノートブックをご覧ください。

対象とする読者

機械学習エンジニア
データサイエンティスト
機械学習・データサイエンスに興味のある学生・エンジニア

前提とする知識

回帰分析がわかること
基本的な数学(sin, logや線形代数)

バージョン情報

Python 3.7.3

赤池情報量基準(AIC)とは
多項式回帰での次数と汎化性能
データの分布が正規分布を仮定した場合のAIC
- AICの定義
- データの分布が正規分布を仮定した場合のAIC
AICを使った多項式回帰の次数の決定
まとめ

赤池情報量基準(AIC)とは

Wikipediaによると

元統計数理研究所所長の赤池弘次が1971年に考案し1973年に発表した
統計モデルの良さを評価するための指標
「モデルの複雑さと、データとの適合度とのバランスを取る」ために使用される
AIC最小のモデルを選択すれば、多くの場合、良いモデルが選択できる

です。また、詳細な意味や議論については赤池氏の記事をぜひ読んでみてください。
よく比較されるBICとの違いについてもご本人が議論されています。

詳細な議論はともかくとして、AICを使うことによって様々に作ったモデルの中から「これが一番データに合ってそうだ」というモデルを決めるのに役立ちます。

多項式回帰での次数と汎化性能

PRMLやはじパタで説明されているように、僕もまずはデータを生成し多項式で回帰して汎化性能を確かめてみました。

作成したデータがこちらです。

f:id:sakumadaisuke32:20210523171654p:plain — 図1. 実験データ. 実線が理想的な正弦波. プロットが実際に訓練データとして使うために正弦波にノイズを与えたデータ.

これに0次から20次までの多項式で回帰した結果がこちらです。

f:id:sakumadaisuke32:20210523172024p:plain — 図2. 多項式で回帰した結果. 黒実線が理想的な正弦波. 緑実線が回帰した結果. 赤点が訓練データ. 青点が検証データ.

6次くらいまでは正弦波っぽく描かれていますが、7次以降は赤点(訓練データ)との距離は縮んでそうですが歪みがあるのと青点(検証データ)との距離は離れてそうです。

実際に訓練データと検証データの回帰線とのRMSEを比較したのがこちらです。

f:id:sakumadaisuke32:20210523172436p:plain — 図3. 多項式の次数とRMSE

3次までは訓練・検証ともRMSEは下がっていますね。
それ以降も横ばいな感じはありますが、訓練データは10次まで、検証データも9次まではRMSEが改善してそうです。
このグラフから得られる結果とグラフの形状から、なんとなく3次～6次の間でモデルを決定したいです。

ここで「なんとなく」ではなく定量的に次数を決めるための指標がAICです。

データの分布が正規分布を仮定した場合のAIC

AICの定義

Wikipediaによると

$\displaystyle{ \mathrm{AIC} = -2 \ln {L} + 2 k }$

$L$ は尤度、 $k$ はパラメータの数(今回は多項式の次数 + 1(定数項を含むため)) です。

尤度を計算するのに確率分布を導入しないといけなくて、このあたりの実装が厄介だったのですがこちらの41ページで正規分布の場合は簡単に整理できるという資料がありましたので、この通りに計算してみたのが次です。

データの分布が正規分布を仮定した場合のAIC

$\displaystyle{ \mathrm{AIC} = n \ln{2 \pi} + n \ln{\mathrm{MSE}} + n + 2(k + 2) }$

$n$ はデータの個数、 $\mathrm{MSE}$ は平均自乗誤差です。実装はこちらになります。

from sklearn.metrics import mean_squared_error as MSE

def AIC(ys_actual, ys_pred, num_feature):
    num_data = len(ys_actual)
    mse = MSE(ys_actual, ys_pred)
    return num_data * np.log(2 * np.pi) + num_data * np.log(mse) + num_data + 2 * (num_feature + 2)

いざ書いてみたらシンプルでした。