この記事を読んで得られること

ダブルクロスバリデーションとは？
他の交差検証手法との比較
ダブルクロスバリデーションの実装方法

「ホールドアウト法」、「交差検証」、「クロスバリデーション」 についても学べると思いますので参考にしてください。

※プログラミング関係の内容を他にも投稿していますので、よろしければこちらの一覧から他の投稿も見て頂けますと幸いです。

ダブルクロスバリデーションとは？

機械学習では未知のデータに対する予測精度を評価するために交差検証を行います。ダブルクロスバリデーションはその手法の1つであり、各交差検証の手法を以下にまとめました。

1. 良くない例

まずは良くない例として、データセットを分割せずに学習用データ＝検証用データの例を示しています。これはモデルの学習に用いたデータを使って、そのモデルを評価するというものです。検証用のデータは全て学習済みのデータ（見たことがあるデータ）であるので、これをやってしまうと未知のデータに対する予測性能（汎化性能）が全く測れません。

2. ホールドアウト検証

こちらは最も手軽な汎化性能の検証手法です。まずはデータセットを学習用データ（訓練データ）と検証用データ（テストデータ）に分割し、学習用のデータをモデルに学習させた上で、検証用のデータを使ってモデルの性能を評価します。検証用のデータはモデルにとって見たことのない未知のデータなので、汎化性能を評価することができます。しかし、学習用データと検証用データに偏りがある際はモデルの性能を上手く評価できないという問題があります。

3. k分割交差検証

ホールドアウト法で問題となるデータの偏りをなくそうとしたものがこのk分割交差検証です。この手法ではデータセットをk分割し（図の例では4分割）、1つを検証用データとして残りのデータを学習用データとしてモデルを評価します。この操作をk回繰り返し、それぞれの検証で得られたモデルの精度（決定係数など）を平均化した値をそのモデルの精度とすることが一般的です。ちなみに分割数kをデータ数まで増やした場合は**Leave-One-Out法（LOO法）**と呼ばれます。

4. グリッドサーチクロスバリデーション

予測モデルによってはハイパーパラメータの調整が必要になります。例えばSVMでは、Cやガウシアンカーネルのγといったハイパーパラメータを決定する必要があります。ハイパーパラメータの決定は一般的にはクロスバリデーションによって決定されます。また、ハイパーパラメータが複数ある場合はグリッドサーチクロスバリデーションによって指定したハイパーパラメータの全ての組み合わせに対して網羅的に学習を行い最も良い精度を示したハイパーパラメータを採用することが一般的です。

5. ダブルクロスバリデーション

Nested Cross Validationとも呼ばれます。"Nested"とは"入れ子になった"という意味です。これと比較する格好でグリッドサーチクロスバリデーションはNon-nested Cross Validationとも呼ばれます。この手法では交差検証を内側のクロスバリデーション(CV)と外側のクロスバリデーションに分け、内側のCVではハイパーパラメータの選択に注力し、選択されたハイパーパラメータを用いで外側のCVでモデルの精度を評価します。図の例では外側のCVの分割数は4なので、4通りのハイパーパラメータ（モデル）が抽出されます。それぞれの予測モデルから得られた予測値からモデル全体の精度を求めます。

どういうときにダブルクロスバリデーションを使うのか？

こちらのページでは「サンプルが少ないときでも、適切にモデルの推定性能を検証するために、ダブルクロスバリデーションが利用されます」との記述があります。
個人的には次の理由も考えています。グリッドサーチクロスバリデーションはモデルの予測精度が過大評価されていると考えます。検証用データに対してハイパーパラメータをグリッドサーチしています。つまり、予測モデルは一度検証用データに触れており、その評価結果がベストになるハイパーパラメータを最終的なモデルのハイパーパラメータとして選択しています。それに対して、ダブルクロスバリデーションでは精度を算出するための検証用データは予測モデルに1回も触れておりません。（ハイパーパラメータは内側のクロスバリデーション内で決定され、モデルの評価は一度も触れていない外側のデータに対して行われます）よって、ダブルクロスバリデーションから算出されたモデルの予測精度はグリッドサーチクロスバリデーションから算出された予測精度よりも厳しいものであり、より現実に即していると考えています。（間違っていたらご指摘ください）ただ、そこはデータ数のバイアスもあるのでダブルクロスバリデーションが必ずしも良いとは言い切れないのも悩ましいです…（こちらの「疑問と悩み」に上手く表現されています）

