機械学習／回帰モデルの性能評価を行ってみる

#1.はじめに
　今回は、機械学習に使われる回帰モデルの性能評価を、コードを作りながら行ってみます。

#2.データセット
　使用するデータセットは、sklearnに付属しているボストン住宅価格です。

```python
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.datasets import load_boston
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression,Ridge
from sklearn.tree import DecisionTreeRegressor
from sklearn.ensemble import RandomForestRegressor,GradientBoostingRegressor
from sklearn.pipeline import Pipeline
from sklearn.metrics import r2_score
from sklearn.model_selection import cross_val_score
from sklearn.utils import shuffle

# -------- データセットの読み込み ---------
dataset = load_boston()
X = pd.DataFrame(dataset.data, columns=dataset.feature_names)
y = pd.Series(dataset.target, name='y')
print('X.shape = ', X.shape)
print(X.join(y).head())
```
![スクリーンショット 2019-11-29 10.00.34.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/209705/c0a64bef-e130-4251-b7ff-49a743750b1f.png)
データは全部で506個、特徴量は13項目、yがターゲットとなる住宅価格です。

#3.回帰モデル
　今回使用する回帰モデルは5つです。後で、使い易い様に、パイプラインの形でまとめます。ハイパーパラメータは、デフォルトです。

```python
# ---------- パイプラインの設定 -----------
pipelines = {
  '1.Linear': Pipeline([('std',StandardScaler()),
                        ('est',LinearRegression())]),
     
  '2.Ridge' : Pipeline([('std',StandardScaler()),
                        ('est',Ridge(random_state=0))]),

  '3.Tree'  : Pipeline([('std',StandardScaler()),
                        ('est',DecisionTreeRegressor(random_state=0))]),

  '4.Random': Pipeline([('std',StandardScaler()),
                        ('est',RandomForestRegressor(random_state=0, n_estimators=100))]),  
     
  '5.GBoost': Pipeline([('std',StandardScaler()),
                        ('est',GradientBoostingRegressor(random_state=0))])
}
```
**1.Linear**
　最小二乗法を使った**線形回帰モデル(Linear)**です。

**2.Ridge**
　線形回帰モデルに、L2正則化項目を追加して、過学習を抑制した**リッジ回帰モデル(Ridge)**です。

**3.Tree**
　**決定木(Decision Tree)**による回帰モデルです。

**4.Random**
　ランダムに選んだ特徴量から複数の決定木を作成し、全ての決定木の予測を平均して出力する、**ランダムフォレスト（Random Forest)**です。

**5.GBoost**
　既存のツリー群が説明しきれない情報（残差）を後続のツリーが説明しようとする形で、予測精度を高める、**勾配ブースティング(Gradient Boosting)**です。

#4.評価指標
　誤差指標は、**R2_score**を使用します。これは、予測と実測の二乗誤差Σが、実測と実測平均の二乗誤差Σに対してどれだけ小さく出来るかというものです。

　予測が実測と全て一致すれば指標は1、予測があまりにも悪いと指標はマイナスにもなり得ます。
![スクリーンショット 2019-11-29 19.04.31.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/209705/4bf4dc84-9ba3-c2e3-2856-889172e8e96b.png)

#5.ホールドアウト法
　まず、モデルの性能評価をする一般的な方法である、ホールドアウト法を行います。データを **学習データ：テストデータ＝8:2** に分割し、学習データで学習後、未知のテストデータで評価することによって、汎化性能をみるわけです。

```python
# ----------- ホールドアウト法 -----------
X_train,X_test,y_train,y_test = train_test_split(X, y, test_size=0.20, random_state=1)  

scores = {}
for pipe_name, pipeline in pipelines.items():
    pipeline.fit(X_train, y_train)
    scores[(pipe_name,'train')] = r2_score(y_train, pipeline.predict(X_train))
    scores[(pipe_name,'test')] = r2_score(y_test, pipeline.predict(X_test))
print(pd.Series(scores).unstack())
```
![スクリーンショット 2019-11-29 10.46.24.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/209705/e10b0269-06dd-2d3f-10a4-2e67e4eab0b0.png)
　未知のテストデータで評価した結果、精度が一番高い（0.924750）のは **5.GBoost**でした。

