前回まではニューラルネットワークの複雑な実装を行っていましたが、今回はアルゴリズムではなく、モデルの評価に重点を置いた回です。自分は統計をかじっていたことがあるので、この回はなるほどと思いました。
統計でもこのへんは大事にするのですが（例えば多重共線性）、変数やモデルが有意か有意でないかのほうに関心が行って、変数やデータ量の調整は割と勘と経験でやることが多い気がします。具体的な指標としては、例えば相関行列を取って変数の取捨選択を人力やる、少しデキるアピールするならAICといった尺度でモデルを選択するなどの方法があります。しかし、統計では正則化自体あまりメジャーなトピックではないような気がする（俺ソース）ので、ここまで理論だって「オーバーフィッティングとは何か」「アンダーフィッティングとは何か」を教えてくれたのは初めてでした。

バイアス（bias）vs分散（variance）、つまりモデルのオーバーフィッティングやアンダーフィッティングを評価し（英単語の韻を踏んだものかと）、よりよい改善戦略を取ることは、アルゴリズム問わず普遍的に使えるので、アルゴリズムそのものに関する知識よりも重要ではないかと思います。いくらいいアルゴリズムを使っても、ここの調整ができないと価値が半減してしまうので。

追記 gistからソースファイルをダウンロードできます
自分で実装　https://gist.github.com/koshian2/31133a1fae0d10e857e28fc9bc556862
組み込み　https://gist.github.com/koshian2/132312dd0039ac64895898a960402701

これまでの目次

理論

バイアスvs分散の説明に感動したので、理論をまとめておきます。

アンダーフィッティング（High bias）とオーバーフィッティング（High variance）の感覚的な理解

ロジスティック回帰による分類のグラフの再掲

アンダーフィッティング（High bias）

bias（切片）の寄与が支配的になり、汎用的な分類はできるものの、分類の精度そのものが落ちる。

オーバーフィッティング（High variance）

訓練データに最適化しすぎて、他の（テスト）データに対しての精度が悪くなる。

その前に前提

データを**訓練データ（Training set）とテストデータ（Test set）**に分ける。この比率は7:3とする¹。
正則化のパラメーターλや、$x^2, x^3, \dots$などモデルの次数を調整する場合は、データを訓練データ、交差検証データ（Cross validation set）、テストデータに分ける²。この場合の比率は6:2:2とする。
モデルのパラメーター$\Theta$（回帰分析の場合は切片や係数の値）は訓練データを用いて計算する。つまり$Θ$は訓練データにおける誤差関数$J_{train}(\Theta)$を最小化するように定められる。
訓練データで求めた$\Theta$を、交差検証データに代入し求めた誤差を$J_{CV}(\Theta)$、テストデータに代入したものを$J_{test}(\Theta)$とする。
$J_{train}(\Theta)$と$J_{CV}(\Theta)$のグラフを描く。正則化や次数のパラメーターの最適化がない場合は、$J_{CV}(\Theta)$を$J_{test}(\Theta)$と読み替えてよい。

誤差関数による定義

$J_{train}(\Theta)$も$J_{CV}(\Theta)$も高い状態→アンダーフィッティング
$J_{train}(\Theta)$は低いが、$J_{CV}(\Theta)$は高い状態→オーバーフィッティング

どちらも低くなるのがより理想的なモデル。誤差関数ではなく、相関係数や（機械学習での³）F値などの精度を尺度とする場合は、低・高を逆に読み替えます。

学習曲線（Learning curve）

縦軸に誤差、横軸に訓練データのサンプル数を取りプロットしたものです。

左のように$J_{CV}(\Theta)$と$J_{train}(\Theta)$が似た値に収束しているものの、高止まりしている状態がアンダーフィッティングです。アンダーフィッティングしている場合、訓練データの数を増やすことが精度の向上にほとんど役に立ちません。

右のように、$J_{CV}(\Theta)$と$J_{train}(\Theta)$の間に大きなギャップがある状態がオーバーフィッティングです。この場合、訓練データを増やすことが精度の向上の一つの選択肢に挙げられます。

逆に言えば、訓練データの数に見合ったサイズのモデルを作りましょうということです。アンダーフィッティングしている簡単なモデルではサンプル数を活かしきれていませんし、逆にオーバーフィッティングするような複雑なモデルではサンプル数が足りていません。

