LoginSignup
13

More than 5 years have passed since last update.

Coursera Machine Learningの課題をPythonで: ex5(正則化パラメータの調整)

Last updated at Posted at 2015-11-09

CourseraのMachine Learningクラス(Andrew Ng先生)のMatlab/Octaveプログラミング課題をPythonで実装するシリーズ。コンセプトは変わらず以下のとおりです:

  • 課題のコードをそのまま再現するのではなく、scikit-learnなどのPythonライブラリを使ってできるだけ効率的に実装する

今週(Week6)は"Advice For Applying Machine Learning"という題で、新しい学習モデルを学ぶのではなく、モデルパラメータのチューニングの方法、モデルの性能の検証方法を学びます。このテーマに1週割り当てるあたりに、このコースの「理論偏重ではなく実践的」という特長が現れているのではないかと思います。

モデルのチューニング方法についての、ざくっとした内容は以下のとおり。

  • データがある場合、訓練データ、交差検定(Cross-validation)データ、テストデータに分ける。Andrew先生の推奨は、6:2:2の割合。
  • 訓練データを用いて、異なるモデルやパラメータで学習する。
  • 交差検定をしてどのモデル・パラメータがいいか決定。その際、Learning Curveを描いて決定する。
  • 最後に決定したモデルの性能をテストデータで測定する。

プログラミング課題もこの手順ですすめていきます。

まず、データの読み込み

scipyのscio.loadmat()でmatlabの.mat形式のデータを読み込めます。

import numpy as np
import matplotlib.pyplot as plt
import scipy.io as scio
from sklearn import linear_model, preprocessing

# scipy.io.loadmat()を使ってmatlabデータを読み込み
data = scio.loadmat('ex5data1.mat')
X = data['X']
Xval = data['Xval']
Xtest = data['Xtest']
y = data['y']
yval = data['yval']
ytest = data['ytest']

今回のデータは X=ダムの水位レベル を用いて、y=ダムから流出する水量 を予測するものだそうです。

まず線形回帰してみる

とりあえず線形回帰して、プロットしてみます。

model = linear_model.Ridge(alpha=0.0)
model.fit(X,y)

px = np.array(np.linspace(np.min(X),np.max(X),100)).reshape(-1,1)
py = model.predict(px)
plt.plot(px, py)
plt.scatter(X,y)
plt.show()

いつも使っているlinear_model.LinearRegression()モデルでもいいのですが、のちのち正則化項を入れるのでRidge()モデルを使っています。このモデルではパラメータのalphaで正則化の強さを指定できますが、alpha=0.0にすると正則化なしになり、LinearRegression()モデルと同じになります。

結果はこちら。

スクリーンショット 2015-11-07 09.49.10.png

見ての通り、直線ではうまく当てはまらないデータです。

それでも線形回帰でLearning Curveを描いてみる

直線では当てはまらないことは知りつつも、訓練データ数を変えて学習曲線を描いてみます。訓練データを1から12個まで変えて線形回帰を行い、訓練データに対するエラーと交差検定(Cross Validation)データに対するエラーをプロットします。「エラー」は以下の式で計算できる2乗誤差です。
$$ \frac{1}{2m} \sum (h_\theta(x^{(i)}) - y^{(i)})^2 $$
コードはこちら。

# 線形回帰でLearning Curveを描いてみる
error_train = np.zeros(11)
error_val = np.zeros(11)
model = linear_model.Ridge(alpha=0.0)
for i in range(1,12):
    # 訓練データのサブセットi個のみで回帰を実施
    model.fit( X[0:i], y[0:i] )
    # その訓練データのサブセットi個でのエラーを計算
    error_train[i-1] = sum( (y[0:i] - model.predict(X[0:i]))**2 ) / (2*i)
    # 交差検定用データでのエラーを計算
    error_val[i-1] = sum( (yval - model.predict(Xval) )**2 ) / (2*yval.size)

px = np.arange(1,12)
plt.plot(px, error_train, label="Train")
plt.plot(px, error_val, label="Cross Validation")
plt.xlabel("Number of training examples")
plt.ylabel("Error")
plt.legend()
plt.show()

結果はこうなります。
スクリーンショット 2015-11-07 17.56.21.png

訓練データを12(全部)まで増やしても、Trainデータ、Cross Validationデータともに誤差が下がりません。線形回帰モデルでは当てはまりが悪いということで、次は多項式フィッティングを試します。

多項式フィッティング

