機械学習で二値分類

# 数値計算やデータフレーム操作に関するライブラリをインポートする
import numpy as np
import pandas as pd
import scipy as sp
from scipy import stats

# URL によるリソースへのアクセスを提供するライブラリをインポートする。
# import urllib # Python 2 の場合
import urllib.request # Python 3 の場合

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

# 機械学習関連のライブラリ群

from sklearn.cross_validation import train_test_split # 訓練データとテストデータに分割
from sklearn.metrics import confusion_matrix # 混合行列

from sklearn.decomposition import PCA #主成分分析
from sklearn.linear_model import LogisticRegression # ロジスティック回帰
from sklearn.neighbors import KNeighborsClassifier # K近傍法
from sklearn.svm import SVC # サポートベクターマシン
from sklearn.tree import DecisionTreeClassifier # 決定木
from sklearn.ensemble import RandomForestClassifier # ランダムフォレスト
from sklearn.ensemble import AdaBoostClassifier # AdaBoost
from sklearn.naive_bayes import GaussianNB # ナイーブ・ベイズ
from sklearn.lda import LDA # 線形判別分析
from sklearn.qda import QDA # 二次判別分析

1. 「ワインの品質」データ読み込み

データは UC Irvine Machine Learning Repository から取得したものを少し改変しました。

• 赤ワイン https://raw.githubusercontent.com/chemo-wakate/tutorial-6th/master/beginner/data/winequality-red.txt

• 白ワイン https://raw.githubusercontent.com/chemo-wakate/tutorial-6th/master/beginner/data/winequality-white.txt

　詳細

1. fixed acidity : 不揮発酸濃度（ほぼ酒石酸濃度）
2. volatile acidity : 揮発酸濃度（ほぼ酢酸濃度）
3. citric acid : クエン酸濃度
4. residual sugar : 残存糖濃度
5. chlorides : 塩化物濃度
6. free sulfur dioxide : 遊離亜硫酸濃度
7. total sulfur dioxide : 亜硫酸濃度
8. density : 密度
9. pH : pH
10. sulphates : 硫酸塩濃度
11. alcohol : アルコール度数
12. quality (score between 0 and 10) : 0-10 の値で示される品質のスコア

# ウェブ上のリソースを指定する
url = 'https://raw.githubusercontent.com/chemo-wakate/tutorial-6th/master/beginner/data/winequality-red.txt'
# 指定したURLからリソースをダウンロードし、名前をつける。
# urllib.urlretrieve(url, 'winequality-red.csv') # Python 2 の場合
urllib.request.urlretrieve(url, 'winequality-red.txt') # Python 3 の場合

# データの読み込み
df1 = pd.read_csv('winequality-red.txt', sep='\t', index_col=0) 

df1.head() # 先頭５行だけ表示

2. ２群に分ける

ここでは、ワインの品質を「６未満（よくない）」と「６以上（よい）」の２群に分けてから、機械学習を用いて、pH や volatile acidity などの変数から品質を予測してみましょう。まずは、２群に分けることから始めます。

簡単な例で説明

データを２群に分けるにあたって、pandasの操作が少し分かりにくいので、簡単な例を用いて説明します。

# 簡単な例
toy_data = pd.DataFrame([[1,4,7,10,13,16],[2,5,8,11,14,27],[3,6,9,12,15,17],[21,24,27,20,23,26]],
                   index = ['i1','i2','i3', 'i4'],
                   columns = list("abcdef"))

toy_data # 中身の確認

# F列の値が 20 未満の列だけを抜き出す
toy_data[toy_data['f'] < 20]

# F列の値が 20 以上の列だけを抜き出す
toy_data[toy_data['f'] >= 20]

# F列の値が 20 以上の列だけを抜き出して、そのB列を得る
pd.DataFrame(toy_data[toy_data['f'] >= 20]['b'])

# classという名の列を作り、F列の値が 20 未満なら 0 を、そうでなければ 1 を入れる
toy_data['class'] = [0 if i < 20 else 1 for i in toy_data['f'].tolist()]

toy_data # 中身を確認

実データに戻ります

以下、quality が６未満のワインと６以上のワインに分け、どのように違うのか調べてみましょう。

# quality が 6 未満の行を抜き出して、先頭５行を表示する
df1[df1['quality'] < 6].head()

# quality が 6 以上の行を抜き出して、先頭５行を表示する
df1[df1['quality'] >= 6].head()

fig, ax = plt.subplots(1, 1)