検証曲線（Validation curve）

検証曲線という訳語が正しいのかわかりませんが、横軸に正則化のパラメーターλを、縦軸に誤差を取ったものをValidation curveといいます。これは適切なλを選ぶためのものです⁴。

グラフの左側は、$J_{CV}(\Theta)$が大きく$J_{train}(\Theta)$は小さいので、オーバーフィッティングとなります。グラフの右側は$J_{CV}(\Theta)$も$J_{train}(\Theta)$と同様に大きいので、アンダーフィッティングとなります。

グラフの中央、ちょうど赤い補助線を引いた、$J_{CV}(\Theta)$を最小化するλがちょうどいい正則化のパラメーターとなります。

モデルの複雑さを横軸に取る場合は、おそらくλの場合とグラフの左右が反転しますが、同様に$J_{CV}(\Theta)$を最小化する点を選べば良いのです。

オーバーフィッティング、アンダーフィッティングの改善戦略

実際にオーバーフィッティング、アンダーフィッティングが起こってしまったときにどうすればいいのでしょうか。講義を元にまとめてみました。

コスト	アンダーフィッティング（High bias）	オーバーフィッティング（High variance）
小	λの値を減らす	λの値を増やす
中	多項式の変数を加える	変数を少なくしてモデルを小さくする
大	変数を増やしてモデルを大きくする（×訓練データを増やす）	訓練データを増やす

コストは講義に出ていたものではなく、自分が追加したもので、作業の労力を表すものです。線形回帰を想定し、使うモデルによってはこれらのコストは前後するかもしれません。アンダーフィッティングで変数を増やす場合、今まで取ったデータに新しいカラムを追加することになるので（縦方向、つまり新しいデータを取ってくるのではなく、横方向に拡張するのがポイント）、コストは大きいと判断しました。

ニューラルネットワークの場合、この他に計算量の問題が出てくるので、計算量に見合った大きめのモデルを作ってあとはλで調整するのが良いとのことです。

◎sklearn.model_selectionの組み込み関数を使う

前置きが長くなってしまいましたが、学習曲線も検証曲線もどちらもscikit-learnの組み込みで書けます。自分で実装するのと書くコードの量はそこまで減らないので、ここでは要点だけにします（ワンライナーまではいきませんが、組み込み使ったほうが明瞭になります）。

学習曲線

データをロードします。今回はサイズが少ないのでハードコーディングしています。

import numpy as np
import matplotlib.pyplot as plt

from sklearn.linear_model import Ridge
from sklearn.model_selection import learning_curve, validation_curve

# 無印＝訓練データ、val＝交差検証データ、test＝テストデータ
X = np.array([-15.9368,-29.1530,36.1895,37.4922,-48.0588,-8.9415,15.3078,-34.7063,1.3892,-44.3838,7.0135,22.7627])
Xtest = np.array([-33.3180,-37.9122,-51.2069,-6.1326,21.2612,-40.3195,-14.5415,32.5598,13.3934,44.2099,-1.1427,-12.7669,34.0545,39.2235,1.9745,29.6218,-23.6696,-9.0118,-55.9406,-35.7086,9.5102])
Xval = np.array([-16.7465,-14.5775,34.5158,-47.0101,36.9751,-40.6861,-4.4720,26.5336,-42.7977,25.3741,-31.1096,27.3118,-3.2639,-1.8183,-40.7197,-50.0132,-17.4118,3.5882,7.0855,46.2824,14.6123])
y = np.array([2.1343,1.1733,34.3591,36.8380,2.8090,2.1211,14.7103,2.6142,3.7402,3.7317,7.6277,22.7524])
ytest = np.array([3.31689,5.39769,0.13043,6.19260,17.08849,0.79951,2.82479,28.62123,17.04639,55.38437,4.07937,8.27040,31.32355,39.15906,8.08728,24.11134,2.47735,6.56606,6.03809,4.69274,10.83005])
yval = np.array([4.1702e+000,4.0673e+000,3.1873e+001,1.0624e+001,3.1836e+001,4.9594e+000,4.4516e+000,2.2276e+001,-4.3874e-005,2.0504e+001,3.8583e+000,1.9365e+001,4.8838e+000,1.1097e+001,7.4617e+000,1.4769e+000,2.7192e+000,1.0927e+001,8.3487e+000,5.2782e+001,1.3357e+001])