上で実施した線形回帰の仮説は
$$ h_\theta(x) = \theta_0 + \theta_1x$$
でしたが、多項式フィッティングはここに$x$の階乗の項を加えていきます。
$$ h_\theta(x) = \theta_0 + \theta_1x + \theta_2x^2 + \theta_3x^3 + ... + \theta_px^p$$
のような式です。
具体的には、特徴量$x$の階乗の数値をあらかじめ計算し、これを$x_1, x_2, x_3 ...$という新たな特徴量とし、このデータを用いて
$$ h_\theta(x) = \theta_0 + \theta_1x_1 + \theta_2x_2 + \theta_3x_3 + ... + \theta_px_p$$
で表現される線形回帰モデルで学習します。

scikit-learnでは、sklearn.preprocessing.PolynomialFeaturesという、この多項式の特徴量を計算・作成してくれるクラスがあるのでこれを利用します。
コードはこちら。

# Xの階乗を計算して新しい特徴量 X_poly とする
# Xは m x 1行列、X_polyは m x 8 行列
poly = preprocessing.PolynomialFeatures(degree=8, include_bias=False)
X_poly = poly.fit_transform(X)

# X_polyを使って線形回帰
model = linear_model.Ridge(alpha=0.0)
model.fit(X_poly,y)

# プロット
px = np.array(np.linspace(np.min(X)-10,np.max(X)+10,100)).reshape(-1,1)
# 今回のモデルはx_polyをインプットとして受け付けるので、プロット用のxも階乗の形に展開
px_poly = poly.fit_transform(px)
py = model.predict(px_poly)
plt.plot(px, py)
plt.scatter(X, y)
plt.show()

フィッティングの結果はこちら。

スクリーンショット 2015-11-07 21.25.37.png

8次多項式でフィッティングするとすべての訓練データに当てはまります。が、これは過学習で、新しいデータにはうまく予測できないモデルになっている可能性があります。今度はこのモデルを交差検定用データで検証しつつ、正則化項を入れて正則化パラメータを調整していきます。

正則化パラメータのチューニング

正則化項を入れることにより、線形回帰のコスト関数は
$$ J = \frac{1}{2m} \sum_{i=1}^m (h_\theta(x^{(i)}) - y^{(i)})^2 + \frac{\lambda}{2m} \sum_{j=1}^n \theta_j^2$$ のようになります。第2項が正則化項で、0から外れたパラメータの値にペナルティを与えることにより、過学習を防ぐことができます。この形の正則化をL2正則化、Ridge回帰、などと呼びます。

第2項の分子にある$\lambda$が正則化の強さを調整するパラメータです。上で見たように、これはlinear_model.Ridge()ではalphaパラメータに対応します。Courseraと同じくこのパラメータを 0.001, 0.003, 0.01, 0.03, 0.1, 0.3, 1, 3, 10 と変えてlearning curveをプロットし、どの$\lambda$がいいか検討します。

コードはこちら。

# Xの階乗を計算して新しい特徴量 X_poly とする
# Xは m x 1行列、X_polyは m x 8 行列
poly = preprocessing.PolynomialFeatures(degree=8, include_bias=False)
X_poly = poly.fit_transform(X) # 訓練データ
Xval_poly = poly.fit_transform(Xval) # Cross Validationデータ

# λを変えてLearning Curveを描いてみる
error_train = np.zeros(9)
error_val = np.zeros(9)
lambda_values = np.array([0.001, 0.003, 0.01, 0.03, 0.1, 0.3, 1.0, 3.0, 10.0])
for i in range(0,9):
    # X_polyを使って線形回帰
    model = linear_model.Ridge(alpha=lambda_values[i]/10, normalize=True ) # 正則化パラメータalphaを変える
    model.fit(X_poly,y)
    # 訓練データでのエラーを計算(正則化項つき)
    error_train[i] = sum( (y - model.predict(X_poly))**2 ) / (2*y.size) + sum(sum( model.coef_**2 )) * lambda_values[i]/(2*y.size)
    # 交差検定用データでのエラーを計算(正則化項つき)
    error_val[i] = sum( (yval - model.predict(Xval_poly) )**2 ) / (2*yval.size) + sum(sum( model.coef_**2 ))* lambda_values[i]/(2*yval.size)

px = lambda_values
plt.plot(px, error_train, label="Train")
plt.plot(px, error_val, label="Cross Validation")
plt.xlabel("Lambda")
plt.ylabel("Error")
plt.legend()
plt.show()

プロットはこのようになり、交差検定でのエラー値が最小となっている$\lambda=3$あたりがよいという結果になりました。

ex5.PNG

おわりに

sklearn.linear_model.Ridge()には交差検定用のsklearn.linear_model.RidgeCV()というモデルもあり、学習させると最適なalphaの数字をいっしょに計算してくれるようです。

参考文献

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
13