Chapter 3 Classification
Hands-On Machine Learning with Scikit-Learn, Keras & Tensorflow 2nd Edition by A. Geron
MNIST
MNISTの説明：手書き数字のデータベースで、機械学習の"hello world"である。
Scikit-Learnからダウンロードできる。
from sklearn.datasets import fetch_openml
mnist = fetch_openml("mnist_784", version=1)
mnist.keys( )
dict_keys(["data", "target", "feature_names", "DESCR", "details", "categories", "url"])
X, y = mnist["data"], mnist["target"]
X.shape
(70000, 784)
y.shape
(70000,)
28 x 28 pixels, from 0 (white) to 255 (black)
画像の表示
import matplotlib as mpl
import matplotlib.pyplot as plt
some_digit = X[0]：後で何度も使われる
some_digit_image = some_digit.reshape(28, 28)
plt.imshow(some_digit_mage, cmap="binary")
plt.axis("off")
plt.show( )
y[0]
'5'
MNISTは、訓練用とテスト用にすでに分けられており、シャッフルされており、cross-validationにも対応（分割してもその中でのラベルの分布は同じ）している。
X_train, X_test, y_train, y_test = X[ : 60000], X[60000 : ], y[ : 60000], y[60000 : ]
Training a Binary Classifier
問題を簡単にするために、1つの数字、例えば５、を認識することにしよう。５の検出、すなわち５と５以外の数字に分類する。
ここで、分類のためのラベル target vector を作る。
y_train_5 = (y_train == 5) # True for all 5s, False for all other digits
y_test_5 = (y_test == 5)
y_train # ラベルを表示する
array([5, 0, 4, ..., 5, 6, 8], dtype=uint8)
y_train_5 # ５と５以外の数字を分類するためのラベルを表示する。
array([ True, False, False, ..., True, False, False])
まずは、Stochastic Gradient Descent (SGD) classifier を使う。このclassifierは、非常に大きなデータセットに対して効率よく動作する。online learningにも適している。
60000の全トレーニングセットに適用する。
from sklearn.linear_model import SGDClassifier
sgd_clf = SGDClassifier(random_state=42)
sgd_clf.fit(X_train, y_train_5)
sgd_clf.predict([some_digit]) # some_digit = X[0], ---> sgd_clf_predict([X[0]])
array([True])
Performance Measures
分類モデルの性能評価は、回帰モデルの性能評価よりも、トリッキーなので、多くのスペースを割いて説明する、とのこと。
Measuring Accuracy Using Cross-Validation
ここでも、cross validation が有効である。
from sklearn.linear_model import SGDClassifier
from sklearn.model_selection import cross_val_score
sgd_clf = SGDClassifier(randomstate=42)
cross_val_score(sgd_clf, X_train, y_train_5, cv=3, scoring="accuracy")
array([0.96355, 0.93795, 0.95615)
凄いな、と思う前に、どれも５ではない、という答えしかできないモデルを使った分類モデルを、cross validationによって評価してみよう。
from sklearn.base import BaseEstimator
class Never5Classifier(BaseEstimator):
def fit(self, X, y=none):
return self
def predict(self, X):
return np.zeros*1
ovr_clf.fit(X_train, y_train)
ovr_clf.predict([some_digit])
array([5], dtype=unit8)
len(ovr_clf.estimators_)
10
SGDClassifierをtraining
sgd_clf.fit(X_train, y_train)
sgd_clf.predict([some_digit])
array([5], dtype=unit8)
SGD classifierは、もともとmultiple classifierなので、OvRやOvOは動作しない。
decision_function( )により、スコアを出力してみよう。
sgd_clf.decision_function([some_digit])
array(-15955, -38080, -13326, 573, -17680, 2412, -25526, -12290, -7946, -10631)
０から９までのスコアであり、５のスコアが最大値を示していることがわかる。
それではcross_val_score( )を計算してみよう。
cross_val_score(sgd_clf, X_train, y_train, cv=3, scoring="accuracy")
array([0.849, 0.871, 0.870])
この１から９までの分類においては、ランダム分類では10％のスコアだから、85％はまずまずの値である。
2章で学んだscaling the inputを行ってみよう。
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler( )
X_train_scaled = scaler.fit_transform(X_train.astype(np.float64))
cross_val_score(sgd_clf, X_train_scaled, y_train, cv=3, scoring="accuracy")
array([0.897, 0.896, 0.907])
入力データをスケーリングするだけで、84％から90％までスコアが上がった。
こういうデータ処理scalingは、忘れないようにしよう！
Error Analysis
もしこれが実際のプロジェクトであれば、チェックリストに従って進めているところだろう。複数のモデルを試して、候補が絞られてきたら、GridSearchCVを使ってhyperparameterをfine-tuningするところだろう。有望なモデルが見つかったら、それをさらに改善するには、どういうエラーをしているかを調べるのも１つの方法である。
confusion matrixを調べてみよう。
cross_val_predict( )で予測してから、confusion_matrix( ) functionを呼び出せばよい。
y_train_pred = cross_val_predict(sgd_clf, X_train_scaled, y_train, cv=3)
conf_mx = confusion_matrix(y_train, y_train_pred)
conf_mx
array([[5577, 0, 22, 5, 8, 43, 36, 6, 225, 1],
[ 0, 6400, 37, 24, 4, 44, 4, 7, 212, 10],
[ 27, 27, 5220, 92, 73, 27, 67, 36, 378, 11],
[ 22, 17, 117, 5227, 2, 203, 27, 40, 403, 73],
[ 12, 14, 41, 9, 5182, 12, 34, 27, 347, 164],
[ 27, 15, 30, 168, 53, 4444, 75, 14, 535, 60],
[ 30, 15, 42, 3, 44, 97, 5552, 3, 131, 1],
[ 21, 10, 51, 30, 49, 12, 3, 5684, 195, 210],
[ 17, 63, 48, 86, 3, 126, 25, 10, 5429, 44], [ 25, 18, 30, 64, 118, 36, 1, 179, 371, 5107]], dtype=int64)
図で表示すると、
plt.matshow(conf_mx, cmap=plt.cm.gray)
plt.show( )
５の位置のイメージが他より若干暗い。 ５のクラスの数が少ないか、または、分類精度が良くない場合に暗くなるのだが、今回は両方に当てはまっている。
誤認にfocusするために、error rateをプロットする。
row_sums = conf_mx.sum(axis=1, keepdims=True)
norm_conf_mx = conf_mx / row_sums
np.fill_diagonal(norm_conf_mx, 0)
plt.matshow(norm_conf_mx, cmap=plt.cm.gray)
plt.show( )
８の列が明るいのは、8に誤分類される割合が高いということを示唆している。
class ８は、特に明るい領域は無いので、適切に分類されているようである。
３と５の明暗が対照的であることから、互いに誤認されていることがわかる。
3に分類されたものと5に分類されたを、クラス３とクラス５について表示する。
cl_a, cl_b = 3, 5
X_aa = X_train[(y_train == cl_a) & (y_train_pred == cl_a)]
X_ab = X_train[(y_train == cl_a) & (y_train_pred == cl_b)]
X_ba = X_train[(y_train == cl_b) & (y_train_pred == cl_a)]
X_bb = X_train[(y_train == cl_b) & (y_train_pred == cl_b)]
plt.figure(figsize=(8,8))
plt.subplot(221); plot_digits(X_aa[:25], images_per_row=5)
plt.subplot(222); plot_digits(X_ab[:25], images_per_row=5)
plt.subplot(223); plot_digits(X_ba[:25], images_per_row=5)
plt.subplot(224); plot_digits(X_bb[:25], images_per_row=5)
save_fig("error_analysis_digits_plot")
plt.show()
左半分が３に分類され、右半分が5に分類されたものである。
大半が、明らかな分類ミスのように見えるので、ここから改善策を見つけるのは難しそうである。
前処理として、文字の傾きや回転を修正する、センタリングする、というのは効果があるかもしれない。
Multilabel Classification
顔認証の場合は、1枚の画像中に、複数あるうちのどの顔があるかを予測する。
例題は、MNISTデータベースの場合で、１つのinstanceの１つの文字に対して、7以上か否かと、奇数か否かを予測するものである。
この２つのケースはまったく異なっているように見えるが、１つのinstanceに対して複数のラベルを与えるということにおいては同じである。
コードをみてみよう。
from sklearn.neighbors import KNeighborsClassifier
y_train_large = (y_train >= 7)
y_train_odd = (y_train % 2 == 1)
y_multilabel = np.c_[y_train_large, y_train_odd]
knn_clf = KNeighborsClassifier( )
knn_clf.fit(X_train, y_multilable)
KNeighborsClassifierは、multiclassificationをサポートしているが、すべてのclassifierがmulticlassificationに対応しているわけではない。
予測してみよう。
knn_clf.predict([some_digit])
array(False, True)
some_difitは５なので、not large (False), odd (True)ということで正しく学習し、予測している。
スコアを算出してみよう。
y_train_knn_pred = cross_val_predict(knn_clf, X_train, y_multilable, cv=3)
f1_score(y_multilable, y_train_knn_pred, average="macro")
0.97641
lableに重みづけしてスコアを算出するには、スコアの計算で、average="weighted"とすればよい。
Multioutput Classification
これは、興味深い。
文字画像に下図に示すようにノイズをのせる。
noise = np.random.randint(0, 100, (len(X_train), 784))
X_train_mod = X_train + noise
noise = np.random.randint(0, 100, (len(X_test), 784))
X_test_mod = X_test + noise
y_train_mod = X_train
y_test_mod = X_test
学習は、KNeighborsClassifierを用い、ノイズの無い画像をラベルにして行う。
knn_clf.fit(X_train_mod, y_train_mod)
clean_digit = knn_clf.predict([X_test_mod[some_index]])
plot_digit(clean_digit
ノイズ除去ができるなんて、驚きだ。
デザイン・設計次第で何でもできそうな気がしてくるのが、機械学習の面白いところだ。ディープラーニングもいいが、機械学習もなかなか面白いものだなと感じてきた。
これまで、機械学習は過去の遺物で、ディープラーニング全盛の時代だと思っているが、まだまだ使えるのかな。それとも、最先端では、・・・。
Exercises
3. Tackle the Titanic dataset.
初心に戻ろうと思って5月8日にKaggleのTitanicに取り組んだ。
チュートリアルに従って、一通りの手順をおさらいした後、notebookを参考に勉強しようとしたが、ついていけるコードは学ぶことが殆どなく、良さそうなnotebookは、Rだったり、C++だったり、Pythonでもコードが理解できなかったりで、いつもの壁にぶち当たった。
そこで、A. Geronさんのテキストに取り組み始めた。
2章をクリヤ（Exercisesは未着手）し、3章もクリヤ（Exercisesは未着手）し、4章に進もうと思ったが、機械学習を習得できるか否かは、A. Geronさんの、このテキストを読み込むことができるか否かにかかっているような気がしてきた。
2章、3章と読み進めるにしたがって、このテキストの良さがわかってきた。
ならば、Exercisesにもちゃんと取り組まないと本当に理解したことにはならないのではないかと思い、Exercisesに取り組むことにした。
という状況で第3章のExercisesを眺めていたら、KaggleのTitanicが取り上げられているではないか。
これは、good timing である。
とりあえず、コードを最初の方から、眺めてみよう。
train_data = load_titanic_data("train.csv")
test_data = load_titanic_data("test.csv")
train_data.head( )
train_data.info( )
train_data.describe( )
train_data["Survived"].value_counts()
0 549
1 342
Name: Survived, dtype: int64
train_data["Pclass"].value_counts( )
3 491
1 216
2 184
Name: Pclass, dtype: int64
train_data["Sex"].value_counts( )
male 577
female 314
Name: Sex, dtype: int64
train_data["Embarked"].value_counts( )
S 644
C 168
Q 77
Name: Embarked, dtype:int64
ここから本格的な前処理が始まる。
ところで、3章で学習したClassifierは、データが単純だった。５かそれ以外の数字か、という分類であった。分類の根拠は、数字のピクセル情報に求めるしかないのだが、モデルが数字をどうやって分類しているのか、その根拠がわからなかった。
Titanicもclassifierだが、分類の根拠となりそうな情報attributeがたくさん示されている。そうすると、2章でやったhousingのように、目的とする情報との相関性を調べることが重要だと思うのだが、それは、ここでのexaciseの主題ではなく、自分で検討しなさいということのようだ。
とりあえず、Exerciseにもどろう。
使われるのは、Pipeline, FeatureUnion, DataFrameSelectorとのこと。
よくわかんないね。
pipelimeってなんだったかな。38ページの囲み記事。
A sequence of data processing components is called a data pipeline.
実例は70ページと71ページ
num_pipeline = Pipeline([
('imputer', SimpleImputer(strategy="median")),
('attribs_adder', CombinedAttributesAdder( )),
('std_scaler', StandardScaler( )),
])
housing_num_tr = num_pipeline.fit_transform(housing_num)
--------------------------------------------
num_attribs = list(housing_num)
cat_attribs = ["ocean_proximity"]
full_pipeline = ColumnTransformer([
("num", num_pipeline, num_attribs),
("cat", OneHotEncoder( ), cat_attribs),
])
housing_prepared = full_pipeline.fit_transform(housing)
ipelineはよく使いそうだが、FeatureUnionとかDataFrameSelectorはどうなんだろう。
Titanicでもpipelineは、大活躍。そこにFeatureUnionとDataFrameSelectorが出てくる。
num_pipeline = PipeLine([
("select_numeric", DataFrameSelector(["Age", "SibSup", "Parch", "Fare"])),
("imputer", SimpleImputer(strategy="median")),
])
cat_pipeline =PipeLine([
("select_cat", DataFrameSelector(["Pclass", "Sex", "Embarked"])),
("imputer", MostFrequentImputer( )),
("cat_encoder", OneHotEncoder(sparse=False)),
])
preprocess_pipeline = FeatureUnion(transformer_list=[
("num_pipeline", num_pipeline),
("cat_pipeline", cat_pipeline),
])
X_train = preprocess_pipeline.fit_transform(train_data)
こうやって並べてみると、pipelineが中心にいるのがよくわかるので、とにかく、いつでもどこでもpipelineということで、pipelineの習得に努めていく。
ここでの、Titanicの課題に対する、モデルの比較を見ておこう。
svm_scores = cross_val_score(svm_clf, X_train, y_train, cv=10)
svm_scores.mean( )
0.7365
forest_scores = cross_val_score(forest_clf, X_train, y_train, cv=10)
forest_scores.mean( )
0.8150
cross validationといっても、元のデータ数が少ないので、overfittingになる可能性があるので、このままテストデータに適用しても、0.80は超えないだろうと推測される。
Scikit-LearnのBaseEstimatorとTransformerMixinがわからない。
Exercise 4. Build a spam classifier (a more challenging exercise)
コードは追えなくても、3章までの知識で、どの程度のことができるのかは知っておきたい。attributeが、メールの構造、単語の種類、数、並びだったりするようで、その解析ために、mailやNLTKが使われる。データの解析、前処理、訓練データの作成、訓練、予測、評価、・・・。
First, let's fetch the data:
import os
import tarfile
import urllib
DOWNLOAD_ROOT = "http://spamassassin.apache.org/old/publiccorpus/"
HAM_URL = DOWNLOAD_ROOT + "20030228_easy_ham.tar.bz2"
SPAM_URL = DOWNLOAD_ROOT + "20030228_spam.tar.bz2"
SPAM_PATH = os.path.join("datasets", "spam")
def fetch_spam_data(spam_url=SPAM_URL, spam_path=SPAM_PATH):
if not os.path.isdir(spam_path):
os.makedirs(spam_path)
for filename, url in *2:
path = os.path.join(spam_path, filename)
if not os.path.isfile(path):
urllib.request.urlretrieve(url, path)
tar_bz2_file = tarfile.open(path)
tar_bz2_file.extractall(path=SPAM_PATH)
tar_bz2_file.close()
fetch_spam_data( )
let's load all the emails:
HAM_DIR = os.path.join(SPAM_PATH, "easy_ham")
SPAM_DIR = os.path.join(SPAM_PATH, "spam")
ham_filenames = [name for name in sorted(os.listdir(HAM_DIR)) if len(name) > 20]
spam_filenames = [name for name in sorted(os.listdir(SPAM_DIR)) if len(name) > 20]
len(ham_filenames)
2500
len(spam_filenames)
500
We ca use Python's email module to parse these emails (this handles headers, encoding, and so on):
import email
import email.policy
def load_email(is_spam, filename, spam_path=SPAM_PATH):
directory = "spam" if is_spam else "easy_ham"
with open(os.path.join(spam_path, directory, filename), "rb") as f:
return email.parser.BytesParser(policy=email.policy.default).parse(f)
ham_emails = [load_email(is_spam=False, filename=name) for name in ham_filenames]
spam_emails = [load_email(is_spam=True, filename=name) for name in spam_filenames]
Let's look at one example of ham and one example of spam, to get a feel of what the data looks like:
print(ham_emails[1].get_content( ).strip( )
print(spam_emails[6].get_content( ).strip())
Some emails are actually multipart, with images and attachiments (which can have their own attachments). Let's look at the various types of structures we have.
メールの構造に違いがありそうだということだが、どうやってメールの構造を解析するのだろう。
def get_email_structure(email):
if isinstance(email, str):
return email
payload = email.get_payload( )
if isinstance(payload, list):
return "multipart({})".format(", ".join([
get_email_structure(sub_email)
for sub_email in payload
]))
else:
return email.get_content_type( )
--------------------------------------------
from collections import Counter
def structures_counter(emails):
structures = Counter( )
for email in emails:
structure = get_email_structure(email)
structures[structure] += 1
return structures
structures_counter(ham_emails).most_common( )
structures_counter(spam_emails).most_common( )
It seems that the ham emails are more often plain text, while spam has quite a lot of HTML. Moreover, quite a few ham emails are sighned using PGP , while no spam is. In short, it seems that the email structure is useful information to have.
こうしている間にも、spam mailが送られてきた。
有用だと思ってこちらからメールアドレスを登録した場合でも、目的外の勧誘メールを送ってくることがある。これもspamに分類したいね。
次はヘッダーを調べている。
for header, value in spam_emails[0].items( ):
print(header, ":", value)
spam_emails[0]["Subject"]
'Life Insurance - Why Pay More?'
yahooなんかも、この種の見出しの広告であふれている。
見ないのが一番なのだが、テレビもバカな広告であふれている。
Okay, before we learn too much about the data, let's not forget to split it to a training set and a test set:
import numpy as np
from sklearn.model_selection import train_test_split
X = np.array(ham_emails + spam_emails)
y = np.array([0] * len(ham_emails) + [1] * len(spam_emaila))
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
HTMLをplain textに変換するには、BeautifulSoupe libraryというものがあるらしいが、ここでは、教育的配慮もあって、カスタムメードでいくようだ。
Okay, let's start writing the preprocessing functions. First, we will need a function to convert HTML to plain text. Arguably the best way to do this would be to use the great BeautifulSoupe library, but I would like to avoid adding another dependency to this project, so let's hack a quick & dirty solution using regular expressions.
The following function first drops the <head> section, then converts all <a> tags to the word HYPERLINK, then it gets rid of all HTML tags, leaving only the plain text. For readability, it also replaces multiple newlines with single newlines, and finally it unescapes html entities (such as > or ):
import re
from html import unescape
def html_to_plain_text(html):
text = re.sub('<head.*?>.*?</head>', '', html, flags=re.M | re.S | re.I)
text = re.sub('<a¥s.*?>', ' HYPERLINK ', text, flags=re.M | re.S | re.I)
text = re.sub('<.*?>', '', text, flags=re.M | re.S)
text = re.sub(r' (¥s*¥n)+' , ' ¥n', text, flags=re.M | re.S)
return unescape(text)
プログラムコードを書く（作る）のは面倒だが、ヒトは、もっと面倒な作業を毎日繰り返している。プログラムを組んでしまえば、一瞬で終わってしまう作業かもしれないのに！
html_spam_emails = [email for email in X_train[y_train==1]
if get_email_structure(email) == "text/html"]
sample_html_spam = html_spam_emails[7]
print(sample_html_spam.get_content( ).strip( )[:1000], "...")
----------------------------------------------------------------------
<HTML><HEAD><TITLE></TITLE><META http-equiv="Content-Type" content="text/html; charset=windows-1252"><STYLE>A:link {TEX-DECORATION: none}A:active {TEXT-DECORATION: none}A:visited {TEXT-DECORATION: none}A:hover {COLOR: #0033ff; TEXT-DECORATION: underline}</STYLE><META content="MSHTML 6.00.2713.1100" name="GENERATOR"></HEAD>
<BODY text="#000000" vLink="#0033ff" link="#0033ff" bgColor="#CCCC99"><TABLE borderColor="#660000" cellSpacing="0" cellPadding="0" border="0" width="100%"><TR><TD bgColor="#CCCC99" valign="top" colspan="2" height="27">
<font size="6" face="Arial, Helvetica, sans-serif" color="#660000">
<b>OTC</b></font></TD></TR><TR><TD height="2" bgcolor="#6a694f">
<font size="5" face="Times New Roman, Times, serif" color="#FFFFFF">
<b> Newsletter</b></font></TD><TD height="2" bgcolor="#6a694f"><div align="right"><font color="#FFFFFF">
<b>Discover Tomorrow's Winners </b></font></div></TD></TR><TR><TD height="25" colspan="2" bgcolor="#CCCC99"><table width="100%" border="0"
--------------------------------------------------------------------------
このHTMLをテキストに変換する。
print(html_to_plain_tezt(sample_html_spam.get_content( ))[:1000], "...")
-----------------------------------------------------------------------
OTC
Newsletter
Discover Tomorrow's Winners
For Immediate Release
Cal-Bay (Stock Symbol: CBYI)
Watch for analyst "Strong Buy Recommendations" and several advisory newsletters picking CBYI. CBYI has filed to be traded on the OTCBB, share prices historically INCREASE when companies get listed on this larger trading exchange. CBYI is trading around 25 cents and should skyrocket to $2.66 - $3.25 a share in the near future.
Put CBYI on your watch list, acquire a position TODAY.
REASONS TO INVEST IN CBYI
A profitable company and is on track to beat ALL earnings estimates!
One of the FASTEST growing distributors in environmental & safety equipment instruments.
Excellent management team, several EXCLUSIVE contracts. IMPRESSIVE client list including the U.S. Air Force, Anheuser-Busch, Chevron Refining and Mitsubishi Heavy Industries, GE-Energy & Environmental Research.
RAPIDLY GROWING INDUSTRY
Industry revenues exceed $900 million, estimates indicate that there could be as much as $25 billi ...
-----------------------------------------------------------------------------
HTMLとtextの相互変換のニーズが高いのか、検索すると、ウエブ上で変換できるサイトがいくつも出てきた。上記のHTMLを試してみたら、textの最初の3行が出力された。
-----------------------------------------------------------------------------
つづく
*1:len(X), 1), dtype=bool)
never_5_clf = Never5Classifier( )
cross_val_score(never_5_clf, X_train, y_train_5, cv=3, scoring="accuracy")
array([0.91125, 0.90855, 0.90915])
90％は、５ではないのだから、予測結果を、全て５ではない、としておけば、約90％の正解率になるのは当然である。
簡単な例であるが、accuracyの数値だけで判断してはいけないことが、よくわかる。
MNISTは各ラベルの分布は均等になるように設計されているが、一般的には偏りがある（ラベルの割合が異なる）ので、ベースラインを抑えることは重要である。
Titanicでは、女性だけが生存すると予測すれば、約75％の正解率になるわけだから、75％より高い値になるモデルを構築しなければならないということになる。
Confusion Matrix
分類の正確さを見積もる良い方法として、confusion matrixの考え方がある。
AがAと判定される場合の数とBと判定される場合の数から評価する方法である。
予測の正確さの評価には予測結果が必要だが、テストデータを使ってはいけない。クロスバリデーションを使う。ただし、cross_val_scoreの代わりに、cross_val_predictを使う。
from sklearn.model_selection import cross_val_predict
y_train_pred = cross_val_predict(sgd_clf, X_train, y_train_5, cv=3)
この予測値y_train_predと、ラベルy_train_5から、confusion matricを算出する。
from sklearn.metrics import confusion_matrix
confusion_matrix(y_train_5, y_train_pred)
array([[53057, 1522],
[ 1325, 4096]])
5であることをpositive、5でない数値であることをnegativeと定義していれば、
53057は、5でない数値を5でないと正しく評価したということで、TN:true negative、
4096は、5を正しく5であると評価したということで、TP:true positiveとなる。
1522は、5でない数値を5であると誤って評価したということで、FP:false positive
1325は、5を5でないと誤って評価したということで、FN:false negativeとなる。
perfect classifierによるconfusion matrixは、
y_train_perfect_predictions = y_train_5 # pretand we reached perfection
confusion_matrix(y_train_5, y_train_perfect_predictions)
array([[54579, 0],
[ 0, 5421]])
このconfusion matrixからも評価できるが、precisionやrecallが以下のように定義されて、用いられている。
precision=TP/(TP+FP)
recall=TP/(TP+FN)
このrecallは、sensitivityとかtrue positive rateと呼ばれている。
Precision and Recall
sklearn.metricsには、precision_scoreとrecall score、さらにfi_scoreを計算する関数がある。Fiはprecisionとrecallのharmonic meanである。
from sklearn.metrics import precision_score, recall_score
precision_score(y_train_5, y_train_pred) # == 4096 / (4096 + 1522)
0.729085
recall_score(y_train_5, y_train_pred) # == 4096 / (4096 + 1325)
0.755580
from sklearn.metrics import f1_score
f1_score(y_train_5, y_train_pred) # == 4096 / (4096 + 1522/2 + 1325/2)
0.742096
どのスコアを選ぶかは、分類の目的次第とのことである。
たとえば、最近の例では、新型コロナの保菌者かどうかを判定するPCR検査がある。保菌者は陽性Posotiveと判定される。陽性と判定されても保菌者でない場合は、偽陽性FPである。保菌者が陰性と判定されると偽陰性FNである。偽陰性では、入院できず十分な治療が受けられないという深刻な事態になる可能性がある。こういう判定の場合は、保菌者が陰性と判定される場合の数を直接反映するrecallを装置の指標にするのがよいということになる。
precisionとrecallは、trade-offの関係にあることが多い。
Precision/Recall Trade-off
MNISTのデータベースを用いて、5とそれ以外の数字に分類する、ということを目的として話が進められている。
この分類には、SGDClassifireが用いられている。
分類は、decision functionに基づいて各instanceのスコアを計算し、しきい値より大きいか小さいかで分類している。94ページに、しきい値を中心に12個のinstanceをスコア順に並べ、しきい値が変化すると、precisionとrecallがどのように変化するかを、視覚に訴えながら説明している。
スコアの小さい左側から、８、７、３、９、５、２、５、５、６、５、５、５、と並んでいる。
しきい値が左方の９と５の間にあれば、precisionは、分母がしきい値より右側のすべてのinstanceだから８個、分子は５の個数だから6個となり、6/8で75％である。このときrecallは、分母が５の個数だから6個で、分子も6個となり、100％になる。
しきい値が、最右方の６と5の間にあれば、precisionは3/3で100％、recallは、3/6で50％となる。
このように、分子はしきい値の右側にある５の数だからprecisionもrecallも同じで、precisionの分母は、しきい値の右側にある全てのinstanceの個数だが、recallの分母は、しきい値の位置によらず、すべての５の個数である。
以上のモデルから、precisionとrecallのtrade-offの関係を、少し理解できた気がする。
Scikit-Learnでは、このしきい値を直接、設定することはできないが、アクセスすることはできる。
decision_function( )を使えば、各instanceのスコアがわかる（出力される）。
y_scores = sgd_cls.decision_function([some_digit]) # some_digit=X[0]
y_scores
array([2412.5317])
これがX[0]のスコアで、他の５について調べてみると、スコアが、3939や4742となっていた。
threshold = 0
y_some_digit_pred = (y_scores > threshold)
y_some_digit_pred
array([true])：しきい値を0とすると、5に分類された。
threshold = 8000
y_some_digit_pred = (y_scores > threshold)
y_some_digit_pred
array([False])：しきい値を8000とすると、５に分類されない。
しきい値によって、precisionとrecallが変化することが分かった。
95ページには、しきい値を横軸に、縦軸にprecisionとrecallをプロットした図がある。precisionとrecallが逆の動きをすることがよくわかる。
さらに、96ページには、横軸にrecallを、縦軸にprecisionをプロットした図が示されている。
The ROC Curve
ROC curveもbinary classificationでよく用いられる道具である。
横軸が false positive rate、縦軸が、true positive rate (recall)である。
この曲線ROCの下の面積AUCは、分類器の性能の指標として用いられる。
ROC曲線は次のようにして描かれる。（一部省略）
from sklearn.metrics import roc_curve
fpr, tpr, thresholds = roc_curve(y_train_5, y_scores)
def plot_roc_curve(fpr, tpr, lable=None):
plt.plot(fpr, tpr, linewidth=2, lable=lable)
plt.plot([0, 1], [0, 1], 'k--') # dashed diagonal
plot_roc_curve(fpr, tpr)
plt.show( )
ランダムな場合に直線となり、性能が高いほど、左上隅に近づく。
曲線下の面積AUCは、ランダムが0.5、完全分類が1.0となる。
from sklearn.metrics import roc_auc_score
roc_auc_score(y_train_5, y_scores)
0.96117788
96ページのPR曲線（Precision/Recall curve）とROC曲線の使い分け。
経験則として、precision/recall:PRは、positive classが少ないとき（5の数字とそれ以外の数字の分類のように）、または、false negativeよりfalse positiveの方が重要な場合に用いる。それ以外の場合はROCを用いる。
さて、RandomForestClassifierモデルを訓練して、ROC曲線とROC AUC scoreをSGDClassifierと比較してみよう。
RandamForestClassifierはdecision_function( ) methodを持っていないが、そのかわり、predict_proba( ) methodを持っている。
from sklearn.ensamble import RandomForestClassifier
forest_clf = RandomForestClassifier(randomstate=42)
y_probas_forest = cross_val_predict(forest_clf, X_train, y_train_5, cv=3,
method="predict_proba")
y_probas_forest = y_probas_forest[ : , 1] # score = proba of positive class
fpr_forest, tpr_forest, threshold_forest = roc_curve(y_train_5, y_scores_forest)
plt.plot(fpr, tpr, "b:", lable="SGD")
plot_roc_curve(fpr_forest, tpr_forest, "Random Forest")
plt.legend(loc="lower right")
plt.show( )
RandomForestClassifierのROC curcurveはSGDClassifierよりもさらに左上隅にある。
ROC AUCを計算してみると、
roc_auc_score(y_train_5, y_scores_forest)
0.99834367
さらに、precisionとrecallは、
y_train_pred_forest = cross_val_predict(forest_clf, X_train, y_train_5, cv=3)
precision_score(y_train_5, y_train_pred_forest)
0.9905083
recall_score(y_train_5, y_train_pred_forest)
0.8662608
悪くない結果だ！
ここまで、binary classifierをみてきた。適切な評価法を選ぶ、cross-validationでclassifierを評価する、precision/recallを使い分ける、ROC curveとROC AUC scoreをモデルを比較するために用いる。次は、５以外の数字を検出しよう。
Multiclass Classification
SGD Classifiers, Random Forest Classifiers, naive Bayesなどは、もともと、multiple classes に対応している。
Logistic RegressionやSupport Vector Machineなどは、binary classifierである。
しかしながら、binary classifierも、multiple binary classifiersとして、multiple classificationは可能である。
one-versus-the-rest:OvR, one-versus-all, one-versus-one:OvO
binary classifierをmulticlass classifierとして使うにはいろいろな手法があるが、
Scikit-Learnは、multi-class classificationにbinary classification algirithmを使っても、algorithmによって、自動的に、OvRとOvOを使い分けている。
実際に、SVMをMNISTの多値分類（0-9）に適用してみよう。
svm_clf.fit(X_train, y_train) # y_train, not y_train_5
svm_clf.predict([some_digit]) # some_digit : X[0] : [5]
array([5], dtype=unit8)
y_trainを使っているので、0 - 9 original target classesである。
内部では、OvO strategyが用いられ、45のbinary classifiersが訓練され、各instanceに対してdecision scoreが計算される。
some_digit_scores = svm_clf.decision_function([some_digit]) # some_digit : [5]
some_digit_scores
array( [ [ 2.92, 7.02 3.94, 0.90, 5.97, 9.50, 1.91, 8.03, -0.13, 4.94] ] )
5のスコア：9.50が最も高くなっていることが確認できる。
念のため、
np.argmax(some_digit_scores)
5
svm_clf.classes_
array( [ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9], dtype=unit8)
svm_clf.classes_[5]
5
Scikit-LearnはOvOかOvRかを自動的に選んでいるが、指定もできる。
from sklearn.multiclass import OneVsRestClassifier
ovr_clf = OneVsRestClassifier(SVC(