組み込みを使うの場合は、訓練データと交差検証データを別の変数に格納する必要はなく、どこからが訓練でどこからが交差検証かを示すインデックスを与えてやればOKです。ShuffleSplitを引数に与えることもできます。ここでは交差検証のジェネレーターを自作しています。ここがちょっと独特。

# trainとCVをくっつける（組み込みでは分ける必要はない）
Xjoin = np.r_[X, Xval].reshape(-1, 1)
yjoin = np.r_[y, yval]
# 交差検証、訓練を分けるインデックスのジェネレーターを生成する関数を作る
def gen_cv(): 
    trainIndices = np.arange(len(y))
    testIndices = np.arange(len(y), len(y)+len(yval))
    yield (trainIndices, testIndices)

Xには12x1次元のデータが格納されています。つまり最もシンプルな1次元の線形回帰を想定します。sklearn.linear_modelにあるLinearRegressionでは正則化が対応できないので、**（L2の）正則化が必要な場合はRidge**を使います。Ridgeでの正則化のパラメーターはαですが、講義でのλと同じと読み替えて差し支えなさそうです。

## 通常のLearning Curve
lambda_ = 0
regr = Ridge(alpha = lambda_)
# train_sizesは割合でも指定できるし、絶対数でもOK
train_sizes, train_scores, test_scores = learning_curve(
        regr, Xjoin, yjoin, cv=gen_cv(), train_sizes=np.arange(len(y))+1)
plt.plot(train_sizes, train_scores, label="Training")
plt.plot(train_sizes, test_scores, label="Cross Validation")
plt.xlabel("Number of training examples")
plt.ylabel("R2 score")
plt.title("Learning curve for linear regression")
plt.legend()
plt.show()

これは最低限のグラフで、もう少しきれいなグラフが見たければ公式ドキュメントを参照してください。そこに載っているコードをコピペすれば信頼区間も出ます。

参考：http://scikit-learn.org/stable/auto_examples/model_selection/plot_learning_curve.html

講義では縦軸に誤差そのものの値を出していましたが、scikit-learnのデフォルトでは**精度（回帰分析の場合は決定係数）**が縦軸になります。もちろん尺度を変えることができますが、このままの設定のほうがやりやすいです。精度が高くなるほど誤差は少なくなるので、誤差のグラフを上下逆さまにしたものと考えましょう。
したがって、オーバーフィッティング＝訓練精度は高いが交差検証での精度は低い、アンダーフィッティング＝訓練精度も交差検証の精度も低いとなります。サンプル数を増やしていくほどジャストフィットしていくのがわかりますね。

オーバーフィッティングを再現するために、多項式の変数を導入します。ここは本質的ではなく、後で記すの自分で実装場合と同様なので省略します。8次のべき乗までの項を作り、訓練データをX_poly、交差検証データをX_poly_val、2つをまとめたものをX_poly_joinの変数にそれぞれ格納するものとします。

# 多項式版の学習曲線
lambda_ = 0
regr = Ridge(alpha = lambda_)
train_sizes, train_scores, test_scores = learning_curve(
        regr, X_poly_join, yjoin, cv=gen_cv(), train_sizes=np.arange(len(y))+1)
plt.plot(train_sizes, train_scores, label="Training")
plt.plot(train_sizes, test_scores, label="Cross Validation")
plt.xlabel("Number of training examples")
plt.ylabel("R2 score")
plt.title(f"Polynomial Regression Fit (lambda = {lambda_})")
plt.legend()
plt.show()

訓練データの精度が高止まりしている（完全にフィットできている）一方で、交差検証データの数が極端に少ないか、増えるほど精度が落ちるのでこれはオーバーフィッティングになります。

検証曲線

検証曲線はvalidation_curve関数を使うことで書くことができます。使い方は学習曲線の場合と同じです。λを10^-6から10^4まで20に分けてプロットしてみます。

# 多項式版の検証曲線
lambda_vec = np.logspace(-6, 4, 20)
train_scores, test_scores = validation_curve(
    regr, X_poly_join, yjoin, "alpha", lambda_vec, cv=gen_cv())
