教師あり機械学習（分類・回帰）


いくつかの説明変数 $X_n$ に対してそれに対応する目的変数 $y$ が分かっているとき、 $y = f(X)$ となる関数 $f$ を発見することを「教師あり機械学習」と言います。その中で最も簡単なのが「線形単回帰」や「線形重回帰」です。

* 回帰 ... 目的変数 $y$ が連続値のとき
* 分類 ... 目的変数 $y$ が離散値のとき

# 分類問題

例として「ピマ・インディアンの糖尿病診断」というデータを使います


```python
# URL によるリソースへのアクセスを提供するライブラリをインポートする。
import urllib.request 
# ウェブ上のリソースを指定する
url = 'https://raw.githubusercontent.com/maskot1977/ipython_notebook/master/toydata/pima-indians-diabetes.txt'
# 指定したURLからリソースをダウンロードし、名前をつける。
urllib.request.urlretrieve(url, 'pima-indians-diabetes.txt') 
```




    ('pima-indians-diabetes.txt', <http.client.HTTPMessage at 0x7fd16c201550>)




```python
# 表計算っぽいデータ処理用ライブラリのインポート
import pandas as pd 
# データを読み込んでデータフレーム形式として保存
df = pd.read_csv('pima-indians-diabetes.txt', delimiter="\t", index_col=0)
# 中身の確認
df
```




<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>NumTimePreg</th>
      <th>OralGluTol</th>
      <th>BloodPres</th>
      <th>SkinThick</th>
      <th>SerumInsulin</th>
      <th>BMI</th>
      <th>PedigreeFunc</th>
      <th>Age</th>
      <th>Class</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>6</td>
      <td>148</td>
      <td>72</td>
      <td>35</td>
      <td>0</td>
      <td>33.6</td>
      <td>0.627</td>
      <td>50</td>
      <td>1</td>
    </tr>
    <tr>
      <th>1</th>
      <td>1</td>
      <td>85</td>
      <td>66</td>
      <td>29</td>
      <td>0</td>
      <td>26.6</td>
      <td>0.351</td>
      <td>31</td>
      <td>0</td>
    </tr>
    <tr>
      <th>2</th>
      <td>8</td>
      <td>183</td>
      <td>64</td>
      <td>0</td>
      <td>0</td>
      <td>23.3</td>
      <td>0.672</td>
      <td>32</td>
      <td>1</td>
    </tr>
    <tr>
      <th>3</th>
      <td>1</td>
      <td>89</td>
      <td>66</td>
      <td>23</td>
      <td>94</td>
      <td>28.1</td>
      <td>0.167</td>
      <td>21</td>
      <td>0</td>
    </tr>
    <tr>
      <th>...</th>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
    </tr>
    <tr>
      <th>765</th>
      <td>5</td>
      <td>121</td>
      <td>72</td>
      <td>23</td>
      <td>112</td>
      <td>26.2</td>
      <td>0.245</td>
      <td>30</td>
      <td>0</td>
    </tr>
    <tr>
      <th>766</th>
      <td>1</td>
      <td>126</td>
      <td>60</td>
      <td>0</td>
      <td>0</td>
      <td>30.1</td>
      <td>0.349</td>
      <td>47</td>
      <td>1</td>
    </tr>
    <tr>
      <th>767</th>
      <td>1</td>
      <td>93</td>
      <td>70</td>
      <td>31</td>
      <td>0</td>
      <td>30.4</td>
      <td>0.315</td>
      <td>23</td>
      <td>0</td>
    </tr>
  </tbody>
</table>
<p>768 rows × 9 columns</p>



上記の「Class」が目的変数 $y$ で、１なら糖尿病、0なら糖尿病ではないという判断になります。それを予測するモデルを作りましょう。


```python
# 説明変数
X = df.iloc[:, :8]
# 最大値を1、最小値を0にするような正規化。
# axis=1 とすれば、列ではなく行単位で正規化します。
X = X.apply(lambda x: (x-x.min())/(x.max() - x.min()), axis=0)
```


```python
# 目的変数
y = df.iloc[:, 8]
```

## 訓練データとテストデータに分ける

機械学習ではその性能評価をするために、既知データを訓練データ（教師データ、教師セットともいいます）とテストデータ（テストセットともいいます）に分割します。訓練データを用いて訓練（学習）して予測モデルを構築し、その予測モデル構築に用いなかったテストデータをどのくらい正しく予測できるかということで性能評価を行ないます。そのような評価方法を「交差検定」(cross-validation)と呼びます。ここでは、

* 訓練データ（全データの６０％）

    * X_train : 訓練データの説明変数
    * y_train : 訓練データの目的変数

* テストデータ（全データの４０％）

    * X_test : テストデータの説明変数
    * y_test : テストデータの目的変数

とし、X_train と y_train の関係を学習して X_test から y_test を予測することを目指します。

Python の機械学習用ライブラリ scikit-learn では、訓練データとテストデータに分割するためのメソッドを提供しています。

* https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.train_test_split.html


```python
# 訓練データとテストデータに分割するメソッドのインポート
from sklearn.model_selection import train_test_split 
# 訓練データ・テストデータへ6:4の比でランダムに分割
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=.4) 
```

## ロジスティック回帰

ロジスティック回帰は、線形重回帰に類似していますが、目的変数 $y$ が 0 か 1 かという離散値に対応したものです。先日の講義では、ロジスティック回帰を scipy のライブラリを用いて実装しましたが、実用的には scikit-learn を使う方が便利です。

以下のサイトから、どのようなメソッドか、そしてどのようなパラメータがあるか確認しましょう。

* https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LogisticRegression.html


```python
from sklearn.linear_model import LogisticRegression # ロジスティック回帰
classifier = LogisticRegression() # 分類器の生成
classifier.fit(X_train, y_train) #学習
```




    LogisticRegression(C=1.0, class_weight=None, dual=False, fit_intercept=True,
                       intercept_scaling=1, l1_ratio=None, max_iter=100,
                       multi_class='warn', n_jobs=None, penalty='l2',
                       random_state=None, solver='warn', tol=0.0001, verbose=0,
                       warm_start=False)



## 時間の計測

学習に必要な時間を計測するために、timeit を用いると便利です。


```python
import timeit # 実行時間を計測するためのライブラリ
timeit.timeit(lambda: classifier.fit(X_train, y_train), number=1)
```




    0.006864218999993454



## 正解率の計算

「学習に用いたデータの分類精度」と「学習に用いなかったデータの分類精度」は区別しなければなりません。前者が極端に高くて後者が低いモデルは汎化性能が低く、「過学習」（overfit）したといわれます。


```python
# 正解率 (train) : 学習に用いたデータをどのくらい正しく予測できるか
classifier.score(X_train,y_train)
```




    0.7800289435600579




```python
# 正解率 (test) : 学習に用いなかったデータをどのくらい正しく予測できるか
classifier.score(X_test,y_test)
```




    0.7402597402597403



## 学習に用いなかったデータの予測と、混同行列