また、こちらのページではクロスバリデーションとダブルクロスバリデーションの違いについて説明しており、クロスバリデーションはハイパーパラメータの決定、ダブルクロスバリデーションは手法の選択(SVRかPLSかを決める)に用いるとされています。

ダブルクロスバリデーションの実装

sklearnの公式を見てもこちらやこちらのサイトに実装方法が載っていますし、その他にも検索すれば様々なページが見つかります。これらのコードでは学習用データもモデルの評価に含まれています。（以下の図の②精度に該当します）しかし、私は検証用データのみの評価（以下の図の①精度に該当します）を積み上げてモデルを評価したいと考えています。

そこで、検証用データのみの評価を積み上げてモデル全体の精度を評価するコードを以下で紹介します。

環境

windows10
anaconda 0.0.1.1
python 3.7.13
scikit-learn 1.0.2

今回のデータセットはこちらのページで紹介したcompound_2_edit.csvを用います。compound_2_edit.csvのデータセットはこちらのGithubに保存しています。

# データの読み込み
import pandas as pd

df = pd.read_csv("compound_2_edit.csv", encoding="shift_jis", index_col=0)
df.head()

df.shape

出力

(229, 8)

229個の化合物に対して、8つの物性値を持ったデータセットです。ここでは、"band_gap"を目的変数として、以下で分析を進めていきます。今回はSVRでモデルを作成しようと思うので、以下で説明変数を標準化していきます。

x = df.iloc[:, 1:]  # 最初の列のbandgapを目的変数とする
y = df.iloc[:,1]
autoscaled_x = (x - x.mean()) / x.std()  # オートスケーリング
autoscaled_x

# 標準化した説明変数xと目的変数yを結合して新たなデータフレームを作る
df = autoscaled_x
df["y"] = y
df

分析していく新たなデータセットができました。以下ではSVRを予測モデルとして、ダブルクロスバリデーションを実行していきます。なお、ハイパーパラメータの探索範囲はこちらを参考にしました。

import numpy as np
from sklearn.model_selection import KFold
from sklearn import svm
from sklearn.model_selection import GridSearchCV

# グリッドサーチのスコア、グリッドサーチの最適モデル、yの予測値、yの観測値、最適モデルでのR2をそれぞれ定義する
gs_scores = []
best_models = []
df_y_preds = pd.DataFrame()
df_y_observes = pd.DataFrame()
r2_scores = []

# 外側と内側のCVの分割数を定義する
outer_cv = 4
inner_cv = 4

# 外側のクロスバリデーションをKFoldで設定する（5分割する）
kfold = KFold(n_splits = outer_cv)
for train, test in kfold.split(df):
    X_train = df.iloc[train,:-1]
    y_train = df.iloc[train, -1]
    X_test = df.iloc[test,:-1]
    y_test = df.iloc[test, -1]
    
    # 内側のクロスバリデーションをGridSearchCVで設定する（4分割する）
    # クロスバリデーションの条件設定
    C_range = list(2 ** np.arange(-5, 11, dtype=float))
    epsilon_range = list(2 ** np.arange(-10, 1, dtype=float))
    gamma_range = list(2 ** np.arange(-20, 11, dtype=float))
    param_grid = {"C":C_range, "epsilon":epsilon_range, 
                  "gamma":gamma_range}
    # モデルの定義
    model = svm.SVR(kernel="rbf")
    # グリッドサーチの設定
    gs = GridSearchCV(estimator=model, param_grid=param_grid, scoring="r2", cv=inner_cv, n_jobs=-1, return_train_score=True)
    # グリッドサーチの実行
    gs.fit(X=X_train, y=y_train)
    # グリッドサーチの最適モデルでの平均スコアの定義
    gs_score = gs.best_score_
    gs_scores.append(gs_score)
    # 最適モデルの定義
    best_model = gs.best_estimator_
    print(best_model)
    best_models.append(best_model)
    # 最適モデルで予測
    y_pred = best_model.predict(X_test)
    df_y_preds = pd.concat([df_y_preds, pd.DataFrame(y_pred)])
    # 教師データ（観測値）
    df_y_observes = pd.concat([df_y_observes, y_test])
    # 構築したモデルのスコア
    r2_score = best_model.score(X_test, y_test)
    r2_scores.append(r2_score)