　学習データでは、精度の一番高い(1.0000) **3.Tree** は、テストデータでは(0.821282)と大きく後退しており、過学習に陥っていることが分かります。

　**2.Ridge** は **1.Linear** の改良版のはずですが、テストデータでは精度が僅かに逆転しています。これはホールドアウト法による精度測定に、ある程度バラツキがあるためで、後でもっと厳密なk-fold法で性能比較を行います。

#6.残差プロット
　モデル性能を可視化するために、残差プロットを行います。これは、横軸に**予測値**、縦軸に**予測値と実際値の差**を取り、学習データとテストデータをプロットするものです。

```python
# ------------- 残差プロット ------------
for pipe_name, est in pipelines.items():
    y_train_pred = est.predict(X_train)
    y_test_pred = est.predict(X_test)
    plt.scatter(y_train_pred, y_train_pred - y_train, c = 'blue', alpha=0.5, marker = 'o', label = 'train')
    plt.scatter(y_test_pred, y_test_pred - y_test, c = 'red', marker ='x', label= 'test' )
    plt.hlines(y = 0, xmin = 0, xmax = 50, color = 'black')
    plt.ylim([-20, 20])
    plt.xlabel('Predicted values')
    plt.ylabel('Residuals')        
    plt.title(pipe_name)
    plt.legend()
    plt.show()
```
　コードの出力としては、1.Linear 〜 5.GBoost まで5つの残差プロットされますが、ここでは代表的なものを3つだけ上げます。

![スクリーンショット 2019-11-29 11.07.29.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/209705/55e00f15-af4f-23d1-e389-a34884819009.png)
　**線形回帰モデル**です。学習データ、テストデータとも、ほぼ同じ様な残差のバラツキに成っています。

![スクリーンショット 2019-11-29 11.07.54.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/209705/6e6a879a-9de2-21b3-0fd8-b25c173335cd.png)
　**決定木**です。学習データでは完全な残差ゼロ（精度100%)な一方で、テストデータは大きく残差がバラついています。典型的な過学習です。

![スクリーンショット 2019-11-29 11.08.09.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/209705/43ac89d0-b3c8-beb4-2f05-b38ad89b89c3.png)
　**勾配ブースティング**です。学習データ、テストデータとも、残差のバラツキが抑えられています。

#7.k-hold法
　ホールドアウト法よりも厳密なモデル評価が出来るのが **k-hold法** (k分割交差検証)です。

　具体的な手順は、まずデータをk個に分割し、1つづつ順番に選んでテストデータとし、残りのk-1個を学習データとします。

　そして、学習データで学習を行い、テストデータで精度測定を行うことをk回繰り返し、得られた精度の平均を取りモデルの精度とするのが **k-hold法** です。ここでは、k=5で行います(cv = 5 で指定)。

　なお、**cross_val_score** は、**train_test_split** の様に、自動でデータをシャッフルしませんので、最初にユーティリティー **shuffle** を使って、データをシャッフルしてから処理を行います。

```python
# -------------- k-fold法 --------------
X_shuffle, y_shuffle =shuffle(X, y, random_state= 1)  # データシャッフル

scores={}
for pipe_name, est in pipelines.items():  
    cv_results = cross_val_score(est, X_shuffle, y_shuffle, cv=5, scoring='r2')    
    scores[(pipe_name,'avg')] = cv_results.mean()
    scores[(pipe_name,'score')] = np.round(cv_results,5)  # np.roundは桁調整
print(pd.Series(scores).unstack())
```
![スクリーンショット 2019-11-29 11.21.38.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/209705/bf85de9d-762a-1a91-4732-98cee11996ee.png)
　**avg**の横にある[  ]に囲まれた5つの数字**score**が、1回毎に算出した精度です。**1.Linear**を見ると、**最低0.64681〜最高0.76342**までバラついていることが分かります。この数字を平均することで、厳密なモデル評価を行うわけです。

　**2.Ridge**は、**1.Linear**より僅かですが、精度が向上していることが分かると思います。

　最終的に、ボストン住宅価格において、最も優秀なモデルは、**5.GBoost** でした。