正解率だけではなく、具体的にどのデータがどちらに分類されるか予測を行うことができます。


```python
# 学習に用いなかったデータを予測する
y_pred = classifier.predict(X_test)
```


```python
y_pred
```




    array([1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,
           0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0,
           0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
           0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0])



正解が分かっている場合は答え合わせができます。その集計に役立つのが混同行列です。


```python
from sklearn.metrics import confusion_matrix # 混同行列を計算するメソッド
# 予測結果と、正解（本当の答え）がどのくらい合っていたかを表す混同行列
pd.DataFrame(confusion_matrix(y_pred, y_test), 
             index=['predicted 0', 'predicted 1'], columns=['real 0', 'real 1'])
```





<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>real 0</th>
      <th>real 1</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>predicted 0</th>
      <td>47</td>
      <td>18</td>
    </tr>
    <tr>
      <th>predicted 1</th>
      <td>2</td>
      <td>10</td>
    </tr>
  </tbody>
</table>



## ROCカーブ・PRカーブ

どのデータがどちらに分類されるかは、それぞれの分類への「自信の強さ」から判断されます。その自信の強い順に並べることで、ROCカーブやPRカーブなどを描いて、その予測モデルの性能を評価することができます。


```python
# 予測結果の自信の強さを計算する
y_proba = classifier.predict_proba(X_test)
```


```python
y_proba
```




    array([[0.21417054, 0.78582946],
           [0.46404957, 0.53595043],
           [0.70401466, 0.29598534],
           [0.75314361, 0.24685639],
           [0.76452966, 0.23547034],
           [0.33685542, 0.66314458],
           [0.76393323, 0.23606677],
           [0.82487552, 0.17512448],
           [0.87720401, 0.12279599],
           [0.83530283, 0.16469717],
           [0.64980016, 0.35019984],
           [0.78574888, 0.21425112],
           [0.51054138, 0.48945862],
           [0.24870259, 0.75129741],
           [0.91082684, 0.08917316],
           [0.86200773, 0.13799227],
           [0.71562431, 0.28437569],
           [0.62886446, 0.37113554],
           [0.63181921, 0.36818079],
           [0.77975231, 0.22024769],
           [0.65396517, 0.34603483],
           [0.81535938, 0.18464062],
           [0.54607196, 0.45392804],
           [0.79688063, 0.20311937],
           [0.80333846, 0.19666154],
           [0.728435  , 0.271565  ],
           [0.36817034, 0.63182966],
           [0.54025915, 0.45974085],
           [0.6614052 , 0.3385948 ],
           [0.74309548, 0.25690452],
           [0.92572332, 0.07427668],
           [0.80406998, 0.19593002],
           [0.61165474, 0.38834526],
           [0.43564389, 0.56435611],
           [0.42922327, 0.57077673],
           [0.61369072, 0.38630928],
           [0.68195508, 0.31804492],
           [0.86971152, 0.13028848],
           [0.81006182, 0.18993818],
           [0.86324924, 0.13675076],
           [0.82269894, 0.17730106],
           [0.48717372, 0.51282628],
           [0.72772261, 0.27227739],
           [0.81581007, 0.18418993],
           [0.54651378, 0.45348622],
           [0.65486361, 0.34513639],
           [0.69695761, 0.30304239],
           [0.50397912, 0.49602088],
           [0.70579261, 0.29420739],
           [0.56812519, 0.43187481],
           [0.28702944, 0.71297056],
           [0.78684682, 0.21315318],
           [0.77913962, 0.22086038],
           [0.20665217, 0.79334783],
           [0.64020202, 0.35979798],
           [0.54394942, 0.45605058],
           [0.74972094, 0.25027906],
           [0.89307226, 0.10692774],
           [0.63129007, 0.36870993],
           [0.775181  , 0.224819  ],
           [0.88651222, 0.11348778],
           [0.83087546, 0.16912454],
           [0.52015754, 0.47984246],
           [0.17895175, 0.82104825],
           [0.68620306, 0.31379694],
           [0.6503939 , 0.3496061 ],
           [0.53702941, 0.46297059],
           [0.74395419, 0.25604581],
           [0.79430285, 0.20569715],
           [0.70717315, 0.29282685],
           [0.74036824, 0.25963176],
           [0.35031104, 0.64968896],
           [0.59128595, 0.40871405],
           [0.62945511, 0.37054489],
           [0.85812094, 0.14187906],
           [0.95492842, 0.04507158],
           [0.82726693, 0.17273307]])




```python
# 図やグラフを図示するためのライブラリをインポートする。
import matplotlib.pyplot as plt
%matplotlib inline
```

ROCカーブと、その下面積であるAUCスコアを取り扱うメソッドはこちらです。


* https://scikit-learn.org/stable/modules/generated/sklearn.metrics.roc_curve.html


* https://scikit-learn.org/stable/modules/generated/sklearn.metrics.auc.html


```python
from sklearn.metrics import roc_curve
from sklearn.metrics import auc

# AUCスコアを出す
fpr, tpr, thresholds = roc_curve(y_test, y_proba[:, 1])
roc_auc = auc(fpr, tpr)
print ("AUC curve : %f" % roc_auc)

# ROC curve を描く
plt.figure(figsize=(4,4))
plt.plot(fpr, tpr, label='ROC curve (area = %0.2f)' % roc_auc)
plt.plot([0, 1], [0, 1], 'k--')
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.0])
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('ROC curve: AUC=%0.2f' % roc_auc)
plt.legend(loc="lower right")
plt.show()
```

    AUC curve : 0.756560


![output_26_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/643271bb-ac62-7d9d-70be-9ff50c284be8.png)


同様に、PRカーブを描くためのメソッドはこちら。

* https://scikit-learn.org/stable/modules/generated/sklearn.metrics.precision_recall_curve.html


```python
from sklearn.metrics import auc
from sklearn.metrics import precision_recall_curve
precision, recall, thresholds = precision_recall_curve(y_test, y_proba[:, 1])
area = auc(recall, precision)
print ("AUPR score: %0.2f" % area)

# PR curve を描く
plt.figure(figsize=(4,4))
plt.plot(recall, precision, label='Precision-Recall curve')
plt.xlabel('Recall')
plt.ylabel('Precision')
plt.ylim([0.0, 1.05])
plt.xlim([0.0, 1.0])
plt.title('Precision-Recall curve: AUPR=%0.2f' % area)
plt.legend(loc="lower left")
plt.show()
```

    AUPR score: 0.68


![output_28_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/b9a620a8-13d3-03b6-f63d-e2954f2ea6d0.png)



## グリッドサーチによるパラメータ・チューニング

機械学習のためのメソッドは多くのパラメータを必要とします。デフォルト（初期設定値）では良い予測はできません。良いパラメータを発見するための方法の一つが「グリッドサーチ」です。GridSearchCV では、学習データをさらに分割して（デフォルトでは３分割クロスバリデーション）、パラメータ候補の全組み合わせを試し、その中から平均的に優れた性能を示すパラメータを探します。