出力

SVR(C=1024.0, epsilon=0.0009765625, gamma=1.52587890625e-05)
SVR(C=512.0, epsilon=0.0009765625, gamma=3.0517578125e-05)
SVR(C=512.0, epsilon=0.0009765625, gamma=6.103515625e-05)
SVR(C=1024.0, epsilon=0.0009765625, gamma=1.52587890625e-05)

ダブルクロスバリデーションで最適化した各モデルのハイパーパラメータを整理してまとめておきます。

# 結果をまとめておく
df_result = pd.concat([pd.DataFrame(gs_scores), pd.DataFrame(best_models).iloc[:,0], pd.DataFrame(r2_scores)], axis=1)
df_result.columns = ["Gridserach_score", "model", "r2"]
# 結果をcsvファイルに書き出す
df_result.to_csv("wcv_result.csv", index=False)
# データフレームで結果の表示
df_result

次に予測値と観測値をまとめておきます。

# df_observesとdf_y_predは分割ごとにindextがリセットされるので連番にしておく
df_y_observes = df_y_observes.reset_index(drop=True)
df_y_preds = df_y_preds.reset_index(drop=True)

# 実測値、予測値の順にデータフレームに格納し、化合物のラベルをふる
df_y_obs_pred = pd.concat([df_y_observes, df_y_preds], axis=1)
df_y_obs_pred.columns = ["y_observed", "y_predicted"]
df_y_obs_pred.index = df.index

df_y_obs_pred

あとでモデルを分析する際に備えてcsvファイルに保存しておくと便利です。

# csvファイルに保存
df_y_obs_pred.to_csv("yosoku_SVR4.csv", index=False)

結果を可視化していきます。可視化はこちらのコードを参考にしています。

from matplotlib import pyplot as plt
from sklearn.metrics import r2_score
from sklearn.metrics import mean_squared_error

# グラフのラベルを定義
graph_title = "target:band_gap,  model:Support Vector reg."

# yyplot 作成関数
def yyplot(y_obs, y_pred):
    yvalues = np.concatenate([y_obs.flatten(), y_pred.flatten()])
    ymin, ymax, yrange = np.amin(yvalues), np.amax(yvalues), np.ptp(yvalues)
    fig = plt.figure(figsize=(8, 8))
    plt.scatter(y_obs, y_pred)
    plt.plot([ymin - yrange * 0.01, ymax + yrange * 0.01], [ymin - yrange * 0.01, ymax + yrange * 0.01])
    plt.xlim(ymin - yrange * 0.01, ymax + yrange * 0.01)
    plt.ylim(ymin - yrange * 0.01, ymax + yrange * 0.01)
    plt.xlabel('y_observed', fontsize=24)
    plt.ylabel('y_predicted', fontsize=24)
    plt.title(graph_title, fontsize=24)
    plt.tick_params(labelsize=16)
    plt.show()

    return fig

# yyplot の実行例
np.random.seed(0)
y_obs = np.random.normal(size=(1000, 1))
y_pred = y_obs + np.random.normal(scale=0.3, size=(1000, 1))
# dataframe→numpy配列にするために.valuesを実行
fig = yyplot(df_y_observes.values, df_y_preds.values)

# 決定係数の算出
r2 = r2_score(df_y_observes.values, df_y_preds.values)
print("決定係数：" + '{:.3f}'.format(r2))

# RMSEの算出
mse = mean_squared_error(y_true=df_y_observes.values, y_pred=df_y_preds.values)
print("RMSE：" + '{:.2f}'.format(np.sqrt(mse)))

このように評価することで、外側のクロスバリデーションからのみでモデルの予測精度を導くことができます。（ただ、今回は予測が容易すぎるデータセットでしたね…反省）