PyData.Tokyoでは毎月開催している中上級者向けの勉強会に加え、初心者の育成を目的としたチュートリアルイベントを開催します。今回のイベントでは下記の項目にフォーカスします。
このチュートリアルでは実際のデータを使ったコーディングを行うことで実践力をつけることを目的とします。扱う事例はタイタニックの乗客データを使った生存者推定モデルの生成です。乗客の年齢、性別、その他の情報を機械学習アルゴリズムに学習させることで、初心者でも80％に近い精度で生存者を当てることができるようになります。
チュートリアル第二部では、Pythonの機械学習ライブラリscikit-learnを使って、次の2つの点について学びます。
タイタニックの乗客データ： Titanic: Machine Learning from Disaster
※データのダウンロードには、Kaggleのアカウントが必要です。
PyData.Tokyo オーガナイザー 田中 秀樹（@atelierhide）
シリコンバレーでPython×Dataの魅力に出会う。その後、ディープラーニングに興味を持ち、PyCon JP 2014に登壇したことがきっかけとなりPyData.Tokyoをスタート。カメラレンズの光学設計エンジニアをする傍ら、画像認識を用いた火星および太陽系惑星表面の構造物探索を行うMarsface Project（@marsfaceproject）に参加。
In [1]:
from IPython.display import Image
Image(url='http://graphics8.nytimes.com/images/section/learning/general/onthisday/big/0415_big.gif')
Out[1]:
最初に、必要なライブラリをインポートしましょう。
In [2]:
%matplotlib inline
In [3]:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
from sklearn.metrics import accuracy_score, classification_report, confusion_matrix
from sklearn.model_selection import train_test_split, cross_val_score, KFold, GridSearchCV
from sklearn.preprocessing import LabelEncoder
from sklearn.linear_model import LogisticRegression
from sklearn.svm import SVC
from sklearn.tree import DecisionTreeClassifier, export_graphviz
from IPython.display import Image
# Pandasの設定をします
pd.set_option('chained_assignment', None)
# matplotlibのスタイルを指定します。これでグラフが少しかっこよくなります。
plt.style.use('ggplot')
plt.rc('xtick.major', size=0)
plt.rc('ytick.major', size=0)
PandasのDataFrameに2つのデータを読込みます。train.csvは乗客の生存情報が付いているトレーニングデータ（教師データ）です。test.csvは生存情報を推定してKaggleに投稿するためのテストデータのため、生存情報が付いていません。
In [4]:
df_train = pd.read_csv('data/train.csv')
df_test = pd.read_csv('data/test.csv')
2つのデータを確認してみましょう。df_trainにのみ生存情報（Survived）があるのがわかります。
In [5]:
df_train.tail()
Out[5]:
In [6]:
df_test.tail()
Out[6]:
前半のチュートリアルのデータ解析で、生存確率は男性より女性の方が高いことが分かりました。先ず、最も単純なモデルとして、性別により生存者を予測するモデル（ジェンダーモデル）を考えてみましょう。
In [7]:
x = df_train['Sex']
y = df_train['Survived']
ジェンダーモデルで生存者を推定します。ジェンダーモデルは、女性は全員が生存（0）、男性は全員が死亡（1）と仮定するモデルです。y_predのpredはpredictionの略です。pandasのmapを使って計算してみましょう。
In [8]:
y_pred = x.map({'female': 1, 'male': 0}).astype(int)
推定したデータを評価します。最初に正解率（Accuracy）を求めましょう。accuracy_scoreで計算します。
In [9]:
print('Accuracy: {:.3f}'.format(accuracy_score(y, y_pred)))
78.7%の正解率が得られました。データを理解して仮説を立てることで、単純なモデルでも高い正解率が得られることが分かります。Kaggleでは、コンペによって使われている指標が異なりますが、タイタニックのコンペでは正解率が指標となっています。
他の指標もscikit-learnで簡単に計算出来ます。Precision、Recall、F1-scoreをclassification_reportで計算してみましょう。
In [10]:
print(classification_report(y, y_pred))
混同行列（Confusion Matrix）は、推定結果を理解するのにとても便利です。scikit-learnのconfusion_matrixで計算し、結果をmatplotlibで可視化してみましょう。
In [11]:
cm = confusion_matrix(y, y_pred)
print(cm)
In [12]:
def plot_confusion_matrix(cm):
fig, ax = plt.subplots()
im = ax.imshow(cm, interpolation='nearest', cmap=plt.cm.Blues)
ax.set_title('Confusion Matrix')
fig.colorbar(im)
target_names = ['not survived', 'survived']
tick_marks = np.arange(len(target_names))
ax.set_xticks(tick_marks)
ax.set_xticklabels(target_names, rotation=45)
ax.set_yticks(tick_marks)
ax.set_yticklabels(target_names)
ax.set_ylabel('True label')
ax.set_xlabel('Predicted label')
fig.tight_layout()
plot_confusion_matrix(cm)
In [13]:
x_test = df_test['Sex']
y_test_pred = x_test.map({'female': 1, 'male': 0}).astype(int)
推定した結果を、Kaggleに投稿するためのCSVファイルを作成します。CSVファイルに記載する必要のあるデータはPassengerIdとSurvived（生存者の推定値）です。pandasで投稿データ用のDataFrameを作成し、to_csvを使ってCSV形式で保存します。
In [14]:
df_kaggle = pd.DataFrame({'PassengerId': df_test['PassengerId'], 'Survived':np.array(y_test_pred)})
df_kaggle.to_csv('kaggle_gendermodel.csv', index=False)
In [15]:
df_kaggle.head()
Out[15]:
作成したkaggle_gendermodel.csvをKaggleに投稿し、スコアと順位を確認してみましょう！これで皆さんもKagglerです！
scikit-learnに実装されている機械学習のアルゴリズムを使うことを学びます。先ずは最も基本的な線形モデルから始めましょう。
In [16]:
X = df_train[['Age', 'Pclass', 'Sex']]
y = df_train['Survived']
特徴量のデータフレームを確認します。
In [17]:
X.tail()
Out[17]:
年齢に欠損値があります。教師データのサイズが十分に大きければ、欠損値を使わなくても問題ありません。今回は教師データがあまり大きくないため、欠損値を埋めて使います。チュートリアル第一部では、欠損値を埋める手法をいくつか紹介しましたが、今回は全体の平均値を使うことにします。
In [18]:
X['AgeFill'] = X['Age'].fillna(X['Age'].mean())
X = X.drop(['Age'], axis=1)
また、性別（Sex）はmaleとfemaleという値が入っていますが、scikit-learnでは、このようなカテゴリー情報を扱うことが出来ません。そのため、female、maleを数値に変換する必要があります。femaleを0、maleを1とし、新しくGenderを作成します。
In [19]:
X['Gender'] = X['Sex'].map({'female': 0, 'male': 1}).astype(int)
In [20]:
X.tail()
Out[20]:
次に、女性（Gender=0）で且つ、乗船クラスのランクが高い（Pclass=1）ほど、生存率が高いという仮説を表す新しい特徴量（Pclass_Gender）を作成します。Pclass_Genderは値が小さいほど生存率が高いことになります。
In [21]:
X['Pclass_Gender'] = X['Pclass'] + X['Gender']
In [22]:
X.tail()
Out[22]:
今回は特徴量としてPclass_GenderとAgeの2つを使います。不要になった特徴量は、dropで削除します。
In [23]:
X = X.drop(['Pclass', 'Sex', 'Gender'], axis=1)
In [24]:
X.head()
Out[24]:
データを可視化して「年齢が若く、女性で且つ、乗船クラスのランクが高いほど、生存率が高い」という仮説が正しいか確認してみましょう。横軸が年齢、縦軸がPclass_Genderを表します。
In [25]:
np.random.seed = 0
xmin, xmax = -5, 85
ymin, ymax = 0.5, 4.5
index_survived = y[y==0].index
index_notsurvived = y[y==1].index
fig, ax = plt.subplots()
cm = plt.cm.RdBu
cm_bright = ListedColormap(['#FF0000', '#0000FF'])
sc = ax.scatter(X.loc[index_survived, 'AgeFill'],
X.loc[index_survived, 'Pclass_Gender']+(np.random.rand(len(index_survived))-0.5)*0.1,
color='r', label='Not Survived', alpha=0.3)
sc = ax.scatter(X.loc[index_notsurvived, 'AgeFill'],
X.loc[index_notsurvived, 'Pclass_Gender']+(np.random.rand(len(index_notsurvived))-0.5)*0.1,
color='b', label='Survived', alpha=0.3)
ax.set_xlabel('AgeFill')
ax.set_ylabel('Pclass_Gender')
ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)
ax.legend(bbox_to_anchor=(1.4, 1.03))
plt.show()
いかがでしょうか？仮説は正しく、グラフの左下で生存者が多くなっていますね。
機械学習では、学習にデータをすべて使ってしまうと、モデルを正しく評価出来ません。そこで、データを学習用と評価用の2つに分割します。分割はscikit-learnのtrain_test_splitで簡単に行うことが出来ます。ここでは、データの80%を学習用、20%を評価用として分割します。x_val、y_valのvalはvalidationの略です。
In [26]:
X_train, X_val, y_train, y_val = train_test_split(X, y, train_size=0.8, random_state=1)
In [27]:
X_train
Out[27]:
In [28]:
print('Num of Training Samples: {}'.format(len(X_train)))
print('Num of Validation Samples: {}'.format(len(X_val)))
線形モデルであるロジスティック回帰（LogisticRegression）を使います。clfは分類器を意味するclassifierの略です。
In [29]:
clf = LogisticRegression()
先ほど分割した学習用データを使います。
In [30]:
clf.fit(X_train, y_train)
Out[30]:
これで、学習は完了です。データ数が少ないため、あっという間に終わります。
次に生存者の推定をしますが、こちらも簡単に行えます。先ほど分割した学習用データと、評価用データの両方について推定します。
In [31]:
y_train_pred = clf.predict(X_train)
y_val_pred = clf.predict(X_val)
結果を評価してみましょう。
In [32]:
print('Accuracy on Training Set: {:.3f}'.format(accuracy_score(y_train, y_train_pred)))
print('Accuracy on Validation Set: {:.3f}'.format(accuracy_score(y_val, y_val_pred)))
In [33]:
cm = confusion_matrix(y_val, y_val_pred)
print(cm)
In [34]:
plot_confusion_matrix(cm)
ロジスティック回帰はどのような役割を確認しましょう。matplotlibで可視化します。
In [35]:
h = 0.02
xmin, xmax = -5, 85
ymin, ymax = 0.5, 4.5
xx, yy = np.meshgrid(np.arange(xmin, xmax, h), np.arange(ymin, ymax, h))
Z = clf.predict_proba(np.c_[xx.ravel(), yy.ravel()])[:, 1]
Z = Z.reshape(xx.shape)
fig, ax = plt.subplots()
levels = np.linspace(0, 1.0, 5)
cm = plt.cm.RdBu
cm_bright = ListedColormap(['#FF0000', '#0000FF'])
contour = ax.contourf(xx, yy, Z, cmap=cm, levels=levels, alpha=0.8)
ax.scatter(X_train.iloc[:, 0], X_train.iloc[:, 1]+(np.random.rand(len(X_train))-0.5)*0.1, c=y_train, cmap=cm_bright)
ax.scatter(X_val.iloc[:, 0], X_val.iloc[:, 1]+(np.random.rand(len(X_val))-0.5)*0.1, c=y_val, cmap=cm_bright, alpha=0.5)
ax.set_xlabel('AgeFill')
ax.set_ylabel('Pclass_Gender')
ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)
fig.colorbar(contour)
x1 = xmin
x2 = xmax
y1 = -1*(clf.intercept_[0]+clf.coef_[0][0]*xmin)/clf.coef_[0][1]
y2 = -1*(clf.intercept_[0]+clf.coef_[0][0]*xmax)/clf.coef_[0][1]
ax.plot([x1, x2] ,[y1, y2], 'k--')
plt.show()
ロジスティック回帰は与えられた特徴量に基づいて、乗客が生存したか、生存しなかったかの境界（グラフの点線）を判断しています。これの境界はHyperplane（超平面）または、Decision Boundary（決定境界）と呼ばれます。機械学習の分類問題の目的は、この境界を求めることとも言えます。アルゴリズムによって、この境界の求め方が異なり、結果も異なります。機械学習の様々な分野で広く使われいているSVM（Support Vector Machines）と比較してみましょう。アルゴリズムの詳細の説明は省略します。
In [36]:
clf_log = LogisticRegression()
clf_svc_lin = SVC(kernel='linear', probability=True)
clf_svc_rbf = SVC(kernel='rbf', probability=True)
titles = ['Logistic Regression', 'SVC with Linear Kernel', 'SVC with RBF Kernel',]
h = 0.02
xmin, xmax = -5, 85
ymin, ymax = 0.5, 4.5
xx, yy = np.meshgrid(np.arange(xmin, xmax, h), np.arange(ymin, ymax, h))
fig, axes = plt.subplots(1, 3, figsize=(12,4))
levels = np.linspace(0, 1.0, 5)
cm = plt.cm.RdBu
cm_bright = ListedColormap(['#FF0000', '#0000FF'])
for i, clf in enumerate((clf_log, clf_svc_lin, clf_svc_rbf)):
clf.fit(X, y)
Z = clf.predict(np.c_[xx.ravel(), yy.ravel()])
Z = Z.reshape(xx.shape)
axes[i].contourf(xx, yy, Z, cmap=cm, levels=levels, alpha=0.8)
axes[i].scatter(X.iloc[:, 0], X.iloc[:, 1], c=y, cmap=cm_bright)
axes[i].set_title(titles[i])
axes[i].set_xlabel('AgeFill')
axes[i].set_ylabel('Pclass_Gender')
axes[i].set_xlim(xmin, xmax)
axes[i].set_ylim(ymin, ymax)
fig.tight_layout()
In [37]:
clf = SVC(kernel='rbf', probability=True)
clf.fit(X_train, y_train)
y_train_pred = clf.predict(X_train)
y_val_pred = clf.predict(X_val)
print('Accuracy on Training Set: {:.3f}'.format(accuracy_score(y_train, y_train_pred)))
print('Accuracy on Validation Set: {:.3f}'.format(accuracy_score(y_val, y_val_pred)))
In [38]:
X_train, X_val, y_train, y_val = train_test_split(X, y, train_size=0.8, random_state=33)
clf = LogisticRegression()
clf.fit(X_train, y_train)
y_train_pred = clf.predict(X_train)
y_val_pred = clf.predict(X_val)
print('Accuracy on Training Set: {:.3f}'.format(accuracy_score(y_train, y_train_pred)))
print('Accuracy on Test Set: {:.3f}'.format(accuracy_score(y_val, y_val_pred)))
どの部分を教師データにするかによって結果は異なります。この課題を解決する方法として交差検証（クロスバリデーション、Cross-validation）という手法があります。ここでは、K-分割交差検証（K-fold cross-validation）を使いましょう。K-分割交差検証は、データをK個に分割し、そのうちK-1個のデータセットを学習に、K個のデータを訓練に用いるということをK回繰り返し、得られた結果の平均を得るというものです。例えば、5-fold cross-validationの場合、5つのデータセットを作成します。各データセットには20%のサンプルが含まれていることになりますが、これを利用し、80%のサンプルで学習、20%のサンプルで評価するということを5回繰り返します。
In [39]:
Image(url='http://scott.fortmann-roe.com/docs/docs/MeasuringError/crossvalidation.png')
Out[39]:
scikit-learnでもcross_validationとして実装されています。K-分割交差検証を関数として定義し、実行してみましょう。
In [40]:
def cross_val(clf, X, y, K=5, random_state=0):
cv = KFold(K, shuffle=True, random_state=random_state)
scores = cross_val_score(clf, X, y, cv=cv)
return scores
In [41]:
cv = KFold(5, shuffle=True, random_state=0)
cv
Out[41]:
In [42]:
clf = LogisticRegression()
scores = cross_val(clf, X, y)
print('Scores:', scores)
print('Mean Score: {0:.3f} (+/-{1:.3f})'.format(scores.mean(), scores.std()*2))
In [43]:
X = df_train[['Age', 'Pclass', 'Sex', 'SibSp', 'Parch', 'Embarked']]
y = df_train['Survived']
X_test = df_test[['Age', 'Pclass', 'Sex', 'SibSp', 'Parch', 'Embarked']]
In [44]:
X.tail()
Out[44]:
In [45]:
X['AgeFill'] = X['Age'].fillna(X['Age'].mean())
X_test['AgeFill'] = X_test['Age'].fillna(X['Age'].mean())
X = X.drop(['Age'], axis=1)
X_test = X_test.drop(['Age'], axis=1)
性別の数値化にはscikit-learnのLabelEncoderを使うことも出来ます。
In [46]:
le = LabelEncoder()
le.fit(X['Sex'])
X['Gender'] = le.transform(X['Sex'])
X_test['Gender'] = le.transform(X_test['Sex'])
classes = {gender: i for (i, gender) in enumerate(le.classes_)}
print(classes)
In [47]:
X.tail()
Out[47]:
性別（Sex）と同様に乗船地（Embarked）もそのままでは使えないため、数値化する必要がありますが、対象となるのはS、C、Qの3種類です。このような場合は、One-hot Encording、またはOne-of-K Encordingという手法を使って、新たな特徴量を作ります。pandasのget_dummiesを使います。
In [48]:
X = X.join(pd.get_dummies(X['Embarked'], prefix='Embarked'))
X_test = X_test.join(pd.get_dummies(X['Embarked'], prefix='Embarked'))
In [49]:
X.tail()
Out[49]:
不要な特徴量は削除しましょう。
In [50]:
X = X.drop(['Sex', 'Embarked'], axis=1)
X_test = X_test.drop(['Sex', 'Embarked'], axis=1)
ロジスティック回帰＋交差検証で評価します。
In [51]:
clf = LogisticRegression()
scores = cross_val(clf, X, y)
print('Scores:', scores)
print('Mean Score: {0:.3f} (+/-{1:.3f})'.format(scores.mean(), scores.std()*2))
スコアが改善しました！
決定木は、機械学習の手法の中でも非常によく用いられるものの一つです。分類を決定づけた要因を木構造で説明することが出来るため、非常に分かりやすいという特徴があります。
In [52]:
clf = DecisionTreeClassifier(criterion='entropy', max_depth=2, min_samples_leaf=2)
scores = cross_val(clf, X, y, 5)
print('Scores:', scores)
print('Mean Score: {0:.3f} (+/-{1:.3f})'.format(scores.mean(), scores.std()*2))
In [53]:
Image(url='https://raw.githubusercontent.com/PyDataTokyo/pydata-tokyo-tutorial-1/master/images/titanic_decision_tree.png')
Out[53]:
In [54]:
clf = DecisionTreeClassifier(criterion='entropy', max_depth=3, min_samples_leaf=2)
scores = cross_val(clf, X, y, 5)
print('Scores:', scores)
print('Mean Score: {0:.3f} (+/-{1:.3f})'.format(scores.mean(), scores.std()*2))
In [55]:
clf = DecisionTreeClassifier(criterion='entropy', max_depth=2, min_samples_leaf=2)
param_grid = {'max_depth': [2, 3, 4, 5], 'min_samples_leaf': [2, 3, 4, 5]}
cv = KFold(5, shuffle=True, random_state=0)
grid_search = GridSearchCV(clf, param_grid, cv=cv, n_jobs=-1, verbose=1,return_train_score=True)
grid_search.fit(X, y)
Out[55]:
ベストなスコアとパラメータの組合せを確認します。
In [56]:
print('Scores: {:.3f}'.format(grid_search.best_score_))
print('Best Parameter Choice:', grid_search.best_params_)
全ての結果を確認することも出来ます。
In [57]:
grid_search.cv_results_
Out[57]:
In [58]:
grid_search.cv_results_['mean_test_score']
Out[58]:
In [59]:
scores = grid_search.cv_results_['mean_test_score'].reshape(4, 4)
fig, ax = plt.subplots()
cm = plt.cm.Blues
mat = ax.matshow(scores, cmap=cm)
ax.set_xlabel('min_samples_leaf')
ax.set_ylabel('max_depth')
ax.set_xticklabels(['']+param_grid['min_samples_leaf'])
ax.set_yticklabels(['']+param_grid['max_depth'])
fig.colorbar(mat)
plt.show()
ベストなパラメータの組合せで推定を行います。
In [60]:
y_test_pred = grid_search.predict(X_test)
Kaggleに投稿するためのCSVファイルを作成し、結果を確認してみましょう。
In [61]:
df_kaggle = pd.DataFrame({'PassengerId': df_test['PassengerId'], 'Survived':np.array(y_test_pred)})
df_kaggle.to_csv('kaggle_decisiontree.csv', index=False)
チュートリアル第二部はこれで終了です。ここで学んだことを活かして、さらに高いスコアを目指してください！