Go to Qiita Advent Calendar 2024 Top
LoginSignup
33
56

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 1 year has passed since last update.

@maskot1977(Ikemen Mas Kot)

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

Last updated at Posted at 2019-12-02

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

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

分類問題

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

# 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>)

# 表計算っぽいデータ処理用ライブラリのインポート
import pandas as pd 
# データを読み込んでデータフレーム形式として保存
df = pd.read_csv('pima-indians-diabetes.txt', delimiter="\t", index_col=0)
# 中身の確認
df
NumTimePreg OralGluTol BloodPres SkinThick SerumInsulin BMI PedigreeFunc Age Class
0 6 148 72 35 0 33.6 0.627 50 1
1 1 85 66 29 0 26.6 0.351 31 0
2 8 183 64 0 0 23.3 0.672 32 1
3 1 89 66 23 94 28.1 0.167 21 0
... ... ... ... ... ... ... ... ... ...
765 5 121 72 23 112 26.2 0.245 30 0
766 1 126 60 0 0 30.1 0.349 47 1
767 1 93 70 31 0 30.4 0.315 23 0

768 rows × 9 columns

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

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

# 目的変数
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 では、訓練データとテストデータに分割するためのメソッドを提供しています。

# 訓練データとテストデータに分割するメソッドのインポート
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 を使う方が便利です。

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

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 を用いると便利です。

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

0.006864218999993454

正解率の計算

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

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

0.7800289435600579

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

0.7402597402597403

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

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

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

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])

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

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

ROCカーブ・PRカーブ

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

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

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]])

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

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

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

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

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

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

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

%%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

評価指標を記録して比較

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

scores = []

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

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

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

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

%%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

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

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

勾配ブースティング

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

%%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

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

%%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

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

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

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')

visualize_classification_result(scores)

output_45_0.png

output_45_1.png

output_45_2.png

多層パーセプトロン

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

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

%%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

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

%%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

性能を比較します。

visualize_classification_result(scores)

output_51_0.png

output_51_1.png

output_51_2.png

課題1

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

# 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 データを分類し、性能を評価しなさい。

# 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 データセットなどのような２クラス分類ではなく、３クラス分類であるため、目的変数を次のように前処理しなければならない。

# 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
Cement Slag Fly ash Water SP Coarse Aggr. Fine Aggr. SLUMP(cm) FLOW(cm) Compressive Strength (28-day)(Mpa)
No
1 273.0 82.0 105.0 210.0 9.0 904.0 680.0 23.0 62.0 34.99
2 163.0 149.0 191.0 180.0 12.0 843.0 746.0 0.0 20.0 41.14
3 162.0 148.0 191.0 179.0 16.0 840.0 743.0 1.0 20.0 41.81
... ... ... ... ... ... ... ... ... ... ...
101 258.8 88.0 239.6 175.3 7.6 938.9 646.0 0.0 20.0 50.50
102 297.1 40.9 239.9 194.0 7.5 908.9 651.8 27.5 67.0 49.17
103 348.7 0.1 223.1 208.5 9.6 786.2 758.1 29.0 78.0 48.77

103 rows × 10 columns

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

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

回帰モデルの性能評価

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

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

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

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

%%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

評価指標を記録して比較

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

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

勾配ブースティング

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

# 訓練データとテストデータに分割するメソッドのインポート
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

%%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

%%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

output_40_1.png

output_40_2.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

%%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

%%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

output_48_1.png

output_48_2.png

課題３

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

# 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
33
56
1

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
Sign upLogin
33
56

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?