* https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.GridSearchCV.html



```python
%%time
# ベストなパラメーターを探し当てるためのグリッドサーチ
from sklearn.model_selection import GridSearchCV

# グリッドサーチを行うためのパラメーター（LogisticRegressionのパラメーター）
parameters = [
    {'solver': ['liblinear', 'saga'], 'penalty':['l1', 'l2'], 'C': [0.1, 1, 10, 100]},
    {'solver': ['newton-cg', 'sag', 'lbfgs' ], 'penalty':['l2'], 'C': [0.1, 1, 10, 100]},
]

#グリッドサーチ実行
classifier = GridSearchCV(LogisticRegression(), parameters, cv=3, n_jobs=-1)
classifier.fit(X_train, y_train)
print("Accuracy score (train): ", classifier.score(X_train, y_train))
print("Accuracy score (test): ", classifier.score(X_test, y_test))
print(classifier.best_estimator_) # ベストのパラメーターを持つ分類器
```

    Accuracy score (train):  0.784370477568741
    Accuracy score (test):  0.6883116883116883
    LogisticRegression(C=1, class_weight=None, dual=False, fit_intercept=True,
                       intercept_scaling=1, l1_ratio=None, max_iter=100,
                       multi_class='warn', n_jobs=None, penalty='l1',
                       random_state=None, solver='saga', tol=0.0001, verbose=0,
                       warm_start=False)
    CPU times: user 192 ms, sys: 25.4 ms, total: 218 ms
    Wall time: 2.52 s


## 評価指標を記録して比較

これから、様々な機械学習法の性能を比較したいと思います。そのために、評価指標を記録するための変数を用意しましょう。


```python
scores = []
```

分類モデルの性能を表す指標には様々なものがあります。主なものは次の通りです。それぞれ、意味を確認しましょう。


* https://scikit-learn.org/stable/modules/generated/sklearn.metrics.accuracy_score.html


* https://scikit-learn.org/stable/modules/generated/sklearn.metrics.precision_score.html


* https://scikit-learn.org/stable/modules/generated/sklearn.metrics.recall_score.html


* https://scikit-learn.org/stable/modules/generated/sklearn.metrics.f1_score.html


* https://scikit-learn.org/stable/modules/generated/sklearn.metrics.matthews_corrcoef.html

これらの評価指標をまとめて計算して、記録用の変数に格納する関数を以下のように作りました。異なるデータ分割に対して学習とクロスバリデーションを繰り返し、その平均的性能と標準偏差と学習時間を記録します。


```python
import timeit
from sklearn import metrics
def record_classification_scores(classifier_name, classifier, iter=5):
    records = []
    for run_id in range(iter):
        X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=.4) 
        print('Run ', run_id + 1)
        seconds = timeit.timeit(lambda: classifier.fit(X_train, y_train), number=1)
        print('    Learning Time (s):', seconds)
        y_pred = classifier.predict(X_test)
        y_proba = classifier.predict_proba(X_test)

        accuracy_score = metrics.accuracy_score(y_test, y_pred)
        precision_score = metrics.precision_score(y_test, y_pred, average='weighted')
        recall_score = metrics.recall_score(y_test, y_pred, average='weighted')
        f1_score = metrics.f1_score(y_test, y_pred, average='weighted')

        try:
            fpr, tpr, thresholds = metrics.roc_curve(y_test, y_proba[:, 1])
            roc_auc = metrics.auc(fpr, tpr)
        except:
            roc_auc = 0

        try:
            pre, rec, thresholds = metrics.precision_recall_curve(y_test, y_proba[:, 1])
            aupr = metrics.auc(rec, pre)
        except:
            aupr = 0

        try:
            mcc = metrics.matthews_corrcoef(y_test, y_pred)
        except:
            mcc = 0

        records.append([classifier_name, accuracy_score, precision_score, recall_score, 
                        f1_score, roc_auc, aupr, mcc, seconds])
    return records
```

それでは、先ほど作成した「ベストなパラメータを持つ分類器」を用いて学習し、性能指標を記録しましょう。


```python
%%time
scores += record_classification_scores('LR', classifier.best_estimator_)
```

    Run  1
        Learning Time (s): 0.004809510999990607
    Run  2
        Learning Time (s): 0.004076423000000773
    Run  3
        Learning Time (s): 0.004598837999992611
    Run  4
        Learning Time (s): 0.004291107000000238
    Run  5
        Learning Time (s): 0.003665049000005638
    CPU times: user 65.8 ms, sys: 3.33 ms, total: 69.1 ms
    Wall time: 67.6 ms


平均的性能とその標準偏差は次のようになります。


```python
df_scores = pd.DataFrame(scores, columns = ['Classifier', 'Accuracy', 'Precision', 'Recall', 
                                            'F1 score', 'ROC AUC', 'AUPR', 'MCC', 'Time'])
df_scores_mean = df_scores.iloc[:, :-1].mean()
df_scores_errors = df_scores.iloc[:, :-1].std()
df_scores_mean.plot(kind='bar', grid=True, yerr=df_scores_errors)
```




    <matplotlib.axes._subplots.AxesSubplot at 0x7fd15943ca90>



![output_38_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/1d6564f7-84b4-dc55-4df0-1de770a4cfe6.png)



## 勾配ブースティング

勾配ブースティング（Gradient Boosting）は、最近注目されている手法です。ここでは詳しくは説明しません。下記のサイトから、必要なパラメータをチェックしてみましょう。


* https://scikit-learn.org/stable/modules/ensemble.html#gradient-tree-boosting


```python
%%time
# ベストなパラメーターを探し当てるためのグリッドサーチ
from sklearn.model_selection import GridSearchCV

# 勾配ブースティング
from sklearn.ensemble import GradientBoostingClassifier

# グリッドサーチを行うためのパラメーター
parameters = [{
    'loss': ['deviance', 'exponential'],
    'learning_rate':[0.1,0.2],
    'n_estimators':[20,100,200],
    'max_depth':[3,5,7,9]
}]

#グリッドサーチ実行
classifier = GridSearchCV(GradientBoostingClassifier(), parameters, cv=3, n_jobs=-1)
classifier.fit(X_train, y_train)
print("Accuracy score (train): ", classifier.score(X_train, y_train))
print("Accuracy score (test): ", classifier.score(X_test, y_test))
print(classifier.best_estimator_) # ベストのパラメーター
```

    Accuracy score (train):  1.0
    Accuracy score (test):  0.7142857142857143
    GradientBoostingClassifier(criterion='friedman_mse', init=None,
                               learning_rate=0.2, loss='deviance', max_depth=7,
                               max_features=None, max_leaf_nodes=None,
                               min_impurity_decrease=0.0, min_impurity_split=None,
                               min_samples_leaf=1, min_samples_split=2,
                               min_weight_fraction_leaf=0.0, n_estimators=100,
                               n_iter_no_change=None, presort='auto',
                               random_state=None, subsample=1.0, tol=0.0001,
                               validation_fraction=0.1, verbose=0,
                               warm_start=False)
    CPU times: user 947 ms, sys: 25.1 ms, total: 973 ms
    Wall time: 25.2 s


