PR曲線とROC曲線を理解する(後編)
こんばんは🐼
今回は、3週間ほど前に書いた「PR曲線とROC曲線を理解する(前編)」という記事の続編です。
前回はPR曲線・ROC曲線について理解するために欠かせない混同行列や、適合率・再現率といった指標について、その意味と共になぜ必要なのかを説明しました。
今回の目的としては、
- PR曲線・ROC曲線とは何か
- ビジネス上どのように役立てることができるのか
ということについて、説明・考察していきたいと思います。
PR曲線・ROC曲線を描く目的
PR曲線・ROC曲線について説明する前に、前回の記事のおさらいも含めてPR曲線(ROC曲線)をなぜ描くのかについてお話します。
それは、「作成した機械学習をビジネス上どのように役立てるか」というところに繋がってきます。
- スパムメールを除去するタスクにおいて、必要なメールを削除してしまうよりも、受信トレイにスパムメールが残ってしまう方がマシ。
- ガンの陽性を発見するタスクにおいて、陽性の患者を陰性と診断してしまうよりも、陰性の患者を陽性と診断してしまう方がマシ。
といったように、現実では同じ予測の誤りでも、その誤りによって生じるコストやリスクが異なるため、精度(accuracy)だけではモデルの良し悪しの判断が付かないケースは往々にしてあります。
この時、「ビジネス上どちらの誤りをどの程度優先して減らすか」ということを含めてモデル評価するために、PR曲線・ROC曲線を描画します。
どちらの誤りも減らすということは根本的なモデルの性能を向上させる他ないので、基本的には片方の誤りを減らすと片方の誤りは増えるというトレードオフの関係にあります。
サンプルデータセット
本記事では、不均衡データに対するモデルの性能評価について簡潔な例を示すために、手頃に手に入るskicit-learnのトイデータセットのラベルを不均衡にして使用しようと思います。
import numpy as np from sklearn.datasets import load_digits from sklearn.model_selection import train_test_split # digitsデータセットの読み込み digits = load_digits() y = digits.target == 8 X_train, X_test, y_train, y_test = train_test_split( digits.data, y, random_state=42, stratify=y)
数字が「8である(1:陽性)」「8でない(0:陰性)」の2種類のラベルを付けることで、不均衡データに対する2値分類タスクの例として使います。
np.bincount(y)
全サンプルのうち陽性クラスが174サンプル(9.7%)、陰性クラスが1623サンプル(90.3%)であり、不均衡データとなっていることがわかります。
モデリング
モデルの性能を比較するために、様々なモデルで訓練を行い、テストデータに対する認識精度(accuracy)を出力します。
ダミー分類器
はじめに問答無用で多数派クラスと予測するダミー分類器を作ります。
from sklearn.dummy import DummyClassifier # strategy = 'most_frequent'で常に多数派クラスを予測 dummy_clf = DummyClassifier(strategy='most_frequent').fit(X_train, y_train) pred_dummy = dummy_clf.predict(X_test) print('test accuracy:', dummy_clf.score(X_test, y_test))
認識精度は90.2%となります。内容を知らない状態で精度が90.2%ということだけを聞くと「性能がいいモデルだ」と感じる人が多いかもしれませんが、実情は多数派クラスを出力しているだけです。多数派クラスが99.9%のデータであれば認識精度は99.9%となります。
この簡単な例だけでも認識精度のみを使ってモデルを評価することの危険さが解ると思います。
ロジスティック回帰
from sklearn.linear_model import LogisticRegression param_grid = [{'C': [0.001, 0.01, 0.1, 1, 10]}] grid_search = GridSearchCV(LogisticRegression(random_state=42, solver='lbfgs', max_iter=10000), param_grid, cv=5) grid_search.fit(X_train, y_train) print(f"Best parameters: {grid_search.best_params_}") print(f"Best validation score {grid_search.best_score_}") logreg_clf = grid_search.best_estimator_ pred_logreg = logreg_clf.predict(X_test) print('test accuracy:', logreg_clf.score(X_test, y_test))
グリッドサーチという手法を使って、適切なハイパーパラメータを探索していますがここではPR曲線・ROC曲線の理解が主テーマなので、パラメータの探索は最低限に行っています。(以降のモデルについても同様) 正則化パラメータC: 0.01のロジスティック回帰分類器モデルはテストデータに対して96.4%の認識精度が得られるようです。
決定木
from sklearn.tree import DecisionTreeClassifier param_grid = [{'max_depth': [1, 3, 5, 7, 9]}] grid_search = GridSearchCV(DecisionTreeClassifier(random_state=42), param_grid, cv=5) grid_search.fit(X_train, y_train) print(f"Best parameters: {grid_search.best_params_}") print(f"Best validation score {grid_search.best_score_}") tree_clf = grid_search.best_estimator_ pred_tree = tree_clf.predict(X_test) print('test accuracy:', tree_clf.score(X_test, y_test))
決定木ではロジスティック回帰よりも認識精度が低く94%となります。しかしここでも注意が必要なのは、常に多数派クラスを出力するダミー分類器ですら90.2%の認識精度を出せるのだから、これだけではまだ「ロジスティック回帰よりも決定木が悪いモデルである」と結論付けることは出来ないということです。
不均衡データでは一般に多数派クラスに寄せて予測すると精度が出やすいので、少数派クラスをどの程度正しく予測できているかどうかは精度からは詳しく判定できないことに注意が必要です。
ランダムフォレスト
from sklearn.ensemble import RandomForestClassifier forest_clf = RandomForestClassifier(n_estimators=400, random_state=42) forest_clf.fit(X_train, y_train) pred_forest = forest_clf.predict(X_test) print('test accuracy:', forest_clf.score(X_test, y_test))
ランダムフォレストの認識精度は96.9%となりました。続いて2値分類タスクのモデルを正しく評価するための有効手段である混同行列と各種性能指標を表示します。
混同行列・各種性能指標の表示
PR曲線は前回の記事で説明した適合率・再現率のトレードオフの関係を、同じくROC曲線は後述する真陽性率と偽陽性率の関係を視覚的に表したものです。
そのためまずは可視化のもととなる混同行列と各種性能指標の表示を行います。これらはsklearn.metricsモジュールのconfusion_matrix関数および、classification_report関数で実現できます。
ダミー分類器
from sklearn.metrics import confusion_matrix, classification_report import matplotlib.pyplot as plt import seaborn as sns targets = ['not 8', '8'] matrix = confusion_matrix(y_test, pred_dummy) # 混同行列の可視化 sns.heatmap(matrix, annot=True, cmap='Blues', cbar=False, xticklabels=targets, yticklabels=targets) # 各種性能指標の表示 print(classification_report(y_test, pred_dummy, target_names=targets))
当然ながら全てを「8以外(not 8)」と予測したので、「8である(8)」と予測したサンプルの数は0になっています。
classification_reportでは「8以外」を陽性としたときの適合率・再現率・F値と、「8」を陽性としたときの適合率・再現率・F値の双方が出力されます。
これにより、少数派クラスについて全く予測ができていない(再現率がゼロである)ことが分かります。現実では少数派クラスが見つけたいもの(陽性)であることがほとんどです。
ロジスティック回帰
他のモデルについても同様に、混同行列と各種指標を表示します。
matrix = confusion_matrix(y_test, pred_logreg) sns.heatmap(matrix, annot=True, cmap='Blues', cbar=False, xticklabels=targets, yticklabels=targets) print(classification_report(y_test, pred_logreg, target_names=targets))
ロジスティック回帰では、「8」を陽性としたときの再現率は0.68まで増加しており、68%の「8」を正しく予測できていることがわかります。また、適合率が0.94と高いことから「8」であると予測したものが、実際には「8以外」であることは少ないといえます。
十分な精度が得られているかは別で検討が必要ですが、幾分意味のあるモデルとなっているといえるでしょう。
決定木
matrix = confusion_matrix(y_test, pred_tree) sns.heatmap(matrix, annot=True, cmap='Blues', cbar=False, xticklabels=targets, yticklabels=targets) print(classification_report(y_test, pred_tree, target_names=targets))
決定木の「8」を陽性としたときの再現率はロジスティック回帰とほぼ同等ですが、適合率は大きく下がっており、全体としてロジスティック回帰より性能が悪いことが分かります。
ランダムフォレスト
matrix = confusion_matrix(y_test, pred_forest) sns.heatmap(matrix, annot=True, cmap='Blues', cbar=False, xticklabels=targets, yticklabels=targets) print(classification_report(y_test, pred_forest, target_names=targets))
ランダムフォレスト分の性能指標はロジスティック回帰とよく似ています。「8」を陽性としたときの適合率が100%であることから、より曖昧な「8」についても「8」と予測することで、適合率と再現率のトレードオフを調整することができます。(そのような用途でPR曲線が便利だと考えられます)
PR曲線・ROC曲線の描画
いよいよPR曲線・ROC曲線の描画です。PR曲線・ROC曲線の解釈については、前編の記事で説明した指標(classification report関数で出力される指標)についての理解が出来ていれば難しいことはありません。
探しているクラス(今回の例でいえば「8」)を陽性とした上で、
- PR曲線 : X軸を再現率(recall)、Y軸を適合率(precision)としてあらゆる閾値についてプロットしたグラフ
- ROC曲線 : X軸を真陽性率(=再現率)、Y軸を偽陽性率としてあらゆる閾値についてプロットしたグラフ
この「あらゆる閾値について」という所がポイントで、各種の分類アルゴリズムは分類の際に、サンプルが陽性に属する確率というものを出力します。
この確率は各分類器のpredict_probaメソッドで出力できます。(但し一部の分類器はpredict_probaメソッドを持たず、decision_functionメソッドで予測の不確実性を評価します)
デフォルトでは陽性の確率が50%以上のサンプルを陽性と予測しますが、この閾値は変更可能です。
閾値を増やすと、より陽性の確率が高いもののみを陽性と判定するため、適合率が高まり、再現率が下がる傾向になります。逆に閾値を減らすと、陽性かどうか曖昧なものまで陽性と判定するため、適合率が下がり、再現率が高まる傾向になります。
PR曲線はこの閾値の変更による適合率と再現率のトレードオフの関係を可視化したものです。
ROC曲線についてもトレードオフの関係を表しています。真陽性率は再現率と同じですが、偽陽性率については「本当は陰性のサンプルの中で陽性と予測された割合」を表しています。
偽陽性率は誤りの割合なので低いほど良い指標です。初見ではやや混乱するので定期的に確認して覚えていきましょう。
偽陽性率は、1 - (陰性を陽性と見た時の再現率)と言い換えることも出来ます。(陰性を陽性と見た時の再現率、つまり本当に陰性のサンプルの中で陰性と予測された割合は「特異度」と呼びます)
そのため自分はROC曲線が何を表しているか混乱した時は、陽性の再現率と陰性の再現率のトレードオフと覚えています。(厳密に言うと再現率は「実際に陽性のサンプルのうち何%を正しく予測できたのか」なので、再現率と特異度のトレードオフですが、特異度という言葉はあまり使われないので)
長くなりましたが、ロジスティック回帰とランダムフォレストについてそれぞれの曲線を描画してみましょう。
sns.set(); from sklearn.metrics import precision_recall_curve, roc_curve fig, axes = plt.subplots(1, 2, figsize=(12, 6)) # PR曲線の描画 (ロジスティック回帰, ランダムフォレスト) precision_logreg, recall_logreg, thresholds_logreg = precision_recall_curve(y_test, logreg_clf.predict_proba(X_test)[:, 1]) precision_forest, recall_forest, thresholds_forest = precision_recall_curve(y_test, forest_clf.predict_proba(X_test)[:, 1]) default_threshold_logreg = np.argmin(np.abs(thresholds_logreg - 0.5)) default_threshold_forest = np.argmin(np.abs(thresholds_forest - 0.5)) axes[0].plot(recall_logreg, precision_logreg, label='LogisticRegression') axes[0].plot(recall_forest, precision_forest, label='RandomForestClassifier', linestyle='--') axes[0].plot(recall_logreg[default_threshold_logreg], precision_logreg[default_threshold_logreg], 'x', c='k', markersize=10, label='default threshold(50%)') axes[0].plot(recall_forest[default_threshold_forest], precision_forest[default_threshold_forest], 'x', c='k', markersize=10) axes[0] axes[0].legend() axes[0].set_title('PR Curve', fontsize=14) axes[0].set_xlabel('recall') axes[0].set_ylabel('precision') # ROC曲線の描画 (ロジスティック回帰, ランダムフォレスト) fpr_logreg, tpr_logreg, thresholds_logreg = roc_curve(y_test, logreg_clf.predict_proba(X_test)[:, 1]) fpr_forest, tpr_forest, thresholds_forest = roc_curve(y_test, forest_clf.predict_proba(X_test)[:, 1]) default_threshold_logreg = np.argmin(np.abs(thresholds_logreg - 0.5)) default_threshold_forest = np.argmin(np.abs(thresholds_forest - 0.5)) axes[1].plot(fpr_logreg, tpr_logreg, label='LogisticRegression') axes[1].plot(fpr_forest, tpr_forest, label='RandomForestClassifier', linestyle='--') axes[1].plot(fpr_logreg[default_threshold_logreg], tpr_logreg[default_threshold_logreg], 'x', c='k', markersize=10, label='default threshold(50%)') axes[1].plot(fpr_forest[default_threshold_forest], tpr_forest[default_threshold_forest], 'x', c='k', markersize=10) axes[1].legend() axes[1].set_title('ROC Curve', fontsize=14) axes[1].set_xlabel('TPR(True Positive Rate)') axes[1].set_ylabel('FPR(False Positive Rate)') plt.show()
X印でデフォルトの閾値の際の値を表しています。PR曲線は右上に広がる(適合率も再現率も高い)ほど、ROC曲線は左上に広がる(真陽性率が高く、偽陽性率が低い)ほど良いモデルであると言えます。
デフォルトの閾値ではあまり差がなかったロジスティック回帰とランダムフォレストですが、ランダムフォレストにはより良い閾値がありその閾値の元でロジスティック回帰よりも良い結果をもたらすことがPR曲線・ROC曲線の描画により初めて分かりました。
AUC
ROC曲線は左上に広がるほど良いモデルであることから、モデルの有効な指標としてROC曲線の下の領域の面積(0から1)が使われ、この指標はAUCスコアと呼ばれます。 AUCスコアはroc_auc_score関数により算出できます。
from sklearn.metrics import roc_auc_score print('AUC score(logreg):', roc_auc_score(y_test, logreg_clf.predict_proba(X_test)[:, 1])) print('AUC score(forest):', roc_auc_score(y_test, forest_clf.predict_proba(X_test)[:, 1]))
AUCスコアにより、閾値を考慮したモデルの性能を捉えられていることがわかります。
閾値の選択
ランダムフォレストにはより良い閾値の選択の余地がありそうだと分かりました。PR曲線から適合率を下げずに0.80程度まで再現率を増やすことができると読み取れます。
いくつかの閾値について、PR曲線にプロットを追加してみましょう。
fig = plt.figure(figsize=(6, 6)) ax = fig.add_subplot(111) precision_forest, recall_forest, thresholds_forest = precision_recall_curve(y_test, forest_clf.predict_proba(X_test)[:, 1]) threshold_125_forest = np.argmin(np.abs(thresholds_forest - 0.125)) threshold_22_forest = np.argmin(np.abs(thresholds_forest - 0.22)) threshold_36_forest = np.argmin(np.abs(thresholds_forest - 0.36)) default_threshold_forest = np.argmin(np.abs(thresholds_forest - 0.5)) ax.plot(recall_forest, precision_forest, label='RandomForestClassifier') ax.plot(recall_forest[default_threshold_forest], precision_forest[default_threshold_forest], 'x', c='k', markersize=10, label='threshold 50%(default)') ax.plot(recall_forest[threshold_36_forest], precision_forest[threshold_36_forest], 'o', c='k', markersize=10, fillstyle="none", label='threshold 36%') ax.plot(recall_forest[threshold_22_forest], precision_forest[threshold_22_forest], '^', c='k', markersize=10, fillstyle="none", label='threshold 22%') ax.plot(recall_forest[threshold_10_forest], precision_forest[threshold_125_forest], 's', c='k', markersize=10, fillstyle="none", label='threshold 12.5%(min)') ax.legend() ax.set_ylim([0, 1.1]) ax.set_title('PR Curve', fontsize=14) ax.set_xlabel('recall') ax.set_ylabel('precision') plt.show()
PR曲線の閾値の最小は必ずしもゼロでないことに注意が必要です。
thresholds_forest.min()
適合率と再現率のトレードオフの選択は、冒頭に述べた通りビジネス上の影響を考慮して決定する必要があります。
上記のプロットからは、「8」を必ず正しく予測することが重要なタスクでは閾値0.36を、「8」の1つ残らず予測することが必要な(予測漏れがない)タスクでは閾値0.125を、バランスをとって閾値0.22といった閾値が考えられます。
例: 閾値を0.36に設定した場合の混同行列と分類レポートの表示
これはテストデータで閾値を決定し、テストデータで評価を行っている悪い例ですが、閾値の変更により指標が変わる様子を示すために記載しておきます。
pred_forest_36 = forest_clf.predict_proba(X_test)[:, 1] >= 0.36 matrix = confusion_matrix(y_test, pred_forest_36) sns.heatmap(matrix, annot=True, cmap='Blues', cbar=False, xticklabels=targets, yticklabels=targets) print(classification_report(y_test, pred_forest_36, target_names=targets))
適合率は1.00のまま、再現率が0.68から0.82に向上していることがわかります。適切な閾値の設定により精度が向上する一例と言えそうです。
テストデータの使用に関する注意点
モデルの評価の際には、基本的に「テストデータを使って性能評価を行うのは最後の1回だけ」を意識しておかないといけません。
テストデータの予測結果を基に描画したPR曲線・ROC曲線からモデルの閾値を決定すれば、それはテストデータに合わせた閾値となってしまうため楽観的な精度が得られます。
実際にはテストデータとは別に検証データを用意し、検証データに対するPR曲線・ROC曲線を基に閾値の選択を行いましょう。
まとめ
書きたいこと書いてると記事が長くなってしまいましたね。
PR曲線やROC曲線は2値分類タスクにおいて、ビジネス上の誤判定の影響を考慮する際には、不均衡データを扱う場合に限らず使用されます。
閾値の変更を行う際は、ハイパーパラメータの決定を行う時と同様に検証データを使いましょう。初心者のうちは、テストデータの情報を使ってまやかしの精度を叩き出してしまわないよう常に気を配らないとやらかします。(自分もよくやらかします)
では、いってきます🐼
【ブログ紹介】機械学習・データサイエンティスト初心者向けのブログ
こんばんは🐼
今回はブログの紹介記事です。機械学習・データサイエンティスト系のブログって最近増えているんですが、その中でもmyフォロワー(@takahashinbun)のお気に入りのブログを1つ紹介したいと思います。
ブログ紹介
machine-learning-stock-prediction.com
このブログの筆者はデータサイエンティスト目指して勉強中ということですが、記事を拝見する限りOJT2ヶ月の私より知識量豊富です()
エンジニアライターということで、難しいことを分かりやすく伝えてくれるのが良いですね。レイアウトが見やすいのもいいです。
株価予測等、時系列分析のアプリケーションには私自身も興味がありますし、そういったアプリケーションを作る時には何かお話したいなーと思ったり。
これから大きなブログになりそうな予感がするので、先取りで紹介しました。
まとめ
Twitterもブログも最近始めましたが、初学者として一緒に頑張ってる仲間には刺激をもらえますね!有名ブログっていっぱいあるんですけど、自分が調べもの以外に能動的にブログを見に行くのって、同じ目線で共感できるブログだったりします。
自分もゆくゆくは一流のデータサイエンティストになって、このブログを本格化したいものです。それまでは、1週間1回投稿は継続していきます。
では、おやすみなさい🐼
【CBT方式】統計検定2級合格!おすすめ統計書籍の紹介
こんばんは🐼
3月1日に統計検定2級に合格しましたので、今回はおすすめ書籍の紹介などをメインにしていきたいと思います。
対象読者はデータサイエンティストを目指す人向けに書いてますが、統計を勉強したいなら書籍はオススメなので参考にしてみてください!
試験概要
統計検定2級・3級はCBT方式の受験があり、試験会場が空いている日であれば年中いつでも受けることができます。
持ち物の電卓だけは忘れないように!
なお電卓以外は持ち込み不可です。筆記用具は用意されています。筆記試験とは少し勝手が違うので、過去問をやってみて少し時間が余るぐらいが理想です。
統計検定2級を受けるメリット
- 統計の基礎力が付いたことの確認ができる。
- 資格専用の勉強が不要。(統計の理解に集中できる)
- メジャーな資格なので、客観的評価を得られる。
- 年中受験可能(CBT)
データサイエンティストを目指す上で、統計の基本的な考え方を身に着けることはマストだと個人的には思います。
統計検定2級は、データサイエンティストに必要な統計スキルを身に着けていく過程で自然に取れるようになっているので、資格専用の勉強はほとんど必要ありません。(ただし時間配分・最終確認に過去問はやりましょう)
資格の合格よりも、統計を楽しむことを考えて勉強に取り組んだ方が後々実りは多いかと思います。
おすすめの統計書籍
コアテキスト統計学
コア・テキスト統計学 (ライブラリ経済学コア・テキスト&最先端)
- 作者: 大屋幸輔
- 出版社/メーカー: 新世社
- 発売日: 2012/01/01
- メディア: 単行本
- 購入: 1人
- この商品を含むブログを見る
東大出版の統計学入門(通称赤本)よりも優しく、網羅性があり一番統計検定2級のレベルとマッチした入門書になっていると思います。
平均・分散・標準偏差等の基本統計量から始まり、分布を学び、ゴールは検定・推定の考え方を理解することです。細かな証明は後で構いません。
正直これ一冊ちゃんと理解すれば、十分合格に届きます。
統計学入門
- 作者: 小島寛之
- 出版社/メーカー: ダイヤモンド社
- 発売日: 2006/09/28
- メディア: 単行本(ソフトカバー)
- 購入: 215人 クリック: 3,105回
- この商品を含むブログ (115件) を見る
上の書籍が難しいと感じた人は、こちらの書籍をオススメします。冗長な説明がまどろっこしく感じる所もあるかもしれませんが、最小限の労力で検定の基本的な考え方を身につけることができます。
この書籍だけでは知識としては不十分ですが、入門者の足がかりには最適かもしれません。
統計学入門(東大出版)
- 作者: 東京大学教養学部統計学教室
- 出版社/メーカー: 東京大学出版会
- 発売日: 1991/07/09
- メディア: 単行本
- 購入: 158人 クリック: 3,604回
- この商品を含むブログ (79件) を見る
統計検定2級を合格する上ではかなりのオーバースペックですが、有名なだけあり、統計の基本的な内容が厳密さを保って綺麗にまとめられています。 最初の一冊にはおすすめしませんが、基礎が身についた上でのステップアップには良いと思います。
番外編(計量経済学)
最後にスタンダートな統計に加えて、学ぶと面白いものとして計量経済学・因果推論を紹介します。
「広告を出すと、売上があがる」「気温があがると、アイスクリームが売れる」「朝ごはんを食べると、成績があがる」などなど、世の中には「Xが上がると、Yも上がる(Yが下がる)」のような、 本当かどうか分からない因果関係が多数存在します。
このような因果関係の正しい見方、実験・評価の方法を学ぶことで、データから正しい因果を導き出すことができます。
データから得られた因果を元に政策・施策等に役立てよう!ということが計量経済学の分野では活発に行われています。
「原因と結果」の経済学
- 作者: 中室牧子,津川友介
- 出版社/メーカー: ダイヤモンド社
- 発売日: 2017/02/17
- メディア: 単行本(ソフトカバー)
- この商品を含むブログ (15件) を見る
データ分析の力 因果関係に迫る思考法
- 作者: 伊藤公一朗
- 出版社/メーカー: 光文社
- 発売日: 2017/04/18
- メディア: 新書
- この商品を含むブログ (13件) を見る
これら2冊は、
- 正しく因果関係を捉えるための視点
- 計量経済学が現実にどのように使われているのか
について実際の事例をもとに説明しています。読み物として読めるので、電車の中で読むのにおすすめです。
計量経済学の第一歩
計量経済学の第一歩 -- 実証分析のススメ (有斐閣ストゥディア)
- 作者: 田中隆一
- 出版社/メーカー: 有斐閣
- 発売日: 2015/12/17
- メディア: 単行本(ソフトカバー)
- この商品を含むブログ (2件) を見る
一歩進んで、数式を交えながらより深く計量経済学を学びたい時の入門書にはこちらの書籍がおすすめです。単回帰分析や重回帰分析といった馴染みのある手法も、因果関係に着目するとより理解が深まります。
書籍のレベル感としては統計検定2級よりも少し高い程度ですが、スタンダートな統計の部分は急ぎ足なので、別書籍で補うのが良いと思います。
Rによる実証分析
- 作者: 星野匡郎,田中久稔
- 出版社/メーカー: オーム社
- 発売日: 2016/10/26
- メディア: 単行本(ソフトカバー)
- この商品を含むブログを見る
計量経済学の第一歩よりも数式を控えめに、R言語(プログラミング言語)でデータを読み込んでシミュレーションすることで、実践を通して理解することができます。
難易度としてはこちらの方がやや優しいので、PC環境があるなら因果推論・実証分析の入門としておすすめできます。
Pythonが機械学習系で主流になっていますが、統計解析・モデリングでR言語を使えたほうが役立つような場面もあるので、データサイエンティストとしてはR言語に慣れておいて損はないと思います。
まとめ
本日はおすすめの書籍の紹介をメインに、統計検定2級合格記を書かせていただきました。
まとめると、「統計を統計検定2級のためだけに勉強するのはもったいないぞ」っていうことです。
来る6月には統計検定準1級があるので、8月頃にそちらの合格記が上がらなかった時は察してください\(^o^)/6月に統計検定受ける方は一緒に頑張りましょう。
では、おやすみなさい🐼
PR曲線とROC曲線を理解する(前編)
こんばんは🐼
今回の記事は、2クラス分類タスクの際にしばしば登場するPR曲線とROC曲線について理解することを目的とします。機械学習を学ぶ上で最初に重要になるのは、結果をしっかりと評価出来るようになることだと思います。
前編ではPR曲線・ROC曲線を理解するために、必要な2クラス分類の性能評価指標のお話をしたいと思います。
認識精度の罠
機械学習の分類タスクでは、認識精度という指標が最も広く使われています。認識精度は「全体の予測のうち何%が正解したか」という最も単純で理解のしやすい指標ですが、認識精度だけで機械学習のモデル性能を評価すると不十分なことがあります。
ガンの診断
このトピックで最もよく使用されているであろう例として、ガンの陽性・陰性を判断するタスクがあります。一般にガンを患っている人はわずかなので、このような分類タスクの正解は、サンプルの選別等を行わない場合、陽性が1%・陰性が99%のように極端に偏ることになります。
ここで2種類のクラス分類器を使って、1000人(陽性10人陰性990人)に対して分類を行うとしましょう。どちらの分類器の方が優れているでしょうか。
分類器1
認識精度 99%のクラス分類器
分類器2
認識精度 95%のクラス分類器
... 何の疑いもなく「分類器1の方が優れている!」と考えてしまうのは危険です。結果を詳細に見れば、認識精度だけで判断するのが不十分なことが一目で分かるでしょう。
分類器1
正しく予想できた人数
- 陽性: 0人/ 10人(0%)
- 陰性: 990人/990人(100%)
- 全体: 990人/1000人(99%)
分類器2
正しく予想できた人数
- 陽性: 8人/ 10人(80%)
- 陰性: 942人/990人(95.2%)
- 全体: 950人/1000人(95%)
このようなタスクでは、全員を「陰性である」と予想すれば高い認識精度を誇ることができます。しかし陽性の患者を1人も発見できないので、全く使い物になりません。
当然といえば当然な気がしますが、初心者は見落としがちなポイントです。
ニ値分類の性能指標
認識精度という指標だけではクラス分類器の性能を評価するのに不十分なことが分かりました。このような時には、複数の指標でモデルを評価することが重要です。
そこでまず、TP・TN・FP・FNの意味を抑えておきましょう。英語の意味を考えれば、何ら混乱することはありません。ちなみに"陽クラス"とは、分類の際に注目しているクラスであって、ガンの陽性に限らず一般のニ値分類タスクで使われる用語です。
- TP(true positive: 真陽性): 正しく予測された陽クラスの観測値の数
- TN(true negative: 真陰性): 正しく予測された陰クラスの観測値の数
- FP(false positive: 偽陽性): 陽クラスだと予測されたが、実際には陰クラスの観測値の数
- FN(false negative: 偽陰性): 陰クラスだと予測されたが、実際には陽クラスの観測値の数
これらの数を使って、さまざまな性能指標を定義することができます。例えば、最もよく用いられる認識精度の指標は、
で表すことができます。
適合率(precision)
適合率は、陽性であると予測されたすべての観測値のうち、実際に陽性であった観測値の割合です。
今回のテーマの1つであるPR曲線のPはPrecision(適合率)を表しています。
適合率が高いほど、陽性であることが確実な場合にしか陽性であると予測しない傾向にあるといえます。適合率が重要なタスクの例として、スパムメールの分類があります。
「確実にスパムメール(陽クラス)である」といえるメール以外をスパムメールと予測すると、それだけスパムでないメール(陰クラス)をスパムメールとして弾いてしまう可能性が高まるからです。
再現率(recall)・真陽性率(TPR : true positive rate)
再現率は、実際に陽性である観測値のうち、陽性であると予想された観測値の割合です。
今回のテーマの1つであるPR曲線のPはRecall(再現率)を表しています。また、ROC曲線では真陽性率(true positive rate)と呼ばれますが、内容は再現率と全く同じです。
再現率が高いほど、陽性であるかどうか分からない場合にも陽性であると予測する傾向にあるといえます。再現率が重要なタスクの例として、ガンのスクリーニングがあります。
ガンの発見のようなタスクでは、健康な患者を陽性であると診断してしまうリスクよりも、ガンを患っている患者を陰性であると診断してしまうリスクの方が圧倒的に高いからです。(陽性であると分類された患者に対しては、医師が直接確認することができる)
F1値(F1 value)
適合率と再現率には後述するようにトレードオフの関係があり、一方の値だけでなくこれらの指標のバランスが重要です。
そこで、適合率と再現率の調和平均をとったF1値という指標が使われることもあります。
適合率と再現率のトレードオフ
適合率と再現率にはトレードオフの関係があります。前述の通り、タスクによって適合率と再現率のどちらを優先すべきかは変わるため、タスク毎に最適なバランスを考慮する必要があります。その際にPR曲線は役立ちます。
前編まとめ
今回は随時例を挙げながら、2クラス分類の性能指標のお話をメインにしました。
現実には完璧な分類器を作ることは不可能なので、2クラス分類では分類器の性能を上げる以外にも、ビジネス的な視点から適合率と再現率のトレードオフについて考える必要があります。
後編では実際にサンプルコードを実行しながら、今回説明した指標の評価を視覚的に行うことができる混同行列・PR曲線・ROC曲線といった道具について理解を深めたいと思います。
では、おやすみなさい🐼
【Pandas】初見データ確認の心得と使えるメソッド!
こんばんは🐼
今日の記事は初見データの確認について。データを扱う以上、どんな分析を行うにも最初はまずデータを眺めることになると思います。
Pandasを使ってcsvファイルやpickleファイルやなんやらからデータフレームに読み込んだはいいものの、初手どうするの!?って人に向けて自分も初心者ながら、参考にしていただきたくこの記事を書きます。
下記のサンプルコードを動作させるには、Pandasライブラリが必要です。
使用するデータ
この記事では、下記のオンライン購買データ(公開データセット)を使用します。下記のコードで23.7MBのファイルを、ネット経由でローカルに保存して使用します。 (実行環境は、Jupyter notebookを想定しています)
import pandas as pd import requests url = "http://archive.ics.uci.edu/ml/machine-learning-databases/00352/Online%20Retail.xlsx" response = requests.get(url) # カレントディレクトリに"onlineRetail.xls"という名前で保存 output = open("onlineRetail.xls", 'wb') output.write(response.content) output.close() # データフレームの作成 transaction = pd.ExcelFile("onlineRetail.xls") transaction_table = transaction.parse('Online Retail')
データの表示
まずはデータフレームにデータが読み込めているか確認しましょう。head()メソッドを使うとデフォルトでデータフレームの先頭5件を表示できます。
# 先頭5件を表示
transaction_table.head()
引数に数を指定すると、その先頭からn行取り出すことができます。
# 先頭10件を表示 transaction_table.head(10)
末尾を取り出したい場合、tail()メソッドを使用しましょう。
# 末尾5件を表示
transaction_table.tail()
データを読み込んだり、変更を加えたりした時は逐一確認することが基本です。
そのための一番手軽な方法が、head()メソッドやtail()メソッドを使うことです。(先頭や末尾を確認するだけでは不十分な場合は、適切に参照を行う必要があります)
メタ情報の表示
info()メソッドを使用することで、データフレームのメタ情報を取得できます。
transaction_table.info()
この出力には、以下の情報が含まれていることがわかります。
- インデックスの総数
- カラム数
- 各カラムの欠損値を除いたデータ数
- 各カラムのデータ型
- データフレームのメモリの使用量
データフレームを読み込んだら、まずは全体の情報から確認して細部を確認するのが基本です。
「そもそもデータは足りているのか?」「データ型はそのままでよいのか?」「欠損はどの程度あるのか?」・・・などは一番最初に確認してしまいましょう。
各カラムの欠損値の数
isnull()メソッドを使用することで、欠損値が含まれているかどうかの真理値を得ることができます。sum()メソッドと組み合わせることで、各カラムの欠損値を取得できます。
transaction_table.isnull().sum()
info()メソッドの結果から計算することもできますが、欠損値の扱いは前処理の中でも特に重要なので、各カラムに欠損がいくつ含まれているか明確に把握しておくのがオススメです。
要約統計量の表示
describe()メソッドを使用することで、各カラムの要約統計量を表示できます。
transaction_table.describe()
デフォルトでは数値型を持つカラムの要約統計量が表示されますが、include='all'オプションを付けることで、全てのカラムの要約統計量を表示することができます。
transaction_table.describe(include='all')
数値型では平均・標準偏差・最小・最大・四分位数といった情報が表示されますが、オブジェクト型(文字列型など)や日付型では、ユニーク数・最頻値とその頻度等の情報が表示されます。
要約統計量で各カラムの「値」のおおまかな特徴が掴めます。
ここでは、CustomerIDが数値として扱われていることに違和感を感じるでしょう。型の変換はastype()メソッドを使うことで行なえます。
astype()メソッドは元のシリーズを書き換えるわけではないので、データフレームに変更を加えたい場合は以下のように代入を行いましょう。
transaction_table['CustomerID'] = transaction_table['CustomerID'].astype('str')
ただしこれでは一つ問題があります。欠損値を表す「NaN」が文字列の"nan"として扱われてしまいます。下記のように欠損値以外をキャストするか、最初から扱いたい型が決まっている場合はデータフレームの読み込み時のオプションで型を指定しましょう。
transaction_table['CustomerID'] = transaction_table[transaction_table['CustomerID'].notnull()]['CustomerID'].astype('str')
データ型の詳細確認
info()メソッドを使うことで、各カラムのデータ型を確認することができましたが、"object"型には複数のデータ型が含まれている可能性があることに注意しましょう。
各カラムのデータ型の詳細情報を確認するのは、以下のような手法で行うことができます。
transaction_table['InvoiceNo'].map(type(str)).value_counts()
InvoiceNoのカラムには、str型とint型が混在していることが分かりました。
先程と同様にキャストを行うことでデータ型を揃えることができますが、このデータフレームではStockCodeやDescriptionも同様なので、データフレーム読み込み時に指定するのが一番楽でしょう。
各カラムのユニークレコードの取得
describe()メソッドの時にも、ユニーク数という言葉が出てきましたがユニークとは「一意な」という意味で、要はカラムに何種類のレコードが存在するかということです。
全体を確認した後は、カラムごとに具体的にどんな値が入っているのかを確認することも重要です。その用途ではunique()メソッドが使用できます。
unique()メソッドを使用することで、重複を除いた一意なレコードを表示できます。
pd.Series(transaction_table['Description'].unique())
一意なレコードと共に出現数まで表示したい場合は、value_counts()を使用します。出現数順にpandas.Series型で結果を返してくれるため、自分はこちらをよく利用します。
まとめ
Pandasを使ったデータの初手確認方法の最低限を紹介しました!勿論これが正解ではありませんし、この後に様々な角度から基礎集計を行う必要があります。
データ確認や基礎集計の技術は、あらゆるデータを扱う場合に必要になってくるので、最初に固めておくと役立つと思います。
特に重要なのは膨大なデータであっても、
- 欠損値や異常値を見つけられること
- データにどのような値が入っているのか、どのような型で扱うのか、など全体を掴むこと
なのだと思います。データの外観を掴まないまま基礎集計等を行っても結局は2度手間を踏むか、最悪誤った分析結果を導きだしてしまうかもしれません。
またデータ確認や基礎集計の過程で、データを壊してしまうのも注意しなくてはなりません。特にデータの変更・削除等の処理を行ったなら、データが意図した通りに変更されているか確認しましょう。
分析の経験のある人には当たり前のことだとは思いますが、初心者視点から書かせていただきました。データの確認の習得は、データ分析をスマートに行う第一歩だと思います!
では、おやすみなさい🐼
正則化の意味を理解する
「過学習」に陥らないようにモデルを構築することは、機械学習(教師有り学習)の基本です。正則化は過学習を防ぐテクニックの一つです。Weight Decayとも呼ばれ、DeepLearningなど最新の手法でも使われています。
正則化は難しい概念ではないので、実際に使ってみるとわかります。 下記のサンプルコードを動作させるには、
- Scikit-learn
- NumPy
- Matplotlib
の3つのライブラリが必要です。
今回は線形回帰を例に、Scikit-learnのライブラリを使って、正則化の挙動を確認します。その前にデータを用意しましょう。今回は描画を行うために、2次元のデータセットで確認します。
但し、一般には特徴量が多い(描画できない)場合に、過学習が起きやすい傾向があることに留意しましょう。
import numpy as np import matplotlib.pyplot as plt plt.style.use('ggplot') rng = np.random.RandomState(1) x = 10 * rng.rand(50) y = np.sin(x) + 0.1 * rng.randn(50) fig = plt.figure() ax = fig.add_subplot(111) ax.scatter(x, y) ax.set_xlabel('x') ax.set_ylabel('y') plt.show()
このトイデータセットを使ってxを説明変数、yを目的変数とする回帰を行います。ここで、過学習が起こる条件の1つは、特徴量が多くモデルの表現力が高いことです。
今回は、xについて基底関数を適用することでモデルを複雑さ(表現力)を意図的に増やして、過学習を起こします。例えば、xに関する高次の項を用意するために、scikit-learnのPolynomialFeatures変換器が利用できます。
のように、xの3次の項まで用意したい場合は、PolynomialFeaturesのdegreeオプションに3を指定してインスタンスを生成すれば良いです。
poly = PolynomialFeatures(degree=3) # 切片項が不要な場合、include_bias=Falseを指定する。 # poly = PolynomialFeatures(degree=3, include_bias=False) x_deg_3 = poly.fit_transform(x[:, np.newaxis])
今回はガウス基底関数の和をモデルにすることで、より顕著な過学習を確認することが出来るため、ガウス基底関数を利用します。ガウス基底関数はscikit-learnに組み込まれていないので、自前で用意する必要があります。 (ガウス分布=正規分布のことです。xの区間をN等分して、N個の中心を持つ正規分布の和をモデルとします。)
下記のコードは、書籍「Pythonデータサイエンスハンドブック」の引用です。
from sklearn.base import BaseEstimator, TransformerMixin class GaussianFeatures(BaseEstimator, TransformerMixin): #1次元入力の等間隔ガウス関数 def __init__(self, N, width_factor=2.0): self.N = N self.width_factor = width_factor @staticmethod def _gauss_basis(x, y, width, axis=None): arg = (x - y) / width return np.exp(-0.5 * np.sum(arg ** 2, axis)) def fit(self, X, y=None): # データ範囲に沿ってN個の中心を作成する self.centers_ = np.linspace(X.min(), X.max(), self.N) self.width_ = self.width_factor * (self.centers_[1] - self.centers_[0]) return self def transform(self, X): return self._gauss_basis(X[:, :, np.newaxis], self.centers_, self.width_, axis=1)
正則化を行わずに、30次のガウス基底関数の和を使って、回帰を行ってみましょう。
from sklearn.pipeline import make_pipeline from sklearn.linear_model import LinearRegression model = make_pipeline(GaussianFeatures(30), LinearRegression()) model.fit(x[:, np.newaxis], y) xfit = np.linspace(0, 10, 1000) # 回帰直線プロット用 y_pred = model.predict(xfit[:, np.newaxis]) fig = plt.figure() ax = fig.add_subplot(111) ax.scatter(x, y) ax.plot(xfit, y_pred) ax.set_xlabel('x') ax.set_ylabel('y') ax2 = fig.add_subplot(212) ax2.plot(np.linspace(0, 10, 30), model.steps[1][1].coef_) ax2.set_xlabel('basis') ax2.set_ylabel('coefficient') plt.show()
上が回帰プロット、下が30次のガウス基底それぞれの係数プロットです。確かに全ての散布図データにフィットするような回帰曲線が引かれていますが、未知のデータに対応できるモデルとは言えないでしょう。機械学習(特に教師有り学習)では、汎化性能のみに興味があり、既に得られているデータへの当てはまりは全く重要でない場合が多いというところがポイントです。これが、過学習が起きている状態です。
係数プロットからはかなり大きい係数の値を取っていることが読み取れます。
このように複雑なモデルでは、大きな係数を取ることで学習データに対してのみ上手くフィットする極端な形のグラフがかけてしまうのです。
そこで、係数にペナルティをかけて、極端な形のグラフをかけないようにしようというのが正則化のアイデアです。
数式的には、損失関数に正規化項を加えます。
L1正則化であれば、
が使われます。
回帰係数は、損失関数を最小にするように求めるため、残差を小さくする他に、正則化項も小さくすること、つまり回帰係数を小さくすることが要求されるのです。
主要な線形回帰モデル
scikit-learnの主要な線形回帰器は4つあります。正則化を行わないLinearRegression、L2正則化を行うRidge、L1正則化を行うLasso、L2正則化とL1正則化のハイブリッドElasticNetです。それぞれ表にまとめます。
モデル | 正則化 | 正則化パラメータ(デフォルト) |
---|---|---|
LinearRegression | 正則化なし | なし (正則化ありモデルのalpha=0に相当) |
Ridge | L2正則化 | alpha=1 |
Lasso | L1正則化 | alpha=1 |
ElasticNet | L1 + L2正則化 | alpha=1, l1_ratio=0.5 |
Ridge
Ridge回帰(L2正則化)は、最も一般的な正則化の形式です。回帰係数の二乗和(L2ノルム)を正則化項に加えることで、係数を抑えます。
この記事では行いませんが、全ての正則化に共通することとして、正則化の強さを決めるハイパーパラメータを最適化することが重要です。
from sklearn.linear_model import Ridge model = make_pipeline(GaussianFeatures(30), Ridge(alpha=0.1)) model.fit(x[:, np.newaxis], y) xfit = np.linspace(0, 10, 1000) # 回帰直線プロット用 y_pred = model.predict(xfit[:, np.newaxis]) fig = plt.figure(figsize=(8, 8)) ax = fig.add_subplot(211) ax.scatter(x, y) ax.plot(xfit, y_pred) ax.scatter(6, model.predict(np.array([[6]])), marker='*', s=200) ax.set_xlabel('x') ax.set_ylabel('y') ax2 = fig.add_subplot(212) ax2.plot(np.linspace(0, 10, 30), model.steps[1][1].coef_) ax2.set_xlabel('basis') ax2.set_ylabel('coefficient') plt.show()
正則化項を追加することで、係数のスケールが大幅に小さくなり、感覚的に正しいグラフが得られました。
Lasso
Lasso回帰(L1正則化)は、回帰係数の絶対値の和(L1ノルム)を正則化項に加えることで、係数を抑えます。
L1正則化の重要な特性として、不要な特徴量の回帰係数を(正則化の強さ次第で)0にできることが挙げられます。そのため、特徴選択の用途で使われることがあります。
from sklearn.linear_model import Lasso model = make_pipeline(GaussianFeatures(30), Lasso(alpha=0.001, max_iter=100000)) model.fit(x[:, np.newaxis], y) xfit = np.linspace(0, 10, 1000) # 回帰直線プロット用 y_pred = model.predict(xfit[:, np.newaxis]) fig = plt.figure(figsize=(8, 8)) ax = fig.add_subplot(211) ax.scatter(x, y) ax.plot(xfit, y_pred) ax.scatter(6, model.predict(np.array([[6]])), marker='*', s=200) ax.set_xlabel('x') ax.set_ylabel('y') ax2 = fig.add_subplot(212) ax2.plot(np.linspace(0, 10, 30), model.steps[1][1].coef_) ax2.set_xlabel('basis') ax2.set_ylabel('coefficient') plt.show()
いくつかの係数が、完全にゼロになっていることが確認できます。
ElasticNet
Elastic Netは、L1正則化とL2正則化の両方を使用することで、より良い正則化を目指したモデルです。
Elastic Netには正則化の強さを決めるパラメータが2つあり、1つは全体としての正則化の強さを決めるalpha、もう1つはL1正則化とL2正則化の強さの比率を決めるl1_ratioがあります。
ハイパーパラメータを上手く調整すれば、Ridge回帰やLasso回帰よりも高い性能が得られることがある一方、適切な正則化パラメータを見つける労力がかかるというデメリットもあります。
from sklearn.linear_model import ElasticNet model = make_pipeline(GaussianFeatures(30), ElasticNet(alpha=0.001, l1_ratio=0.5, max_iter=100000)) model.fit(x[:, np.newaxis], y) xfit = np.linspace(0, 10, 1000) # 回帰直線プロット用 y_pred = model.predict(xfit[:, np.newaxis]) fig = plt.figure(figsize=(8, 8)) ax = fig.add_subplot(211) ax.scatter(x, y) ax.plot(xfit, y_pred) ax.scatter(6, model.predict(np.array([[6]])), marker='*', s=200) ax.set_xlabel('x') ax.set_ylabel('y') ax2 = fig.add_subplot(212) ax2.plot(np.linspace(0, 10, 30), model.steps[1][1].coef_) ax2.set_xlabel('basis') ax2.set_ylabel('coefficient') plt.show()
Lasso回帰と比較して、いくつかの特徴量の係数が0ではなくなっていることがわかります。
まとめ
今回は復習も兼ねて、正則化についてまとめました。実際のデータ分析で過学習が起きる状況は、特徴量が多数あり、少なくともそのままでは回帰の様子を可視化できないことがほとんどだと思います。そのため、2次で過学習のイメージを持っておくことは、正則化の直感的な理解に役立つと思います。
この記事では回帰の例を使って説明しましたが、ロジスティック回帰などの線形分類の正則化項についても同様の考え方で理解することができます。ただし、ロジスティック回帰の正則化パラメータのCは、の逆数であるので、Cを大きくするほど正則化が弱まることには注意が必要です。
(線形サポートベクターマシンについては、Cはソフトマージンのパラメータなので、正則化とは異なりますが、モデルのバイアス・バリアンスを制御するという意味では同じです。)
ディープラーニングについても、損失関数を最小化するように係数(重み)を更新するという点については同じであり、全ての重みからなるノルムを損失関数に加えることで、重みにペナルティをかけて過学習を防ぐことを期待しています。
では、おやすみなさい🐼
データサイエンティストへの道
どこぞの無名新米データサイエンティストに興味を持つ人はいないと思うので、初投稿では明瞭簡潔にこのブログの趣旨を書こうと思います。
当ブログで扱う内容は、
まとめると、データサイエンティストに必要な技術的スキルセットです。
ブログの目的は、
- 自分が学習したことをわかりやすく丁寧にアウトプットすること。
- オリジナルの課題を見つけ、分析結果や機械学習システムを公開するまでの過程を残すこと。
ってとこです。
コンテンツで有名になれたら、その時は自己紹介します。
では、おやすみなさい🐼