# quality が 6 未満の行を抜き出して、x軸を volatile acidity 、 y軸を alcohol として青色の丸を散布する
df1[df1['quality'] <  6].plot(kind='scatter', x=u'volatile acidity', y=u'alcohol', ax=ax, 
                              c='blue', alpha=0.5)

# quality が 6 以上の行を抜き出して、x軸を volatile acidity 、 y軸を alcohol として赤色の丸を散布する
df1[df1['quality'] >= 6].plot(kind='scatter', x=u'volatile acidity', y=u'alcohol', ax=ax, 
                              c='red', alpha=0.5, grid=True, figsize=(5,5))
plt.show()

# quality が 6 未満のものを青色、6以上のものを赤色に彩色して volatile acidity の分布を描画
df1[df1['quality'] <  6]['volatile acidity'].hist(figsize=(3, 3), bins=20, alpha=0.5, color='blue')
df1[df1['quality'] >= 6]['volatile acidity'].hist(figsize=(3, 3), bins=20, alpha=0.5, color='red')

上図のように、qualityが６未満のワインと６以上のワインは volatile acidity の分布が異なるように見えます。

同様に、qualityが６未満のワインと６以上のワインでは pH の分布が異なるか調べてみましょう。

# quality が 6 未満のものを青色、6以上のものを赤色に彩色して pH の分布を描画
df1[df1['quality'] <  6]['pH'].hist(figsize=(3, 3), bins=20, alpha=0.5, color='blue')
df1[df1['quality'] >= 6]['pH'].hist(figsize=(3, 3), bins=20, alpha=0.5, color='red')

分類を表す列を追加する

quality が 6 未満のワインを「0」、6以上のワインを「1」とした class 列を追加しましょう。

df1['class'] = [0 if i <= 5 else 1 for i in df1['quality'].tolist()]

df1.head() # 先頭５行を表示

class 列が 0 なら青色、1 なら赤色に彩色します。

# それぞれに与える色を決める。
color_codes = {0:'#0000FF', 1:'#FF0000'}
colors = [color_codes[x] for x in df1['class'].tolist()]

その彩色で散布図行列を描きましょう。

plotting.scatter_matrix(df1.dropna(axis=1)[df1.columns[:11]], figsize=(20, 20), color=colors, alpha=0.5) 
plt.show()

上図から、各変数と、quality の良し悪しとの関係がボンヤリとつかめてきたのではないでしょうか。続いて主成分分析をしてみます。

dfs = df1.apply(lambda x: (x-x.mean())/x.std(), axis=0).fillna(0) # データの正規化

pca = PCA()
pca.fit(dfs.iloc[:, :11])
# データを主成分空間に写像 = 次元圧縮
feature = pca.transform(dfs.iloc[:, :11])
#plt.figure(figsize=(6, 6))
plt.scatter(feature[:, 0], feature[:, 1], alpha=0.5, color=colors)
plt.title("Principal Component Analysis")
plt.xlabel("The first principal component")
plt.ylabel("The second principal component")
plt.grid()
plt.show()

分かったような分からないような結果ですね。quality の良し悪しを分類・予測するのは簡単ではなさそうです。

3. 説明変数と目的変数に分ける

ここまでで、ワインの品質を２群に分けました。次は、目的変数（ここでは、品質）と説明変数（それ以外の変数）に分けましょう。

X = dfs.iloc[:, :11] # 説明変数
y = df1.iloc[:, 12] # 目的変数

X.head() # 先頭５行を表示して確認

pd.DataFrame(y).T # 目的変数を確認。縦に長いと見にくいので転置して表示。

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

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

• 訓練データ（全データの６０％）
  • X_train : 訓練データの説明変数
  • y_train : 訓練データの目的変数
• テストデータ（全データの４０％）
  • X_test : テストデータの説明変数
  • y_test : テストデータの目的変数

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

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=.4) # 訓練データ・テストデータへのランダムな分割

X_train.head() # 先頭５行を表示して確認

pd.DataFrame(y_train).T # 縦に長いと見にくいので転置して表示。

5. ロジスティック回帰

機械学習モデルとして有名なものの一つとして、ロジスティック回帰があります。線形回帰分析が量的変数を予測するのに対して、ロジスティック回帰分析は発生確率を予測する手法です。基本的な考え方は線形回帰分析と同じなのですが、予測結果が 0 から 1 の間を取るように、数式やその前提に改良が加えられています。