得られた「ベストのパラメーターを持つ分類器」を使って学習し性能を記録しましょう。


```python
%%time
scores += record_classification_scores('GB', classifier.best_estimator_)
```

    Run  1
        Learning Time (s): 0.3291641410000068
    Run  2
        Learning Time (s): 0.31575948799999765
    Run  3
        Learning Time (s): 0.3144692120000059
    Run  4
        Learning Time (s): 0.3252903609999862
    Run  5
        Learning Time (s): 0.3103595519999942
    CPU times: user 1.64 s, sys: 7.28 ms, total: 1.65 s
    Wall time: 1.65 s


## 複数の手法の性能比較結果表示

複数の分類手法の性能比較を可視化する関数を作りました。


```python
import matplotlib.pyplot as plt
def visualize_classification_result(scores):
    df_scores = pd.DataFrame(scores, columns = ['Classifier', 'Accuracy', 'Precision', 'Recall', 
                                            'F1 score', 'ROC AUC', 'AUPR', 'MCC', 'Time'])
    df_scores_mean = df_scores.iloc[:, :-1].groupby('Classifier').mean()
    df_scores_errors = df_scores.iloc[:, :-1].groupby('Classifier').std()
    df_scores_mean.T.plot(kind='bar', grid=True, yerr=df_scores_errors.T, 
                          figsize=(12, 2), legend=False)
    plt.legend(loc = 'right', bbox_to_anchor = (0.7, 0.5, 0.5, 0.0))
    df_scores_mean.plot(kind='bar', grid=True, yerr=df_scores_errors, 
                        figsize=(12, 2), legend=False)
    plt.legend(loc = 'right', bbox_to_anchor = (0.7, 0.5, 0.5, 0.0))

    df_time_mean = df_scores.iloc[:, [0, -1]].groupby('Classifier').mean()
    df_time_errors = df_scores.iloc[:, [0, -1]].groupby('Classifier').std()
    df_time_mean.plot(kind='bar', grid=True, yerr=df_time_errors, 
                        figsize=(12, 2), legend=False)
    plt.yscale('log')
```


```python
visualize_classification_result(scores)
```

![output_45_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/aa00637d-4df3-2079-3d44-138c4ae7ec07.png)


![output_45_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/7f390d8e-cd2f-4613-6d24-7d4101f5d13c.png)


![output_45_2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/3de9c0f3-cc9c-76ec-15c7-d20c57f96c80.png)




## 多層パーセプトロン

多層パーセプトロン（Multi-Layer Perceptron）はディープラーニング（深層学習）の中で最も簡単なモデルで、scikit-learn にも実装があります。


以下のサイトから、どのようなメソッドか、そしてどのようなパラメータがあるか確認しましょう。

* https://scikit-learn.org/stable/modules/generated/sklearn.neural_network.MLPClassifier.html


```python
%%time
# ベストなパラメーターを探し当てるためのグリッドサーチ
from sklearn.model_selection import GridSearchCV

# 多層パーセプトロン
from sklearn.neural_network import MLPClassifier
# グリッドサーチを行うためのパラメーター
parameters = [{'hidden_layer_sizes': [8, (8, 8), (8, 8, 8)], 
               'solver': ['sgd', 'adam', 'lbfgs'],
                     'activation': ['logistic', 'tanh', 'relu'],
              'learning_rate_init': [0.1, 0.01, 0.001]}]
#グリッドサーチ実行
classifier = GridSearchCV(MLPClassifier(max_iter=10000, early_stopping=True), 
                          parameters, cv=3, n_jobs=-1)
classifier.fit(X_train, y_train)
print("Accuracy score (train): ", classifier.score(X_train, y_train))
print("Accuracy score (test): ", classifier.score(X_test, y_test))
print(classifier.best_estimator_) # ベストのパラメーターを持つ分類器
```

    Accuracy score (train):  0.7930535455861071
    Accuracy score (test):  0.7272727272727273
    MLPClassifier(activation='relu', alpha=0.0001, batch_size='auto', beta_1=0.9,
                  beta_2=0.999, early_stopping=True, epsilon=1e-08,
                  hidden_layer_sizes=8, learning_rate='constant',
                  learning_rate_init=0.1, max_iter=10000, momentum=0.9,
                  n_iter_no_change=10, nesterovs_momentum=True, power_t=0.5,
                  random_state=None, shuffle=True, solver='lbfgs', tol=0.0001,
                  validation_fraction=0.1, verbose=False, warm_start=False)
    CPU times: user 1.15 s, sys: 39.8 ms, total: 1.19 s
    Wall time: 2min 29s



得られた「ベストのパラメーターを持つ分類器」を使って学習し性能を記録しましょう。


```python
%%time
scores += record_classification_scores('MLP', classifier.best_estimator_)
```

    Run  1
        Learning Time (s): 0.4756240830000138
    Run  2
        Learning Time (s): 0.34581674499997916
    Run  3
        Learning Time (s): 0.15651393699999971
    Run  4
        Learning Time (s): 0.14490434999999025
    Run  5
        Learning Time (s): 0.005184319999955278
    CPU times: user 1.16 s, sys: 3.54 ms, total: 1.17 s
    Wall time: 1.17 s


性能を比較します。


```python
visualize_classification_result(scores)
```

![output_51_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/3f122347-13b5-2b29-6643-a34c125333c4.png)

![output_51_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/495ad18d-466e-99dd-33bd-d1be4d46b95c.png)



![output_51_2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/afd19fbe-dbbf-0998-2621-40b0490f2030.png)




## 課題1

scikit-learn は、機械学習の練習用データとして breast_cancer データセットを提供している。以下のようにしてデータセットを説明変数と目的変数にわけ、 GridSearchCV でパラメーターチューニングをしながら MLPClassifier で breast_cancer データを分類し、性能を評価しなさい。


```python
# https://scikit-learn.org/stable/modules/generated/sklearn.datasets.load_breast_cancer.html
from sklearn.datasets import load_breast_cancer
breast_cancer = load_breast_cancer()
X = pd.DataFrame(breast_cancer.data)
y = pd.DataFrame(breast_cancer.target.ravel())
```


## 課題2

scikit-learn は、機械学習の練習用データとして wine データセットを提供している。以下のようにしてデータセットを説明変数と目的変数にわけ、 GridSearchCV でパラメーターチューニングをしながら MLPClassifier で wine データを分類し、性能を評価しなさい。


```python
# https://scikit-learn.org/stable/modules/generated/sklearn.datasets.load_wine.html
from sklearn.datasets import load_wine
wine = data = load_wine()
X = pd.DataFrame(wine.data)
y = pd.DataFrame(wine.target)
```

ただし、 wine データセットは、 breast_cancer データセットなどのような２クラス分類ではなく、３クラス分類であるため、目的変数を次のように前処理しなければならない。


