線形基底関数モデル⑤ K-分割交差検証

線形基底関数モデル⑤ K-分割交差検証

K-分割交差検証

 前回記事で扱ったホールドアウト検証は、入力データと目標データを一つの方法で分割するので、結果に偏りが生じる可能性があります。そこで今回は下図のようにデータを 4 等分して、テストデータと訓練データの役割を変えながら、近似曲線とテストデータの誤差を計算し、最後に 4 個の標準偏差の平均を求めて、これをモデルの評価基準とします。

 機械学習 交差検証の概念図

 標準偏差の平均が最小となる基底関数の組合わせが決定したら、すべての入力データ x と目標データ y を使ってパラメータベクトルを計算し直して、これを最終的なモデルとして採用します。

 一般に、データを K 分割して誤差を平均して、ばらつきの影響を軽減する手法を K-分割交差検証 (K-fold cross validation) といいます。
 

テストデータと訓練データに分割する方法

 交差検証のプログラムを書く前に、配列を分割する方法について解説します。
 numpy.array_split(x, K) を使うと、配列 x を K 個に分けることができます。x の要素数が K で割り切れるときは K 等分し、割り切れない場合は余りを先頭から振り分けます。入力データ x を 4 等分してみましょう。

# https://python.atelierkobato.com/validation/

# 線型基底関数モデル リストM3-C-1

# NumPyをインポート
import numpy as np

# 入力データ
x = np.array([35, 16, 22, 43, 5,
              66, 20, 13, 52, 1,
              39, 62, 45, 33, 8,
              28, 71, 24, 18, 3])

# 分割数K
K = 4

# 配列XをK分割
xd = np.array_split(x, 4)

print(xd)
[array([35, 16, 22, 43,  5]), array([66, 20, 13, 52,  1]),
array([39, 62, 45, 33,  8]), array([28, 71, 24, 18,  3])]

 返り値は配列を格納したリストなので、たとえば 3 番目の配列は xd[2] で参照できます。xd[2] をテストデータとして選択した場合、訓練データは元の配列 x から xd[2] の要素を除いたものになります。numpy.delete(x, slice(a, b)) を使うと、[a:b] でスライスした要素を配列 x から削除することができます。

# https://python.atelierkobato.com/validation/

# 線型基底関数モデル リストM3-C-2
# 動作条件:リストM3-C-1の実装

# xの要素数
N = x.shape[0]

# 3番目の配列をテストデータとして選択
x_test = xd[2]

# 訓練データの作成
x_train = np.delete(x, slice(2*int(N/K), 3*int(N/K)))

print("テストデータ : {}".format(x_test))
print("訓練データ : {}".format(x_train))
テストデータ : [39 62 45 33  8]
訓練データ : [35 16 22 43 5 66 20 13 52  1 28 71 24 18 3]

 今の場合、int(N/K) は int(20/4) = 5 なので、numpy.delete() の第 2 引数には [10:15] を指定していることになります。int() を使っているのは、N/K が非整数になる場合にも処理できるようにするためです。一般に xd[j] をテストデータに選んだ場合、訓練データは

np.delete(x, slice(j * int(N/K), (j + 1) * int(N/K)))

という記述で作成できます。for文を使って、j を 0 から K - 1 まで回せば、すべての組合わせを得ることができます。
 

K-分割交差検証によるパラメータの最適化

 K-分割交差検証のコードです。
 Fit_funcクラス (リストM3-B-1) の実装が前提となっています。

# https://python.atelierkobato.com/validation/

# 線型基底関数モデル リストM3-B-7
# 動作条件:リストM3-B-1の実装

# 入力データ
x = np.array([35, 16, 22, 43, 5,
              66, 20, 13, 52, 1,
              39, 62, 45, 33, 8,
              28, 71, 24, 18, 3])

# 目標データ
y = np.array([85.19, 58.93, 64.27, 68.91, 21.27,
              68.88, 60.07, 55.18, 88.08, 8.89,
              82.31, 81.18, 80.76, 78.98, 37.55,
              75.9, 72.39, 69.51, 62.04, 12.47])

# ★★★★★★★★★★

# データの数
N = x.shape[0]

# 分割数を設定
K = 4

# データをK分割
xd = np.array_split(x, K)
yd = np.array_split(y, K)

# ★★★★★★★★★★

# ガウス関数を定義
def gauss(x, mu = 0, sigma = 1):
    return np.exp(-(x - mu)**2 / (2 * sigma**2))

