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

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形式のデータを読み込めます。

```py
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=ダムから流出する水量 を予測するものだそうです。

#### まず線形回帰してみる

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

```py
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()`モデルと同じになります。

結果はこちら。

<img width="551" alt="スクリーンショット 2015-11-07 09.49.10.png" src="https://qiita-image-store.s3.amazonaws.com/0/96897/eebe313b-445b-db1e-ba4e-7e1700c53941.png">

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

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

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

```py
# 線形回帰で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()
```

結果はこうなります。
<img width="581" alt="スクリーンショット 2015-11-07 17.56.21.png" src="https://qiita-image-store.s3.amazonaws.com/0/96897/4562b72e-1bb3-31fa-23b7-9d2c6841570d.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`という、この多項式の特徴量を計算・作成してくれるクラスがあるのでこれを利用します。
コードはこちら。

```py
# 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()
```

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

<img width="561" alt="スクリーンショット 2015-11-07 21.25.37.png" src="https://qiita-image-store.s3.amazonaws.com/0/96897/b87c8bc6-bb26-7d1d-491a-be1fbea854cc.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$$ のようになります。第２項が正則化項で、0から外れたパラメータの値にペナルティを与えることにより、過学習を防ぐことができます。この形の正則化をL2正則化、Ridge回帰、などと呼びます。

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

コードはこちら。

```py
# 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](https://qiita-image-store.s3.amazonaws.com/0/96897/11705875-61e5-0494-0ceb-89c2c437c518.png)


#### おわりに

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

参考文献

- [scikit-learnの解説](http://scikit-learn.org/stable/modules/linear_model.html#ridge-regression)
- [回帰モデルにおけるL1正則化とL2正則化の効果](http://breakbee.hatenablog.jp/entry/2015/03/08/041411)