clf = LogisticRegression() #モデルの生成
clf.fit(X_train, y_train) #学習

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

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

正解率の数字を出すだけなら以上でおしまいですが、具体的な予測結果を確認したい場合は次のようにします。

y_predict = clf.predict(X_test)

pd.DataFrame(y_predict).T

# 予測結果と、正解（本当の答え）がどのくらい合っていたかを表す混合行列
pd.DataFrame(confusion_matrix(y_predict, y_test), index=['predicted 0', 'predicted 1'], columns=['real 0', 'real 1'])

6. いろんな機械学習手法を比較する

scikit-learn が用意している機械学習手法（分類器）はロジスティック回帰だけではありません。有名なものは SVM （サポートベクターマシン）などがあります。いろいろ試して、ベストなものを選択してみましょう。

まず、いろんな分類器を classifiers という名のリストの中に収納します。

names = ["Logistic Regression", "Nearest Neighbors", 
         "Linear SVM", "Polynomial SVM", "RBF SVM", "Sigmoid SVM", 
         "Decision Tree","Random Forest", "AdaBoost", "Naive Bayes", 
         "Linear Discriminant Analysis","Quadratic Discriminant Analysis"]

classifiers = [
    LogisticRegression(),
    KNeighborsClassifier(),
    SVC(kernel="linear"),
    SVC(kernel="poly"),
    SVC(kernel="rbf"),
    SVC(kernel="sigmoid"),
    DecisionTreeClassifier(),
    RandomForestClassifier(),
    AdaBoostClassifier(),
    GaussianNB(),
    LDA(),
    QDA()]

さきほど作成した教師データを使って、これらの分類器で順番に予測して、正解率（train）と正解率（test）を計算してみましょう。

result = []
for name, clf in zip(names, classifiers): # 指定した複数の分類機を順番に呼び出す
    clf.fit(X_train, y_train) # 学習
    score1 = clf.score(X_train, y_train) # 正解率（train）の算出
    score2 = clf.score(X_test, y_test) # 正解率（test）の算出
    result.append([score1, score2]) # 結果の格納

# test の正解率の大きい順に並べる
df_result = pd.DataFrame(result, columns=['train', 'test'], index=names).sort('test', ascending=False)
# もし上のコードが動かない場合、以下のコーを試してみてください。
#df_result = pd.DataFrame(result, columns=['train', 'test'], index=names).sort_values('test', ascending=False)

df_result # 結果の確認

# 棒グラフの描画
df_result.plot(kind='bar', alpha=0.5, grid=True)

訓練データの作成はランダムに行なうので、作成のたびに正解率の数字は変わります。場合によっては、分類器の順序が前後することもあります。それでは適切な性能評価がしにくいので、教師データを何度も作り直して正解率を計算してみましょう。

result = []
for trial in range(20): # 20 回繰り返す
    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=.4) # 訓練データ・テストデータの生成
    for name, clf in zip(names, classifiers): # 指定した複数の分類機を順番に呼び出す
        clf.fit(X_train, y_train) # 学習
        score1 = clf.score(X_train, y_train) # 正解率（train）の算出
        score2 = clf.score(X_test, y_test) # 正解率（test）の算出
        result.append([name, score1, score2]) # 結果の格納

df_result = pd.DataFrame(result, columns=['classifier', 'train', 'test']) # 今回はまだ並べ替えはしない

df_result # 結果の確認。同じ分類器の結果が複数回登場していることに注意。

# 分類器 (classifier) 毎にグループ化して正解率の平均を計算し、test の正解率の平均の大きい順に並べる
df_result_mean = df_result.groupby('classifier').mean().sort('test', ascending=False)
# もし上のコードが動かない場合、以下のコーを試してみてください。
#df_result_mean = df_result.groupby('classifier').mean().sort_values('test', ascending=False)

df_result_mean # 結果の確認

# エラーバーの表示に用いる目的で、標準偏差を計算する
errors = df_result.groupby('classifier').std()

errors # 結果の確認

# 平均値と標準偏差を用いて棒グラフを描画
df_result_mean.plot(kind='bar', alpha=0.5, grid=True, yerr=errors)

以上、様々な分類器を用いて、ワインの品質の善し悪しを予測しました。それぞれの分類器にはそれぞれのパラメーターがありますが、上の例では全てデフォルト値を使っています。上手にパラメーターをチューニングすれば、もっと良い予測性能が出せるかもしれません。ですが今回はここまでとさせていただきます。興味があったらぜひ調べてみてください。