# 標準偏差の平均値sd_meanの初期値
sd_mean = 10.0

# 乱数シードを固定
np.random.seed(10)

# 試行回数(Number of trials)を設定
nt = 8000

# ★★★★★★★★★★

# 乱数を使った基底関数の調整
for k in range(nt):
    rdmu = 20 * np.random.rand(5) - 10
    rdsg = 10 * np.random.rand(5) - 5
    
    def phi_0(x):
        return gauss(x, rdmu[0], 10 + rdsg[0])
    
    def phi_1(x):
        return gauss(x, 20 + rdmu[1], 10 + rdsg[1])
    
    def phi_2(x):
        return gauss(x, 40 + rdmu[2], 10 + rdsg[2])

    def phi_3(x):
        return gauss(x, 60 + rdmu[3], 10 + rdsg[3])

    def phi_4(x):
        return gauss(x, 80 + rdmu[4], 10 + rdsg[4])

    def phi_5(x):
        return 1

    # 基底関数のリスト
    func_list = [phi_0, phi_1, phi_2, phi_3, phi_4, phi_5]

    # sd_testの値を初期化
    sd_test = 0

    for j in range(K - 1):
        x_test = xd[j]
        x_train = np.delete(x, slice(j * int(N/K), (j + 1) * int(N/K)))
        y_test = yd[j]
        y_train = np.delete(y, slice(j * int(N/K), (j + 1) * int(N/K)))

        # Fit_funcオブジェクトを作成
        w = Fit_func(x_train, y_train, func_list)

        # 現在選択しているテストデータの平均2乗誤差を計算
        mse_test = np.mean((w.line(x_test) - y_test)**2)

        # 現在選択しているテストデータの標準偏差を加算
        sd_test += np.sqrt(mse_test)

    # 標準偏差の平均値を計算
    sd_m = np.round(sd_test/K, 3)

    # 標準偏差の平均値と基底関数リストの更新
    if sd_m < sd_mean:
        sd_mean = sd_m
        bf_list = func_list
        
# ★★★★★★★★★★

# 採用された基底関数でFit_funcオブジェクトを作成
z = Fit_func(x, y, bf_list)

# Figureを作成
fig = plt.figure(figsize = (8, 6))

# 標準偏差
sd = np.round(z.sd(), 3)

# ★★★★★★★★★★

# FigureにAxesを1つ追加
ax = fig.add_subplot(111)

# Axesのタイトルを設定
ax.set_title("Random Gaussian base\nK-fold cross validation", fontsize = 16)

# テキストボックス書式辞書を作成
boxdic = {"facecolor" : "white", "edgecolor" : "gray",}

# テキストボックスを表示
ax.text(48, 10, "NT = {}\nSD = {}\nSD_mean = {}"
        .format(nt, sd, sd_mean),
        size = 14, linespacing = 1.5, bbox = boxdic)

# 目盛線を表示
ax.grid()

# 軸範囲を設定
ax.set_xlim([0, 70])
ax.set_ylim([0, 100])

# 軸ラベルを設定
ax.set_xlabel("Age", fontsize = 14)
ax.set_ylabel("Weight [kg]", fontsize = 14)

# 身長と体重データの散布図
ax.scatter(x, y, color = "blue")

# line()メソッドを用いて近似曲線を描画
x2 = np.linspace(0, 70, 100)
y2 = z.line(x2)
ax.plot(x2, y2, color = "red")

plt.show()

K分割交差検証 K-fold cross validation

 試行回数は NT = 8000 としています。乱数を使っているので、シード値や試行回数の設定によっては(極端に偏った基底関数が選択されるために)、稀におかしな結果が出力されることもありますが、概ね試行回数 20000 回以上で安定した結果が得られます。一般には乱数を用いない 普通のガウス基底モデル で交差検証を行なっても十分な精度が得られるので、処理速度が気になる人は、そちらも試してみてください。
 

リーブ・ワン・アウト交差検証

 K-分割交差検証において、K = N としたとき、すなわち N 個のデータを N 分割して 1 個ずつのデータをテストデータとする手法を リーブ・ワン・アウト交差検証 (LOOCV : leave-one-out cross-validation) とよびます。当サイトでは具体的に扱いませんが、要するに上のコードで K = N と書き直せば、リーブ・ワン・アウトになります。データ数が少ない時に効果を発揮する手法です。

 ≫ 非線形モデル① 非線形ガウス基底モデル