AICを使った変数選択　ステップワイズ

AICを使った変数選択の方法として、ステップワイズ法を学習します。
ステップワイズ法は、「たくさんの説明変数の中から汎化性能の高い説明変数の組み合わせを求める」アルゴリズムで、勘によらずに良いモデルを求められます。

ステップワイズ法は、以下の種類があります。

変数増加法：空の状態から変数を追加していく方法
変数減少法：全て選んだ状態から削除していく方法
変数増減法：空の状態から追加と削除を繰り返す方法
変数減増法：全て選んだ状態から削除と追加を繰り返す方法
ここでは、よく使われる変数増減法をやってみます。

ここで紹介するステップワイズ法では、AICを評価規準にします（他の評価規準を使う方法もあります）。
また、評価値の計算では、クロスバリデーションの平均値を使います。
同じ説明変数の組合せ（sel）に対し、calc_aic(sel)を実行するので、lru_cacheを指定してキャッシュします。

# 前々回のプログラムの読込
%run 1.ipynb

%matplotlib inline
import matplotlib.pyplot as plt
from sklearn.model_selection import cross_val_score
from functools import lru_cache

# ステップワイズ法のサブルーチン
def _calc_best(calc_aic, best, N, is_add):
    found = False
    for i in range(N):
        key = (best[1] | {i}) if is_add else (best[1] - {i})
        v = calc_aic(tuple(key))
        if v < best[0]:
            found, best = True, (v, key)
    return found, best

# ステップワイズ法
def stepwise(X, y, fit_intercept=True, cv=None):
    @lru_cache(maxsize=1024)
    def calc_aic(sel):
        X_ = X[:, sel]
        return cross_val_score(model.fit(X_, y), X_, y,
                               scoring=AIC, cv=cv).mean()
    model = LinearRegression(fit_intercept=fit_intercept)
    N = X.shape[1]  # 説明変数数
    best = 1e308, set()  # 最良のAICと選択
    found = True
    while found and len(best[1]) < N:
        found, best = _calc_best(calc_aic, best, N, True)
        found2 = True
        while found2 and len(best[1]) > 1:
            found2, best = _calc_best(calc_aic, best, N, False)
            found |= found2
    return best[0], list(best[1])

ステップ 1. 説明変数候補を空にします。
ステップ 2. 説明変数候補に新たな説明変数を1つ追加したとき、最も評価値（AIC）が良いものがあれば、追加します。
ステップ 3. 説明変数候補から説明変数を1つ削除したとき、最も評価値（AIC）が良いものがあれば、削除します。
ステップ 4. 削除が可能な間、ステップ3を繰り返します。
ステップ 5. 更新がないか、全ての説明変数を選んだ場合は終了します。そうでない場合、ステップ2へ行きます。

関数calc_bestは、ステップワイズ法のサブルーチンで、上記ステップの2と3に対応します。頭に「」をつけた名前にして、非公開であることを表しています。

関数stepwiseが、ステップワイズ法に対応します。引数は以下の通りです。

X：全ての説明変数（numpy.ndarray）
y：目的変数
fit_intercept：y切片の有無（LinearRegressionの引数）
cv：クロスバリデーションの分割数（cross_val_scoreの引数）
返り値は、AICの値と選択された説明変数のインデックスのリストです。

# ボストン市住宅価格データで変数選択
aic, sel = stepwise(X.to_numpy(), y, cv=6)
aic, sel

(517.1135458604272, [10, 11, 12, 5])

aic, sel = stepwise(X.to_numpy(), y, cv=6)で、ボストン市住宅価格データに対し、ステップワイズ法を実行しています。

# 真値と予測値の散布図
X_sel = X.iloc[:, sel]
lr = LinearRegression().fit(X_sel, y)
y_pred = lr.predict(X_sel)
plt.xlabel('y_pred')
plt.ylabel('y_true')
plt.plot(y_pred, y, '.');

# 選択された説明変数の列名
X_sel.columns

Index(['PTRATIO', 'B', 'LSTAT', 'RM'], dtype='object')

散布図で見ると、良い結果になっていることが確認できます。
X_sel.columnsで、選択された説明変数の列名がわかります。PTRATIO, B, LSTAT, RMの4つが選ばれていることがわかります。