```python
# https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.OneHotEncoder.html
import numpy as np
from sklearn.preprocessing import OneHotEncoder
encoder = OneHotEncoder(categories="auto", sparse=False, dtype=np.float32)
y = pd.DataFrame(encoder.fit_transform(y))
```

また、MLPのノード数（ニューロン数）が少なすぎると分類性能が低いので、適当にノード数（ニューロン数）を増やすこと。



# 回帰問題

ここまで分類問題を解いてきましたが、次は回帰問題を解いてみましょう。題材として、コンクリートの組成と強度との関係を取り上げます。まず、コンクリートのデータを入手します。

データの取得元は Concrete Slump Test Data Set https://archive.ics.uci.edu/ml/datasets/Concrete+Slump+Test です。


```
# URL によるリソースへのアクセスを提供するライブラリをインポートする。
import urllib.request 
# ウェブ上のリソースを指定する
url = 'https://raw.githubusercontent.com/maskot1977/ipython_notebook/master/toydata/slump_test.data'
# 指定したURLからリソースをダウンロードし、名前をつける。
urllib.request.urlretrieve(url, 'slump_test.data') 
```




    ('slump_test.data', <http.client.HTTPMessage at 0x7ff02ed82518>)




```
import pandas as pd
df = pd.read_csv('slump_test.data', index_col=0)
df
```




<div>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>Cement</th>
      <th>Slag</th>
      <th>Fly ash</th>
      <th>Water</th>
      <th>SP</th>
      <th>Coarse Aggr.</th>
      <th>Fine Aggr.</th>
      <th>SLUMP(cm)</th>
      <th>FLOW(cm)</th>
      <th>Compressive Strength (28-day)(Mpa)</th>
    </tr>
    <tr>
      <th>No</th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>1</th>
      <td>273.0</td>
      <td>82.0</td>
      <td>105.0</td>
      <td>210.0</td>
      <td>9.0</td>
      <td>904.0</td>
      <td>680.0</td>
      <td>23.0</td>
      <td>62.0</td>
      <td>34.99</td>
    </tr>
    <tr>
      <th>2</th>
      <td>163.0</td>
      <td>149.0</td>
      <td>191.0</td>
      <td>180.0</td>
      <td>12.0</td>
      <td>843.0</td>
      <td>746.0</td>
      <td>0.0</td>
      <td>20.0</td>
      <td>41.14</td>
    </tr>
    <tr>
      <th>3</th>
      <td>162.0</td>
      <td>148.0</td>
      <td>191.0</td>
      <td>179.0</td>
      <td>16.0</td>
      <td>840.0</td>
      <td>743.0</td>
      <td>1.0</td>
      <td>20.0</td>
      <td>41.81</td>
    </tr>
    <tr>
      <th>...</th>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
    </tr>
    <tr>
      <th>101</th>
      <td>258.8</td>
      <td>88.0</td>
      <td>239.6</td>
      <td>175.3</td>
      <td>7.6</td>
      <td>938.9</td>
      <td>646.0</td>
      <td>0.0</td>
      <td>20.0</td>
      <td>50.50</td>
    </tr>
    <tr>
      <th>102</th>
      <td>297.1</td>
      <td>40.9</td>
      <td>239.9</td>
      <td>194.0</td>
      <td>7.5</td>
      <td>908.9</td>
      <td>651.8</td>
      <td>27.5</td>
      <td>67.0</td>
      <td>49.17</td>
    </tr>
    <tr>
      <th>103</th>
      <td>348.7</td>
      <td>0.1</td>
      <td>223.1</td>
      <td>208.5</td>
      <td>9.6</td>
      <td>786.2</td>
      <td>758.1</td>
      <td>29.0</td>
      <td>78.0</td>
      <td>48.77</td>
    </tr>
  </tbody>
</table>
<p>103 rows × 10 columns</p>
</div>



データの取得元の説明を読んでみましょう。ここでは、左７列を説明変数とし、一番右の列を目的変数と見なすことにします。


```
X = df.iloc[:, :-3].apply(lambda x: (x-x.min())/(x.max() - x.min()), axis=0)
y = df.iloc[:, -1]
```

## 訓練データとテストデータに分ける


```
# 訓練データとテストデータに分割するメソッドのインポート
from sklearn.model_selection import train_test_split 
# 訓練データ・テストデータへ6:4の比でランダムに分割
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=.4)
```

## 線形重回帰

最も簡単な回帰モデル、線形重回帰 （Multiple Linear Regression） を試してみましょう。


```
from sklearn.linear_model import LinearRegression 
regressor = LinearRegression() # 線形重回帰
regressor.fit(X_train, y_train) # 学習
```




    LinearRegression(copy_X=True, fit_intercept=True, n_jobs=None, normalize=False)



## 回帰モデルの性能評価

scikit-learn の回帰モデルでは、 .score() というメソッドは決定係数(R2値)を算出します。R2値は、

* 目的変数の本当の値と予測値の間のズレの二乗和をとる。
* それを、目的変数の本当の値とその平均値の間のズレの二乗和で割る。
* それを、1から引く。

というふうに計算します。したがって、

* 予測値のズレが小さくなるほど R2値は１に近づく。
* 予測値のズレが大きくなり、目的変数とその平均値とのズレに近づくにつれ、R2値は０に近づく。
* 予測値のズレがさらに大きくなり、目的変数とその平均値とのズレよりも大きくなると、R2値は負の値をとる。

という特徴があります。

それでは、教師セットで学習した後、テストセットで性能を確認しましょう。


```
regressor.score(X_train, y_train), regressor.score(X_test, y_test)
```




    (0.9224703183565424, 0.8177828980042425)



## 学習に用いなかったデータの予測

得られた回帰モデルに、学習に用いなかったデータを代入することで予測ができます。

目的変数とその予測値をプロットしたものを俗に y-y プロットと呼びます（が正式名称ではないようです）。プロットが対角線上にあるほど、良い回帰モデルができたと言えます。


```
import numpy as np
import sklearn.metrics as metrics

y_pred = regressor.predict(X_test) # 学習に用いなかったデータを代入

print("R2=", metrics.r2_score(y_test, y_pred))
print("RMSE=", np.sqrt(metrics.mean_squared_error(y_test, y_pred)))
print("MAE=", metrics.mean_absolute_error(y_test, y_pred))

%matplotlib inline
import matplotlib.pyplot as plt
plt.figure(figsize=(4,4))
plt.scatter(y_test, y_pred, alpha=0.2, c="blue")
plt.plot([y.min(), y.max()], [y.min(), y.max()], c="black")
plt.grid()
plt.xlabel("Real Y")
plt.ylabel("Predicted Y")
plt.show()
```

    R2= 0.8177828980042425
    RMSE= 3.0139396734524633
    MAE= 2.4622169354183447


![output_12_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/fa153f4e-0437-6605-fefb-30dd5cee92d3.png)