plt.semilogx(lambda_vec, train_scores, label="Training")
plt.semilogx(lambda_vec, test_scores, label="Cross Validation")
plt.xlabel("lamdbda")
plt.ylabel("R2 score")
plt.legend()
plt.show()

λが10^0～10^1あたり（真の値は3前後）で交差検証の精度が最も高くなっています。片対数グラフとしてプロットしています。オーバーフィット（λが小さい）とアンダーフィット（λが大きい）がλを媒介変数して綺麗に変化していることを確認できます。

参考：http://scikit-learn.org/stable/auto_examples/model_selection/plot_validation_curve.html

・自分で実装

今回は組み込みと自分で実装する場合さほど変わりません。回帰分析の最適化そのものは今回本質的ではないので、組み込みにまかせました。組み込みで飛ばしたところを若干解説していきます。

1.データの読み込みとプロット

import numpy as np
import matplotlib.pyplot as plt

from sklearn.linear_model import Ridge

# 無印＝訓練データ、val＝交差検証データ、test＝テストデータ
X = np.array([-15.9368,-29.1530,36.1895,37.4922,-48.0588,-8.9415,15.3078,-34.7063,1.3892,-44.3838,7.0135,22.7627])
Xtest = np.array([-33.3180,-37.9122,-51.2069,-6.1326,21.2612,-40.3195,-14.5415,32.5598,13.3934,44.2099,-1.1427,-12.7669,34.0545,39.2235,1.9745,29.6218,-23.6696,-9.0118,-55.9406,-35.7086,9.5102])
Xval = np.array([-16.7465,-14.5775,34.5158,-47.0101,36.9751,-40.6861,-4.4720,26.5336,-42.7977,25.3741,-31.1096,27.3118,-3.2639,-1.8183,-40.7197,-50.0132,-17.4118,3.5882,7.0855,46.2824,14.6123])
y = np.array([2.1343,1.1733,34.3591,36.8380,2.8090,2.1211,14.7103,2.6142,3.7402,3.7317,7.6277,22.7524])
ytest = np.array([3.31689,5.39769,0.13043,6.19260,17.08849,0.79951,2.82479,28.62123,17.04639,55.38437,4.07937,8.27040,31.32355,39.15906,8.08728,24.11134,2.47735,6.56606,6.03809,4.69274,10.83005])
yval = np.array([4.1702e+000,4.0673e+000,3.1873e+001,1.0624e+001,3.1836e+001,4.9594e+000,4.4516e+000,2.2276e+001,-4.3874e-005,2.0504e+001,3.8583e+000,1.9365e+001,4.8838e+000,1.1097e+001,7.4617e+000,1.4769e+000,2.7192e+000,1.0927e+001,8.3487e+000,5.2782e+001,1.3357e+001])

# 訓練データのプロット
plt.plot(X, y, "rx", ms=5)
plt.xlabel("Change in water level (x)")
plt.ylabel("Water flowing out of the dam (y)")

# 回帰分析
# Ridgeを使うと正規化ありの回帰分析ができる
regr = Ridge(alpha = 0)#αはλと同じ
regr.fit(X.reshape(-1, 1), y)
plt.plot(X, regr.intercept_ + X * regr.coef_, "b--", linewidth=0.5)
plt.show()

Ridge（Regression）という言葉に馴染みがなかったのですが、この関数は$\theta_{j}^2$（L2）正則化が使える回帰分析で、alpha=0とすれば正則化なしの回帰分析（LinearRegression関数）と同じになるようです。使い方はLinearRegression同じ。

プロットしているデータはダムの水位の変化と放水量をプロットしたものです。0を中心に打ち切られているので、普通の線形回帰が苦手とする問題です。こういう打ち切りデータって統計だとトービットモデルを使いましょうと言われます。

2.学習曲線