## 回帰モデルの性能評価

scikit-learnでは、回帰モデルの性能評価に次のような指標を計算するメソッドが用意されています。

* sklearn.metrics.r2_score
  * 先ほど説明した R2値です。大きな値（１に近い値）をとる回帰モデルほど、良いモデルと言えます。
    * https://scikit-learn.org/stable/modules/generated/sklearn.metrics.r2_score.html

* sklearn.metrics.mean_squared_error
  * 平均二乗誤差は、目的変数と予測値とのズレの二乗和です。これの平方根を取ったものが RMSE (Root Mean Square Error) です。それぞれ、分散、標準偏差と似たような式になります。小さな値（０に近い値）をとる回帰モデルほど、良いモデルと言えます。
    * https://scikit-learn.org/stable/modules/generated/sklearn.metrics.mean_squared_error.html

* sklearn.metrics.mean_absolute_error
  * 平均絶対誤差は、目的変数と予測値とのズレの絶対値の和です。平均二乗誤差は誤差を二乗するため、大きな誤差を強調する傾向があるのに比べ、平均絶対誤差はそのような傾向がありません。小さな値（０に近い値）をとる回帰モデルほど、良いモデルと言えます。 
    * https://scikit-learn.org/stable/modules/generated/sklearn.metrics.mean_absolute_error.html

# PLS (Partial Least Squares)

部分的最小二乗回帰 (PLS) は線形回帰手法の一種であり、説明変数を互いに無相関になるように線形変換した変数(潜在変数)を用いて回帰します。通常の線形重回帰と比べて、

* 説明変数を互いに無相関な潜在変数に変換して回帰するため、説明変数同士の相関が強い場合に回帰平面が一意に決まりにくくなるという共線性の問題が生じない。
* サンプル数が説明変数よりも少ない場合でも、潜在変数の数を減らすことで、全ての説明変数を考慮に入れて回帰できる。
* 潜在変数のうち重要なもののみを使用することで、予測に貢献しないノイズなどを除外できる。

という長所を持ちます。


```
# 訓練データとテストデータに分割するメソッドのインポート
from sklearn.model_selection import train_test_split 
# 訓練データ・テストデータへ6:4の比でランダムに分割
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=.4)
```


```
from sklearn.cross_decomposition import PLSRegression # PLS回帰を行うメソッド
regressor = PLSRegression() # 回帰器の生成
regressor.fit(X_train, y_train) # 学習
```




    PLSRegression(copy=True, max_iter=500, n_components=2, scale=True, tol=1e-06)




```
import numpy as np
import sklearn.metrics as metrics

y_pred = regressor.predict(X_test)
print("R2=", metrics.r2_score(y_test, y_pred))
print("RMSE=", np.sqrt(metrics.mean_squared_error(y_test, y_pred)))
print("MAE=", metrics.mean_absolute_error(y_test, y_pred))

%matplotlib inline
import matplotlib.pyplot as plt
plt.figure(figsize=(4,4))
plt.scatter(y_test, y_pred, alpha=0.2, c="blue")
plt.plot([y.min(), y.max()], [y.min(), y.max()], c="black")
plt.grid()
plt.xlabel("Real Y")
plt.ylabel("Predicted Y")
plt.show()
```

    R2= 0.7649827980366806
    RMSE= 3.523505754210429
    MAE= 2.8524226793359198


![output_17_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/49cc11fc-82e8-8399-3dd0-95696fcf3eca.png)


# グリッドサーチによるパラメーター・チューニング

パラメーター・チューニングをして、よりよいモデルを作りましょう。


```
%%time
# ベストなパラメーターを探し当てるためのグリッドサーチ
from sklearn.model_selection import GridSearchCV

# グリッドサーチを行うためのパラメーター
parameters = [
    {'n_components': [2, 3, 4, 5, 6], 'scale':[True, False], 'max_iter': [1000]},
]

#グリッドサーチ実行
regressor = GridSearchCV(PLSRegression(), parameters, cv=3, n_jobs=-1)
regressor.fit(X_train, y_train)
print("R2 (train): ", regressor.score(X_train, y_train))
print("R2 (test): ", regressor.score(X_test, y_test))
print(regressor.best_estimator_) # ベストのパラメーターを持つ回帰モデル
```

    R2 (train):  0.9225247360485105
    R2 (test):  0.8162623239997147
    PLSRegression(copy=True, max_iter=1000, n_components=3, scale=False, tol=1e-06)
    CPU times: user 114 ms, sys: 30.4 ms, total: 144 ms
    Wall time: 2.18 s



得られた回帰モデルを使って予測してみましょう。予測性能は改善されたでしょうか。


```
import numpy as np
import sklearn.metrics as metrics

y_pred = regressor.predict(X_test)
print("R2=", metrics.r2_score(y_test, y_pred))
print("RMSE=", np.sqrt(metrics.mean_squared_error(y_test, y_pred)))
print("MAE=", metrics.mean_absolute_error(y_test, y_pred))

%matplotlib inline
import matplotlib.pyplot as plt
plt.figure(figsize=(4,4))
plt.scatter(y_test, y_pred, alpha=0.2, c="blue")
plt.plot([y.min(), y.max()], [y.min(), y.max()], c="black")
plt.grid()
plt.xlabel("Real Y")
plt.ylabel("Predicted Y")
plt.show()
```

    R2= 0.8162623239997147
    RMSE= 3.115474987155132
    MAE= 2.3005236909984426


![output_21_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/197c48a1-ec2c-5952-93c9-a13aa948769e.png)


# 評価指標を記録して比較

これから、様々な回帰モデルの性能を比較したいと思います。そのために、評価指標を記録するための変数を用意しましょう。


```
scores = []
```

回帰モデルの評価指標をまとめて計算して、記録用の変数に格納する関数を以下のように作りました。異なるデータ分割に対して学習とクロスバリデーションを繰り返し、その平均的性能と標準偏差と学習時間を記録します。


```
import timeit
import numpy as np
from sklearn import metrics

def record_regression_scores(regressor_name, regressor, iter=5):
    records = []
    run_id = 0
    successful = 0
    max_trial = 100
    while successful < iter:
        run_id += 1
        if run_id >= max_trial:
            break

        X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=.4) 
        print('Run ', run_id)
        seconds = timeit.timeit(lambda: regressor.fit(X_train, y_train), number=1)
        print('    Learning Time (s):', seconds)
        y_pred = regressor.predict(X_test)
        r2_score = metrics.r2_score(y_test, y_pred)
        rmse_score = np.sqrt(metrics.mean_squared_error(y_test, y_pred))
        mae_score = metrics.mean_absolute_error(y_test, y_pred)

        if r2_score < 0:
            print("\t\t encountered negative r2_score")
            continue
        else:
            successful += 1

        records.append([regressor_name, r2_score, rmse_score, mae_score, seconds])
    return records
```

それでは、先ほど作成した「ベストなパラメータを持つ回帰モデル」を用いて学習し、性能評価指標を記録しましょう。


```
%%time
scores += record_regression_scores("PLS", regressor)
```

    Run  1
        Learning Time (s): 2.0181297670001186
    Run  2
        Learning Time (s): 1.9526900320001914
    Run  3
        Learning Time (s): 1.9921050099997046
    Run  4
        Learning Time (s): 2.0573012720001316
    Run  5
        Learning Time (s): 1.979584856999736
    CPU times: user 552 ms, sys: 101 ms, total: 653 ms
    Wall time: 10 s



平均的性能とその標準偏差は次のようになります。


```
df_scores = pd.DataFrame(scores, columns = ['Regressor', 'R2', 'RMSE', 'Mae', 'Time'])
df_scores_mean = df_scores.iloc[:, :-1].mean()
df_scores_errors = df_scores.iloc[:, :-1].std()
df_scores_mean.plot(kind='bar', grid=True, yerr=df_scores_errors)
```




    <matplotlib.axes._subplots.AxesSubplot at 0x7ff0177f0550>



![output_29_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/6146449e-ce20-70a3-4e27-631821888e8d.png)


# 勾配ブースティング

先ほどは勾配ブースティングによる分類をしましたが、勾配ブースティングによる回帰もできます。


```
# 訓練データとテストデータに分割するメソッドのインポート
from sklearn.model_selection import train_test_split 
# 訓練データ・テストデータへ6:4の比でランダムに分割
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=.4)
```


```
from sklearn.ensemble import GradientBoostingRegressor
regressor = GradientBoostingRegressor() # 勾配ブースティング
regressor.fit(X_train, y_train) # 学習
```




    GradientBoostingRegressor(alpha=0.9, criterion='friedman_mse', init=None,
                              learning_rate=0.1, loss='ls', max_depth=3,
                              max_features=None, max_leaf_nodes=None,
                              min_impurity_decrease=0.0, min_impurity_split=None,
                              min_samples_leaf=1, min_samples_split=2,
                              min_weight_fraction_leaf=0.0, n_estimators=100,
                              n_iter_no_change=None, presort='auto',
                              random_state=None, subsample=1.0, tol=0.0001,
                              validation_fraction=0.1, verbose=0, warm_start=False)




```
regressor.score(X_train, y_train), regressor.score(X_test, y_test)
```




    (0.9996743754326906, 0.7386973055974495)




```
import numpy as np
import sklearn.metrics as metrics

y_pred = regressor.predict(X_test)

print("R2=", metrics.r2_score(y_test, y_pred))
print("RMSE=", np.sqrt(metrics.mean_squared_error(y_test, y_pred)))
print("MAE=", metrics.mean_absolute_error(y_test, y_pred))

%matplotlib inline
import matplotlib.pyplot as plt
plt.figure(figsize=(4,4))
plt.scatter(y_test, y_pred, alpha=0.2, c="blue")
plt.plot([y.min(), y.max()], [y.min(), y.max()], c="black")
plt.grid()
plt.xlabel("Real Y")
plt.ylabel("Predicted Y")
plt.show()
```

    R2= 0.7386973055974495
    RMSE= 4.012901982806575
    MAE= 3.0486670616108


![output_34_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/662178ac-4848-60a3-c0f2-fe71f91d31da.png)



```
%%time
# ベストなパラメーターを探し当てるためのグリッドサーチ
from sklearn.model_selection import GridSearchCV

# 勾配ブースティング
from sklearn.ensemble import GradientBoostingRegressor

# グリッドサーチを行うためのパラメーター
parameters = [{
    'learning_rate':[0.1,0.2],
    'n_estimators':[20,100],
    'max_depth':[3,5]
}]

#グリッドサーチ実行
regressor = GridSearchCV(GradientBoostingRegressor(), parameters, cv=3, n_jobs=-1)
regressor.fit(X_train, y_train)
print("R2 (train): ", regressor.score(X_train, y_train))
print("R2 (test): ", regressor.score(X_test, y_test))
print(regressor.best_estimator_) # ベストのパラメーターを持つ回帰モデル
```

    R2 (train):  0.9996743754326906
    R2 (test):  0.7195388936429337
    GradientBoostingRegressor(alpha=0.9, criterion='friedman_mse', init=None,
                              learning_rate=0.1, loss='ls', max_depth=3,
                              max_features=None, max_leaf_nodes=None,
                              min_impurity_decrease=0.0, min_impurity_split=None,
                              min_samples_leaf=1, min_samples_split=2,
                              min_weight_fraction_leaf=0.0, n_estimators=100,
                              n_iter_no_change=None, presort='auto',
                              random_state=None, subsample=1.0, tol=0.0001,
                              validation_fraction=0.1, verbose=0, warm_start=False)
    CPU times: user 130 ms, sys: 15.4 ms, total: 145 ms
    Wall time: 2.26 s






```
import numpy as np
import sklearn.metrics as metrics

y_pred = regressor.predict(X_test)
print("R2=", metrics.r2_score(y_test, y_pred))
print("RMSE=", np.sqrt(metrics.mean_squared_error(y_test, y_pred)))
print("MAE=", metrics.mean_absolute_error(y_test, y_pred))

%matplotlib inline
import matplotlib.pyplot as plt
plt.figure(figsize=(4,4))
plt.scatter(y_test, y_pred, alpha=0.2, c="blue")
plt.plot([y.min(), y.max()], [y.min(), y.max()], c="black")
plt.grid()
plt.xlabel("Real Y")
plt.ylabel("Predicted Y")
plt.show()
```

    R2= 0.7195388936429337
    RMSE= 4.15741070004397
    MAE= 3.0704656592653317


![output_36_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/1fec43bf-a368-efb4-508e-95407994d059.png)



```
%%time
scores += record_regression_scores("GB", regressor.best_estimator_)
```

    Run  1
        Learning Time (s): 0.027196151000225655
    Run  2
        Learning Time (s): 0.01961764999987281
    Run  3
        Learning Time (s): 0.01894888400011041
    Run  4
        Learning Time (s): 0.019140249999964
    Run  5
        Learning Time (s): 0.020592135999777383
    CPU times: user 123 ms, sys: 2.75 ms, total: 126 ms
    Wall time: 128 ms


## 複数の手法の性能比較結果表示

複数の回帰手法の性能比較を可視化する関数を作りました。


```
def visualize_regression_result(scores):
    df_scores = pd.DataFrame(scores, columns =['Regressor', 'R2', 'RMSE', 'MAE', 'Time'])
    df_scores_mean = df_scores.iloc[:, :-1].groupby('Regressor').mean()
    df_scores_errors = df_scores.iloc[:, :-1].groupby('Regressor').std()
    df_scores_mean.T.plot(kind='bar', grid=True, yerr=df_scores_errors.T, 
                          figsize=(12, 2), legend=False)
    #plt.yscale('log')

    plt.legend(loc = 'right', bbox_to_anchor = (0.7, 0.5, 0.5, 0.0))
    df_scores_mean.plot(kind='bar', grid=True, yerr=df_scores_errors, 
                        figsize=(12, 2), legend=False)
    #plt.yscale('log')

    plt.legend(loc = 'right', bbox_to_anchor = (0.7, 0.5, 0.5, 0.0))
    df_time_mean = df_scores.iloc[:, [0, -1]].groupby('Regressor').mean()
    df_time_errors = df_scores.iloc[:, [0, -1]].groupby('Regressor').std()
    df_time_mean.plot(kind='bar', grid=True, yerr=df_time_errors, 
                        figsize=(12, 2), legend=False)
    plt.yscale('log')
```