# Learning Curve
def learning_curve(X, y, Xval, yval, lambda_):
    m = X.shape[0]
    error_train = np.zeros(m)
    error_val = np.zeros(m)
    regr = Ridge(alpha = lambda_)
    for i in range(m):
        if len(X.shape) == 1:
            regr.fit(X[:(i+1)].reshape(-1, 1), y[:(i+1)])
            error_train[i] = sum((regr.intercept_ + X[:(i+1)] * regr.coef_ - y[:(i+1)]) ** 2) / 2 / (i+1)
            error_val[i] = sum((regr.intercept_ + Xval * regr.coef_ - yval) ** 2) / 2 / len(yval)
        else:
            regr.fit(X[:(i+1)], y[:(i+1)])
            error_train[i] = sum((regr.intercept_ + np.dot(X[:(i+1), :], regr.coef_) - y[:(i+1)]) ** 2) / 2 / (i+1)
            error_val[i] = sum((regr.intercept_ + np.dot(Xval, regr.coef_) - yval) ** 2) / 2 / len(yval)
    return error_train, error_val


lambda_ = 0
m = X.shape[0]
error_train, error_val = learning_curve(X, y, Xval, yval, lambda_)
plt.plot(np.arange(1, m+1), error_train, label="Train")
plt.plot(np.arange(1, m+1), error_val, label="Cross Validation")
plt.xlabel("Number of training examples")
plt.ylabel("Error")
plt.title("Learning curve for linear regression")
plt.legend()
plt.show()

Xがベクトルの場合か行列の場合かで計算が変わるので、ちょっと汚いコードになってしまいました。ベクトルの場合はreshape(-1,1)などで行列に変換するとわかりやすいコードになるかと思います。

Θを訓練データから計算して、それを交差検証データに適用して誤差を求めているだけです。組み込みの場合は縦軸が精度ですが、こちらでは誤差を取っています。精度の良さと誤差の大きさは逆の関係なので、グラフは上下逆になります。注意してください。

3.多項式の変数を追加

擬似的にオーバーフィッティングの状況を作り出すために、多項式（$x_2, x_3, \dots$）の変数を作ります。

# 変数を多項式に拡張
def map_feature(X, p):
    m = X.shape[0]
    ret = np.zeros((m, p))
    for i in range(p):
        ret[:, i] = X ** (i+1)
    # 標準化
    mu = np.mean(ret, axis=0)
    sigma = np.std(ret, axis=0, ddof=1)#不偏分散で推定
    norm = (ret - mu) / sigma
    return norm, mu, sigma, ret
def map_feature_with_normalize(X, mu, sigma):
    p = len(mu)
    norm, d1, d2, poly = map_feature(X, p) #poly以外いらないので捨てる
    return (poly - mu) / sigma

# 訓練データ、交差検証データ、テストデータを多項式化
p = 8
X_poly, mu, sigma, dummy = map_feature(X, p)
X_poly_test = map_feature_with_normalize(Xtest, mu, sigma)
X_poly_val = map_feature_with_normalize(Xval, mu, sigma)

ただし、次数単位でスケールが変わってしまうので、同じ次数（カラム）どうしで標準化しています。標準化は平均を引いて標準偏差で割るというよくあるタイプ。ただ、Numpyのstdの計算がデフォルトで不偏分散ではないので（不偏分散にするにはddof=1オプションを入れる）、そこを無視していると他の計算ソフトと数値の差が出ることがあります。

回帰（曲）線をプロットしてみると次のようになります。

def plot_fit(min_x, max_x, mu, sigma, intercept, coef):
    x = np.arange(min_x - 15, max_x + 25, 0.05)
    x_poly = map_feature_with_normalize(x, mu, sigma)
    plt.plot(x, np.dot(x_poly, coef) + intercept, "b--", linewidth=1)

# 多項式によるオーバーフィッティングをプロット
lambda_ = 0
regr = Ridge(alpha=lambda_)
regr.fit(X_poly, y)

plt.plot(X, y, "rx", ms=5)
plot_fit(min(X), max(X), mu, sigma, regr.intercept_, regr.coef_)
plt.xlabel("Change in water level (x)")
plt.ylabel("Water flowing out of the dam (y)")
plt.title(f"Polynomial Regression Fit (lambda = {lambda_})")
plt.show()

見事にオーバーフィットしてます。

4.多項式の学習曲線

# 多項式の学習曲線
error_train, error_val = learning_curve(X_poly, y, X_poly_val, yval, lambda_)
plt.plot(np.arange(1, m+1), error_train, label="Train")
plt.plot(np.arange(1, m+1), error_val, label="Cross Validation")
plt.xlabel("Number of training examples")
plt.ylabel("Error")
plt.title(f"Polynomial Regression Fit (lambda = {lambda_})")
plt.show()