```
visualize_regression_result(scores)
```

![output_40_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/9383d869-7308-d55b-654f-e8789a58a37d.png)


![output_40_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/e94aabe3-471b-3c94-ce06-b54edffd29f0.png)

![output_40_2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/153a0cf4-ec2c-a3d4-7151-99d0ae9b9dec.png)




# 多層パーセプトロン

さきほどは多層パーセプトロンによる分類をしましたが、多層パーセプトロンによる回帰もできます。


```
# 訓練データとテストデータに分割するメソッドのインポート
from sklearn.model_selection import train_test_split 
# 訓練データ・テストデータへ6:4の比でランダムに分割
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=.4)
```


```
from sklearn.neural_network import MLPRegressor
regressor = MLPRegressor() # 回帰器の生成
regressor.fit(X_train, y_train) # 学習
```

    /usr/local/lib/python3.6/dist-packages/sklearn/neural_network/multilayer_perceptron.py:566: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (200) reached and the optimization hasn't converged yet.
      % self.max_iter, ConvergenceWarning)





    MLPRegressor(activation='relu', alpha=0.0001, batch_size='auto', beta_1=0.9,
                 beta_2=0.999, early_stopping=False, epsilon=1e-08,
                 hidden_layer_sizes=(100,), learning_rate='constant',
                 learning_rate_init=0.001, max_iter=200, momentum=0.9,
                 n_iter_no_change=10, nesterovs_momentum=True, power_t=0.5,
                 random_state=None, shuffle=True, solver='adam', tol=0.0001,
                 validation_fraction=0.1, verbose=False, warm_start=False)




```
import numpy as np
import sklearn.metrics as metrics

y_pred = regressor.predict(X_test)
print("R2=", metrics.r2_score(y_test, y_pred))
print("RMSE=", np.sqrt(metrics.mean_squared_error(y_test, y_pred)))
print("MAE=", metrics.mean_absolute_error(y_test, y_pred))

%matplotlib inline
import matplotlib.pyplot as plt
plt.figure(figsize=(4,4))
plt.scatter(y_test, y_pred, alpha=0.2, c="blue")
plt.plot([y.min(), y.max()], [y.min(), y.max()], c="black")
plt.grid()
plt.xlabel("Real Y")
plt.ylabel("Predicted Y")
plt.show()
```

    R2= -4.0467411800805415
    RMSE= 19.21649631146132
    MAE= 17.449687389239205


![output_44_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/16217684-da86-ae22-63a8-1050a06f2651.png)



```
%%time
# ベストなパラメーターを探し当てるためのグリッドサーチ
from sklearn.model_selection import GridSearchCV

# グリッドサーチを行うためのパラメーター
parameters = [{
    'hidden_layer_sizes': [10, (10, 10)],
    'solver': ['sgd', 'adam', 'lbfgs'],
    #'solver': ['lbfgs'],
    #'activation': ['logistic', 'tanh', 'relu']
    'activation': ['relu']
}]

#グリッドサーチ実行
regressor = GridSearchCV(MLPRegressor(max_iter=10000, early_stopping=True), 
                         parameters, cv=3, n_jobs=-1)
regressor.fit(X_train, y_train)
print("R2 (train): ", regressor.score(X_train, y_train))
print("R2 (test): ", regressor.score(X_test, y_test))
print(regressor.best_estimator_) # ベストのパラメーターを持つ回帰モデル
```

    R2 (train):  0.9742637037080083
    R2 (test):  0.9562295568855493
    MLPRegressor(activation='relu', alpha=0.0001, batch_size='auto', beta_1=0.9,
                 beta_2=0.999, early_stopping=True, epsilon=1e-08,
                 hidden_layer_sizes=(10, 10), learning_rate='constant',
                 learning_rate_init=0.001, max_iter=10000, momentum=0.9,
                 n_iter_no_change=10, nesterovs_momentum=True, power_t=0.5,
                 random_state=None, shuffle=True, solver='lbfgs', tol=0.0001,
                 validation_fraction=0.1, verbose=False, warm_start=False)
    CPU times: user 222 ms, sys: 17.5 ms, total: 239 ms
    Wall time: 7.86 s





```
import numpy as np
import sklearn.metrics as metrics

y_pred = regressor.predict(X_test)
print("R2=", metrics.r2_score(y_test, y_pred))
print("RMSE=", np.sqrt(metrics.mean_squared_error(y_test, y_pred)))
print("MAE=", metrics.mean_absolute_error(y_test, y_pred))

%matplotlib inline
import matplotlib.pyplot as plt
plt.figure(figsize=(4,4))
plt.scatter(y_test, y_pred, alpha=0.2, c="blue")
plt.plot([y.min(), y.max()], [y.min(), y.max()], c="black")
plt.grid()
plt.xlabel("Real Y")
plt.ylabel("Predicted Y")
plt.show()
```

    R2= 0.9562295568855493
    RMSE= 1.789613149534058
    MAE= 1.3873465536350154


![output_46_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/f8b11316-c739-24b6-2fd1-0d903be3960e.png)



```
%%time
scores += record_regression_scores("MLP", regressor.best_estimator_)
```

    Run  1
        Learning Time (s): 0.06779548599979535
    Run  2
        Learning Time (s): 0.1298420270004499
    Run  3
        Learning Time (s): 0.1824235089998183
    Run  4
        Learning Time (s): 0.43246253200004503
    Run  5
        Learning Time (s): 0.22879209799975797
    CPU times: user 1.06 s, sys: 3.13 ms, total: 1.06 s
    Wall time: 1.07 s



```
visualize_regression_result(scores)
```

![output_48_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/f67aa90f-2012-7756-9822-81a8960962d7.png)

![output_48_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/276608e4-2a91-2329-f28f-baa1d8943260.png)

![output_48_2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/3a5db65a-d16e-ec89-4833-c96593dd3657.png)




## 課題３

scikit-learn は、機械学習の練習用データとして diabetes データセットを提供している。以下のようにしてデータセットを説明変数と目的変数にわけ、 GridSearchCV でパラメーターチューニングをしながら MLPRegressor や GradientBoostingRegressor で diabetes データを回帰し、性能を比較しなさい。


```python
# https://scikit-learn.org/stable/modules/generated/sklearn.datasets.load_diabetes.html
from sklearn.datasets import load_diabetes
diabetes = load_diabetes()
X = diabetes.data
y = diabetes.target
```