1次元の多項式拡張だったのでグラフに回帰（曲）線を引けばすぐオーバーフィットはわかりますが、変数の次元が大量になったときも学習曲線を見れば一発でわかるのがとても良いです（なぜなら、交差検証データでの誤差が大きく、訓練データでの誤差が小さいどころかほぼゼロだから）。この概念統計に逆輸入されればいいのに。

5.検証曲線

# 検証曲線
def validation_curve(X, y, Xval, yval):
    lambda_vec = np.array([0,0.001,0.003,0.01,0.03,0.1,0.3,1,3,10])
    error_train = np.zeros(len(lambda_vec))
    error_val = np.zeros(len(lambda_vec))
    for i in range(len(lambda_vec)):
        regr = Ridge(alpha=lambda_vec[i])
        regr.fit(X, y)
        error_train[i] = np.sum((np.dot(X, regr.coef_) + regr.intercept_ - y) ** 2) / 2 / X.shape[0]
        error_val[i] = np.sum((np.dot(Xval, regr.coef_) + regr.intercept_ - yval) ** 2) / 2 / Xval.shape[0]
    return lambda_vec, error_train, error_val

lambda_vec, error_train, error_val = validation_curve(X_poly, y, X_poly_val, yval)

print()
print("lambda\t\tTrain Error\tValidation Error")
for l, et, ev in zip(lambda_vec, error_train, error_val):
    print(f"{l}\t{et}\t{ev}")

plt.plot(lambda_vec, error_train, label="Train")
plt.plot(lambda_vec, error_val, label="Cross Validation")
plt.xlabel("lamdbda")
plt.ylabel("Error")
plt.legend()
plt.show()

最適なλを選択するために、λの値を複数とり、正則化のパラメーターを買えながら回帰分析を行います。学習曲線でもそうでしたが、誤差をプロットする際は、正則化の項を外して構わないそうです⁵。

λ＝3のときに$J_{CV}(\Theta)$が最小になっています。λの選び方がわかりましたね。

出力

lambda          Train Error     Validation Error
0.0     0.02888975680422322     53.878324256520806
0.001   0.1126917209209884      9.839886446279376
0.003   0.17091634131481395     16.29734597988366
0.01    0.2214976601857385      16.945976952467824
0.03    0.2818487032185391      12.828901598311766
0.1     0.4593151824432713      7.587108155937583
0.3     0.9217549102562467      4.636926931005302
1.0     2.0761807140822617      4.260709395467283
3.0     4.9013458026611465      3.8229413339059226
10.0    16.092216544442497      9.945473577441716

Validation Errorの具体的な値はOctaveでやったときとだいぶ値が違いましたが、多分どこかでミスってるのではないかなと思います。でも結論が変わらないから気にしなくていいかなと（どうせ組み込み使うし）

以上です。オーバーフィッティング、アンダーフィッティングが明瞭に理解できてとても楽しかったです。

ニューラルネットワークの複雑な概念をはじめ細かい計算を「無理に理解しようとしなくてもいいよ、私もきっちりわかっていないから」とよく言っていたAndrew Ng先生が、「オーバーフィット、アンダーフィットに対する改善策がなぜそうなるのかわからなかったら、完璧に理解できるまで時間をかけて考えたほうがいい」と言っていたのがとても印象的でした。それだけ重要な概念なのでしょう。

テストデータが交差検証データの文脈で使われていることもあるので、ここらへんは曖昧っぽい ↩
これらの値がハイパーパラメータだからだと思われる。テストデータにフィットさせてしまうと、λの選択自体がオーバー／アンダーフィットしてしまうから。 ↩
統計で使われる回帰分析でのf値（F検定）というものもありますが、全く別のものです ↩
同様にモデルに含む変数の数も検証曲線と同じアプローチで最適化できそうですが、一般に検証曲線と言ったときは、λに関する最適化を行うものだそうです。 ↩
同一のλのときには、正則化の項が定数項としてしか働かないため、グラフが上下に動くだけかと思われる。Θを求めるときに機能していればよい。 ↩