Python 機械学習 Scikit-learnによる決定木学習の実践

今回は決定木によるIrisの花のデータセットの分類を実践します。サンプルプログラムをまずは動かしてみることを念頭にしています。

※Pythonの導入方法http://costudyinfodatabase.nagoya/2017/01/05/%e3%81%82-2/

※Iris のデータセットとは
Irisのデータセットは「Iris setosa」「 Iris versicolor」「 Iris virginica 」の三種のアヤメの花について、sepal length(がくの長さ), sepal width(がくの幅), petal length (花弁の長さ)そして petal width(花弁の幅) の測定値が特徴量として格納されている150個 のデータセットです。機械学習用の標準データセットとして、下記のサンプルコードで簡単にインストールできます。

◆決定木とは

決定木は、分類器の一つ。n次元の特徴量からなるデータの多クラス分類が可能です。分類のための判断プロセスが可視化できるため、意思決定などに多く用いられます。

例を見てみましょう。「服を買うか否か」の分類問題で、特徴量は「ほしい服があるか(2値)」「店までの距離は何分か(連続値)」「天候が晴れか(2値)」であったとします。決定木は、分類が可能となる”質問”の仕方(フロー)であり、決定木学習ではこの質問の最適な仕方・順序を学習します。
最適な質問の仕方の決定には情報利得の最大化という考え方を用います。具体的には、ある特徴量について、境界線(ある値以上or以下、など)を引いてみたときに、分類が最もうまくいくものを探す、という手法です。分類がうまくいっているかどうかについては、分類されたものの純度の高さ、情報量の高さで評価します。
評価指標としては「ジニ不純度」・「エントロピー」・「分類誤差」などがありますが、ここでは詳細を割愛します。

決定木では、特徴量が、数値特徴量(距離など)でも、順序特徴量(大・中・小など)でも、名義特徴量(アメリカ・日本、、など)でも分類が可能です。

 

 

◆サンプルプログラムの大枠の流れ

Irisのデータセットから二つの特徴量、petal width, petal lengthを持つ3種類の花のデータセット

1.irisデータセットから二つの特徴量、petal width, petal lengthを持つ3種類の花のデータセットをインポートし、学習データと検証データに切り分けます。

2.学習データを用いて、scikit-learnにある決定木を適用し、分類器(質問を学習させます)を学習させます。ここでは純度の評価指標として「エントロピー」を用います。決定木の「深さ」も設定することができ、ここではmax_depth=3としています。深すぎるときに過学習に陥りやすくなるため注意が必要です。

3.検証用データを学習させた分類器で分類し、正答率を検証します。

4.「Graphviz」を用いて、学習させた決定木を可視化します。
※Graphvizの公式ダウンロードページhttp://www.graphviz.org/Download.php

 

 

#=========================================================
## data import
from distutils.version import LooseVersion as Version
from sklearn import __version__ as sklearn_version
from sklearn import datasets
import numpy as np

iris = datasets.load_iris()
X = iris.data[:, [2, 3]]
y = iris.target

#data plot
import matplotlib.pyplot as plt
fig = plt.figure()
plt.scatter(np.ravel(X[y==0,0]),np.ravel(X[y==0,1]),c='b', marker='o',label='1')
plt.scatter(np.ravel(X[y==1,0]),np.ravel(X[y==1,1]),c='r', marker='o',label='2')
plt.scatter(np.ravel(X[y==2,0]),np.ravel(X[y==2,1]),c='y', marker='o',label='3')

plt.xlabel('petal length [cm]')
plt.ylabel('petal width [cm]')
plt.legend(loc='upper left')
plt.tight_layout()
#=========================================================

#=========================================================
# split data into training data and test data
if Version(sklearn_version) < '0.18':
    from sklearn.cross_validation import train_test_split
else:
    from sklearn.model_selection import train_test_split

##Split the data into trainning data and test data at a certain ratio
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.3, random_state=0)
#=========================================================

#standardization
#from sklearn.preprocessing import StandardScaler

#sc = StandardScaler()
#sc.fit(X_train)
#X_train_std = sc.transform(X_train)
#X_test_std = sc.transform(X_test)


from sklearn.tree import DecisionTreeClassifier

tree = DecisionTreeClassifier(criterion='entropy', max_depth=3, random_state=0)
tree.fit(X_train, y_train)

y_pred_dt=tree.predict(X_test)

plt.figure(figsize=(8, 6), dpi=80)
# training data
plt.scatter(np.ravel(X_train[y_train==0,0]),np.ravel(X_train[y_train==0,1]),c='b', marker='o',label='train 1')
plt.scatter(np.ravel(X_train[y_train==1,0]),np.ravel(X_train[y_train==1,1]),c='r', marker='o',label='train2')
plt.scatter(np.ravel(X_train[y_train==2,0]),np.ravel(X_train[y_train==2,1]),c='y', marker='o',label='train3')
#predction data
plt.scatter(np.ravel(X_test[y_pred_dt==0,0]),np.ravel(X_test[y_pred_dt==0,1]),c='b', marker='o',label='prediction 1',s=60)
plt.scatter(np.ravel(X_test[y_pred_dt==1,0]),np.ravel(X_test[y_pred_dt==1,1]),c='r', marker='o',label='prediction 2',s=60)
plt.scatter(np.ravel(X_test[y_pred_dt==2,0]),np.ravel(X_test[y_pred_dt==2,1]),c='y', marker='o',label='prediction 3',s=60)
#test data
plt.scatter(np.ravel(X_test[y_test==0,0]),np.ravel(X_test[y_test==0,1]),c='b', marker='o',label='test 1')
plt.scatter(np.ravel(X_test[y_test==1,0]),np.ravel(X_test[y_test==1,1]),c='r', marker='o',label='test 2')
plt.scatter(np.ravel(X_test[y_test==2,0]),np.ravel(X_test[y_test==2,1]),c='y', marker='o',label='test 3')

plt.xlabel('petal length [cm]')
plt.ylabel('petal width [cm]')
plt.legend(loc=2)
plt.tight_layout()

plt.show()

# prediction accuracy
from sklearn.metrics import accuracy_score
print('Misclassified samples: %d' % (y_test != y_pred_dt).sum())
print('Accuracy : %.2f' % accuracy_score(np.rint(y_pred_dt), y_test))


## export dot file
from sklearn.tree import export_graphviz

export_graphviz(tree, 
                out_file='tree_result.dot', 
                feature_names=['petal length', 'petal width'])


## export dot file
#from sklearn.tree import export_graphviz


◆結果

Misclassified samples: 1
Accuracy : 0.98

すべてのデータとその花の種類、検証データに対する予測した花の種類(大きな〇印)をプロットしています。黄色印と赤印の境界に誤分類が発生しているのがわかります。
次に、決定木の構造を、出力データから読み取ってみましょう。
tree_result.dotというファイルが、作業ディレクトリに保存されています。
GVEdit.exeファイルを実行して、 File → Openからtree_result.dotを開いてください。
すると以下のような決定木の構造が表示されます。

分類のための質問に番号を赤字で振っています。(Question1,Question2,,,)
各質問の内容が書かれており、境界の値として何を使っているかが明示されています。

それぞれの質問の境界線を載せてみます。学習として、まずQuestion1:petal width で青を一気に分類できる境界線を引きます。次にQestion2: petal lengthで黄と赤をざっくり分けて行きます。決定木の末端であるQuestion 3,4で赤と黄をさらに分類していきます。

このように、この決定木学習の境界面は、多くの直線(平面)で分断した構造を持っています。

Python 機械学習 Scikit-learnによるサポートベクターマシン(SVM)による学習・分類の実践

今回はSVM(サポートベクターマシン)によるIrisの花のデータセットの分類を実践します。サンプルプログラムをコピーペーストでまずは動かしてみることを念頭にしています。

※Pythonの導入方法http://costudyinfodatabase.nagoya/2017/01/05/%e3%81%82-2/

※Iris のデータセットとは
Irisのデータセットは「Iris setosa」「 Iris versicolor」「 Iris virginica 」の三種のアヤメの花について、sepal length(がくの長さ), sepal width(がくの幅), petal length (花弁の長さ)そして petal width(花弁の幅) の測定値が特徴量として格納されている150個 のデータセットです。機械学習用の標準データセットとして、下記のサンプルコードで簡単にインストールできます。

◆SVMとは

SVMは、分類器の一つ。n次元の特徴量からなるデータ空間に、トレーニングデータで学習して、分類境界線を引きます。そして、その境界線を境にして、未知のデータがどのクラスに属するかを判定します。

線形SVM
特徴量が2次元、クラスが二つ(下では、青とオレンジで色分け)の場合のイメージを以下に示します。2つのクラスの間に境界面を引き、2つのクラスの内、この境界面に最も近いサンプルと境界との距離(マージン)が最も大きくなるようにな境界面を決定境界とするのが線形SVMです。

 

非線形SVM
上記のように、境界面が直線(平面)であるようなSVMを線形SVMと呼んでいます。それに対し、決定境界に曲面を用いる非線形SVMというものがあります。下のようなクラスの分布を持っている場合はどうでしょうか。青とオレンジの2つのクラスを直線で分離するのが困難であり、下図のような曲面の境界が有効であることがわかります。
(曲面の関数を引いてやる手法はここでは割愛。非線形SVMで検索すると大量にヒットします)

いずれの手法を用いるにしても、設定した特徴量空間でクラスが分類できるような集合を成しているか(集合がごちゃ混ぜで境界もなにも引けないような状態になっていないか)が高い識別率を左右します。

◆サンプルプログラムの大枠の流れ

1.irisのデータセットから3次元データ(3つの特徴量を有するデータ)をインポートし、学習データと検証データに切り分けます。

2.学習データを用いて、scikit-learnにある非線形SVMを適用し、分類器(どんな境界線を引くか)を学習させます。

3.検証用データを学習させた分類器で分類し、正答率を検証します。

 

 

irisのデータセットから3次元データ(3つの特徴量を有するデータ)をインポートし、学習データと検証データに切り分けます。切り分けたデータには標準化を行います。

※標準化:データの平均値=0、分散=1になるようにデータを整形すること。

#=========================================================
## data import
from distutils.version import LooseVersion as Version
from sklearn import __version__ as sklearn_version
from sklearn import datasets
import numpy as np

iris = datasets.load_iris()
data1=1
data2=2
data3=3

X = iris.data[:, [data1,data2,data3]]#sepal width/petal length/petal length
y = iris.target

print('Class labels:', np.unique(y))

if Version(sklearn_version) &lt; '0.18':
from sklearn.cross_validation import train_test_split
else:
from sklearn.model_selection import train_test_split

##Split the data into trainning data and test data at a certain ratio
X_train, X_test, y_train, y_test = train_test_split(
X, y, test_size=0.3, random_state=0)

from sklearn.preprocessing import StandardScaler

sc = StandardScaler()
sc.fit(X_train)
X_train_std = sc.transform(X_train)
X_test_std = sc.transform(X_test)
#==========================================================

#==========================================================
#data plot
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

fig = plt.figure()
ax = Axes3D(fig)
#ax.plot_wireframe(X,Y,Z)
#ax.scatter3D(np.ravel(X1),np.ravel(X2),np.ravel(X3),c='b', marker='x',label='1')
ax.scatter3D(np.ravel(X[y==0,0]),np.ravel(X[y==0,1]),np.ravel(X[y==0,2]),c='b', marker='x',label='1')
ax.scatter3D(np.ravel(X[y==1,0]),np.ravel(X[y==1,1]),np.ravel(X[y==1,2]),c='r', marker='o',label='2')
ax.scatter3D(np.ravel(X[y==2,0]),np.ravel(X[y==2,1]),np.ravel(X[y==2,2]),c='y', marker='o',label='3')
ax.set_xlabel(iris.feature_names[data1])
ax.set_ylabel(iris.feature_names[data2])
ax.set_zlabel(iris.feature_names[data3])
plt.legend(loc='upper left')

plt.show()
#========================================================

実行結果

3次元の特徴量空間に3種類の花が色分けしてプロットされています。

 

次に、scikit-learnライブラリからSVMをインポートします。

SVMとして、非線形SVMを選択します。kernel=’rbf’ : 動径基底関数カーネルを用いた非線形SVM。’linear’を選択すると線形SVMを利用できます。ハイパーパラメータであるγを小さくするとサンプルデータの寄与が大きくなり、決定境界がより学習サンプルに追従するように学習されます。詳しくは割愛。

学習させた分類器にテストデータを分類させ、クラス(花の種類)を予測させます。実際の花の種類と、予測された種類を比較し、正答率を表示します。

次に、全データと、テストデータをその予測された花の種類の色で塗りつぶしたものを大きめの〇で囲ってプロットするようにしています。

from sklearn.svm import SVC

## apply non-linear svm to iris
svm = SVC(kernel='rbf', random_state=0, gamma=0.2, C=1.0)
svm.fit(X_train_std, y_train)
y_pred_svm=svm.predict(X_test_std)

# prediction accuracy
from sklearn.metrics import accuracy_score
print('Misclassified samples: %d' % (y_test != y_pred_svm).sum())
print('Accuracy : %.2f' % accuracy_score(np.rint(y_pred_svm), y_test))

#test data plot
fig1 = plt.figure()
ax1 = Axes3D(fig1)
ax1.scatter3D(np.ravel(X_test[y_pred_svm==0,0]),np.ravel(X_test[y_pred_svm==0,1]),np.ravel(X_test[y_pred_svm==0,2]),c='b', marker='x',label='1')
ax1.scatter3D(np.ravel(X_test[y_pred_svm==1,0]),np.ravel(X_test[y_pred_svm==1,1]),np.ravel(X_test[y_pred_svm==1,2]),c='r', marker='o',label='2')
ax1.scatter3D(np.ravel(X_test[y_pred_svm==2,0]),np.ravel(X_test[y_pred_svm==2,1]),np.ravel(X_test[y_pred_svm==2,2]),c='y', marker='o',label='3')
ax1.set_xlabel(iris.feature_names[data1])
ax1.set_ylabel(iris.feature_names[data2])
ax1.set_zlabel(iris.feature_names[data3])
plt.legend(loc='upper left')

#============================================================
#all data plot
fig = plt.figure()
ax2 = Axes3D(fig)
#train data
ax2.scatter3D(np.ravel(X[y==0,0]),np.ravel(X[y==0,1]),np.ravel(X[y==0,2]),c='b', marker='x',label='1')
ax2.scatter3D(np.ravel(X[y==1,0]),np.ravel(X[y==1,1]),np.ravel(X[y==1,2]),c='r', marker='o',label='2')
ax2.scatter3D(np.ravel(X[y==2,0]),np.ravel(X[y==2,1]),np.ravel(X[y==2,2]),c='y', marker='o',label='3')
#test data
ax2.scatter3D(np.ravel(X_test[y_pred_svm==0,0]),np.ravel(X_test[y_pred_svm==0,1]),np.ravel(X_test[y_pred_svm==0,2]),c='b', marker='x',label='1')
ax2.scatter3D(np.ravel(X_test[y_pred_svm==1,0]),np.ravel(X_test[y_pred_svm==1,1]),np.ravel(X_test[y_pred_svm==1,2]),c='r', marker='o',label='2')
ax2.scatter3D(np.ravel(X_test[y_pred_svm==2,0]),np.ravel(X_test[y_pred_svm==2,1]),np.ravel(X_test[y_pred_svm==2,2]),c='y', marker='o',label='3')
#highlight the test data
ax2.scatter3D(np.ravel(X_test[y_pred_svm==0,0]),np.ravel(X_test[y_pred_svm==0,1]),np.ravel(X_test[y_pred_svm==0,2]),c='', marker='o',label='1', s=60)
ax2.scatter3D(np.ravel(X_test[y_pred_svm==1,0]),np.ravel(X_test[y_pred_svm==1,1]),np.ravel(X_test[y_pred_svm==1,2]),c='', marker='o',label='2', s=60)
ax2.scatter3D(np.ravel(X_test[y_pred_svm==2,0]),np.ravel(X_test[y_pred_svm==2,1]),np.ravel(X_test[y_pred_svm==2,2]),c='', marker='o',label='3', s=60)

ax2.set_xlabel(iris.feature_names[data1])
ax2.set_ylabel(iris.feature_names[data2])
ax2.set_zlabel(iris.feature_names[data3])
#============================================================

実行結果

 

Misclassified samples: 1
Accuracy : 0.98

以上のように、98%の正答率で分類できています。おおむねうまく分類されていますが、クラスの集合の境目付近は誤分類が発生しています。

 

Python 機械学習 Scikit-learnによるクラスタリング分析(k-means法)の実践

Pythonでk-means法を用いてクラスタリングを実行します。とりあえず動かせるように、サンプルコードを用意しています。

※Pythonの導入方法(http://costudyinfodatabase.nagoya/2017/01/05/%e3%81%82-2/

今回は、irisのデータセットから三次元データ(3つの特徴量を有するデータ)をインポートし、それに対して、scikit-learnにあるk-means法を適用してクラスタリング分析を実行します。可視化のために3次元データを使っていますが、原理上、次元に制約はありません。
また、クラスタ数の最適化の可視化手段であるエルボー法を用いてみます。

 

◆クラスタリング・k-means法とは?

クラスタリングとは教師なし学習法の一つです。データの集合体から似た者同士をグルーピングしていきます。教師あり学習法と異なり、クラスタリングされた結果の解釈は後から行うことになります。

k-means法は、クラスタリングの手法の一つで、計算効率が良く、産業界・学術界で広く用いられています。映画・音楽等のリコメンデーションエンジンなどに応用されています。

クラスタリング手法は大枠以下のように分類できます。

  1. プロトタイプベース クラスタリング
  2. 階層的クラスタリング
  3. 密度ベースクラスタリング

k-means法はこのうち1のプロトタイプベース クラスタリングにあたります。非階層的クラスタリングとも言います。以下に載せている手法では、クラスタ毎にセントロイドという代表点を設けて、サンプル各点のセントロイドに対する近さでグルーピングを行います。具体的には

  • Step1 クラスタ数kを決定し、その代表点であるk個のセントロイドを全サンプルからランダムに抽出。
  • Step2 各サンプルを最も距離が近いセントロイドに割り当てる。(各サンプルが所属するクラスタを仮決め)
  • Step3 各クラスタのサンプルの重心にセントロイドを移す。

セントロイドの位置が動かなくなるまで2-3を繰り返し、セントロイドの位置、およびクラスタを決定するという流れです。

 

◆サンプルコード

実際に動かしてみましょう。

まずはscikit-learn標準のIrisのデータを分類することを念頭に、データをインポートし可視化します。

Irisのデータセットは「Iris setosa」「 Iris versicolor」「 Iris virginica 」の三種のアヤメの花について、sepal length(がくの長さ), sepal width(がくの幅), petal length (花弁の長さ)そして petal width(花弁の幅) の測定値が特徴量として格納されている150個 のデータセットです。

今回は3つの特徴量、sepal width(がくの幅), petal length (花弁の長さ)そして petal width(花弁の幅)のみ抽出します。

## data import
from distutils.version import LooseVersion as Version
from sklearn import __version__ as sklearn_version
from sklearn import datasets
import numpy as np

iris = datasets.load_iris()
names=iris.target_names
data1=1
data2=2
data3=3

X = iris.data[:, [data1,data2,data3]]#sepal width/petal length/petal length
y = iris.target


print('Class labels:', np.unique(y))

if Version(sklearn_version) &amp;amp;amp;amp;lt; '0.18':
    from sklearn.cross_validation import train_test_split
else:
    from sklearn.model_selection import train_test_split

##Split the data into trainning data and test data at a certain ratio
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.3, random_state=0)

from sklearn.preprocessing import StandardScaler

sc = StandardScaler()
sc.fit(X_train)
X_train_std = sc.transform(X_train)
X_test_std = sc.transform(X_test)

#==========================================================
#data plot
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

fig = plt.figure()
ax = Axes3D(fig)
ax1=ax.scatter3D(np.ravel(X[y==0,0]),np.ravel(X[y==0,1]),np.ravel(X[y==0,2]),c='b', marker='o')
ax2=ax.scatter3D(np.ravel(X[y==1,0]),np.ravel(X[y==1,1]),np.ravel(X[y==1,2]),c='r', marker='o')
ax3=ax.scatter3D(np.ravel(X[y==2,0]),np.ravel(X[y==2,1]),np.ravel(X[y==2,2]),c='y', marker='o')
ax.set_xlabel(iris.feature_names[data1])
ax.set_ylabel(iris.feature_names[data2])
ax.set_zlabel(iris.feature_names[data3])
plt.legend((ax1, ax2, ax2),
           (names[0], names[1], names[2]),
           scatterpoints=1,
           loc='upper left',
           ncol=3,
           fontsize=8)
plt.show()

#========================================================

・実行結果

次に、k-means法を用いてクラスタリングし、可視化します。

from sklearn.cluster import KMeans

km = KMeans(n_clusters=3, 
            init='random', 
            n_init=10, 
            max_iter=300,
            tol=1e-04,
            random_state=0)
y_km = km.fit_predict(X)

#Predicted data plot
fig = plt.figure()
ax = Axes3D(fig)
#ax.plot_wireframe(X,Y,Z)
#==============================================================================
# #original data 
#ax1=ax.scatter3D(np.ravel(X[y==0,0]),np.ravel(X[y==0,1]),np.ravel(X[y==0,2]),c='b', marker='o')
#ax2=ax.scatter3D(np.ravel(X[y==1,0]),np.ravel(X[y==1,1]),np.ravel(X[y==1,2]),c='r', marker='o')
#ax3=ax.scatter3D(np.ravel(X[y==2,0]),np.ravel(X[y==2,1]),np.ravel(X[y==2,2]),c='y', marker='o')
 #ax.scatter3D(np.ravel(X1),np.ravel(X2),np.ravel(X3),c='b', marker='x',label='1')
#==============================================================================
ax_km1=ax.scatter3D(np.ravel(X[y_km==0,0]),np.ravel(X[y_km==0,1]),np.ravel(X[y_km==0,2]),c='y', marker='x')
ax_km2=ax.scatter3D(np.ravel(X[y_km==1,0]),np.ravel(X[y_km==1,1]),np.ravel(X[y_km==1,2]),c='r', marker='x')
ax_km3=ax.scatter3D(np.ravel(X[y_km==2,0]),np.ravel(X[y_km==2,1]),np.ravel(X[y_km==2,2]),c='b', marker='x')
ax_km4=ax.scatter3D(km.cluster_centers_[:, 0],km.cluster_centers_[:, 1],km.cluster_centers_[:, 2],c='lightgreen', marker='s')
ax.set_xlabel(iris.feature_names[data1])
ax.set_ylabel(iris.feature_names[data2])
ax.set_zlabel(iris.feature_names[data3])
plt.legend((ax_km1, ax_km2, ax_km3,ax_km4),
           (names[0], names[1], names[2],"Centroid"),
           scatterpoints=1,
           loc='upper left',
           ncol=3,
           fontsize=8)
plt.show()

 

・実行結果

以下のように、ほぼ元データと同様の分類ができています。

また、黄緑色でセントロイドもプロットしています。

(誤分類率を算出してみましょう)

 

以上の方法ではクラスタ数は自分で設定することになりますが、最適なクラスタ数はどう評価するでしょうか?

簡単な手法であるエルボー法で評価してみましょう。振り分けるクラスタ数を横軸、縦軸にクラスタ内誤差平方和(SSE)をプロットし、SSEが急降下/飽和する点が最適なクラスタ数というわけです。

クラスタ内誤差平方和とは、クラスタ内のサンプルとセントロイドとの距離の総和であり、「クラスタ全体のまとまり度合」のようなものです。

 

##Elbow method###
#Elbow method is to find the saturation point ("Elbow") with changing the cluster numbers
print("Elbow method")
print('Distortion: %.2f' % km.inertia_)
distortions = []
for i in range(1, 11):
    km = KMeans(n_clusters=i, 
                init='k-means++', 
                n_init=10, 
                max_iter=300, 
                random_state=0)
    km.fit(X)
    distortions.append(km.inertia_)
plt.plot(range(1, 11), distortions, marker='o')
plt.xlabel('Number of clusters')
plt.ylabel('Distortion')
plt.tight_layout()
#plt.savefig('./figures/elbow.png', dpi=300)
plt.show()


・実行結果

3個のクラスタで十分であることがわかります。

パーセプトロン・ロジスティック回帰の実践 Python 機械学習 Scikit-learn

Python 標準機械学習ツールであるScikit-learnを用いて、動脈硬化症の発症予測をします。今回もサンプルコードをコピペでとりあえず動かせるようにしています。

※Pythonの導入方法はhttp://costudyinfodatabase.nagoya/2017/01/05/%e3%81%82-2/

元となるデータ(※)は、「動脈硬化症の発生有無」「脂質異常スコア」「性別」「年齢」から構成されます。動脈硬化症の発生有無が目的変数(結果)で、それ以外が説明変数(入力値・特徴量)です。

本データを用いて、パーセプトロンとロジスティック回帰を用いて学習を行い、テストデータの入力に対して動脈硬化の有無を予測し、正答率を見てみましょう。

※出典:http://www.snap-tck.com/room04/c01/stat/stat99/stat0206.pdf

No 動脈硬化症(0=あり) 脂質異常スコア 性(0=男) 年齢
1 0 0 0 36
2 0 0 0 55
3 0 0 1 27
4 0 0 1 42
5 0 1 0 35
6 0 1 0 39
7 0 1 0 41
8 0 1 0 45
9 0 1 1 32
10 0 1 1 42
11 0 1 1 51
12 0 1 1 53
13 0 2 0 43
14 0 2 0 47
15 0 2 1 52
16 1 1 0 46
17 1 1 1 24
18 1 1 1 38
19 1 1 1 58
20 1 2 0 21
21 1 2 0 30
22 1 2 0 37
23 1 2 1 24
24 1 2 1 56
25 1 2 1 58

 

サンプルコード

CSVデータからデータを抽出し、パーセプトロンとロジスティック回帰でトレーニングを行うプログラムです。

ロジスティック回帰では、ある特徴量(X)に対して結果(Y)が2値(0または1)に分類されるときに有効な手法です。

「1」となる事象(正事象)の確率をpとし、あるデータの確率pが0.5を下回るなら「0」、上回るなら「1」に分類します。確率pが適切な予測となるように重みwnを探ります。

今回は以下の流れで実行していきます。

手順

  • Step1 CSVデータをインポートし、トレーニング用データセットと、トレーニング結果検証用データセットに切り分け分けます。※CSVデータは自分で作成してみてください。
  • Step2 トレーニング用にデータの標準化を行います。
  • Step3 Perceptron / Logistic Regression(ロジスティック回帰)各々の手法について、誤分類率を表示します。
  • Step4 Logistic回帰による各データの確率をライブラリを用いて示します。それが上記式から計算した確率に等しいことを確認します。

 

 

import numpy as np

&nbsp;

##csv resd using numpy

csvdata1=np.loadtxt("Arteriosclerosis.csv",delimiter=",")

&nbsp;

X1=csvdata1[:,2] #Lipid abnormality score

X2=csvdata1[:,3] #sex

X3=csvdata1[:,4] #age

Y=csvdata1[:,1] #Arteriosclerosis

X=np.c_[X1,X2,X3] #to match the shape of fitting function's iput

&nbsp;

#==============================================================================

from sklearn.linear_model import LogisticRegression

from distutils.version import LooseVersion as Version

from sklearn import __version__ as sklearn_version

from sklearn.metrics import accuracy_score

&nbsp;

&nbsp;

if Version(sklearn_version) &lt; '0.18':

from sklearn.cross_validation import train_test_split

else:

from sklearn.model_selection import train_test_split

&nbsp;

##Split the data into trainning data and test data at a certain ratio(with shuffuling the data)

X_train, X_test, Y_train, Y_test = train_test_split(

X, Y, test_size=0.2, random_state=0)

&nbsp;

print('Expected value:',Y_test)

print('------------------------------------------------------')

&nbsp;

from sklearn.preprocessing import StandardScaler

&nbsp;

sc = StandardScaler()

sc.fit(X_train)

X_train_std = sc.transform(X_train)

X_test_std = sc.transform(X_test)

&nbsp;

# perceptron

from sklearn.linear_model import Perceptron

&nbsp;

ppn = Perceptron(n_iter=40, eta0=0.1, random_state=0)

ppn.fit(X_train_std, Y_train)

&nbsp;

#print("shape of y_test",Y_test.shape)

&nbsp;

y_pred_perceptron = ppn.predict(X_test_std)

print('Classified by perceptron:',y_pred_perceptron)

print('Misclassified samples: %d' % (Y_test != y_pred_perceptron).sum())

print('Accuracy of Perceptron: %.2f' % accuracy_score(Y_test, y_pred_perceptron))

print('Bias:',ppn.intercept_)

print('Weights after fitting:',ppn.coef_)

print('------------------------------------------------------')

&nbsp;

&nbsp;

#Logistic Regression

lr = LogisticRegression(C=1000.0, random_state=0)

lr.fit(X_train_std, Y_train)

y_pred_lr=lr.predict(X_test_std)

print('Classified by Logistic Regression:',y_pred_lr)

print('Misclassified samples: %d' % (Y_test != y_pred_lr).sum())

print('Bias:',lr.intercept_)

print('Weights after fitting:',lr.coef_)

print('Accuracy of Logistic Regression: %.2f' % accuracy_score(Y_test, y_pred_lr))

&nbsp;

##additional This is to see the probability calculated based on the relation

## equals to the probability calculated in the library

print('probability p,1-p of all test data:\n',lr.predict_proba(X_test_std))

print('probability p of the 1st test datum',(1 / (1 + np.exp(np.dot(X_test_std,-lr.coef_.T) -lr.intercept_)))[0])

&nbsp;

実行結果

以下のような結果となるはずです。いかがでしょうか?この場合正答率は80%となります。

 

Expected value: [ 0.  0.  1.  1.  0.]——————————————————

Classified by perceptron: [ 0.  0.  1.  0.  0.]

Misclassified samples: 1

Accuracy of Perceptron: 0.80

Bias: [-0.1]

Weights after fitting: [[ 0.17902872 -0.1         0.01490374]]

——————————————————

Classified by Logistic Regression: [ 0.  0.  1.  0.  0.]

Misclassified samples: 1

Bias: [-0.58884717]

Weights after fitting: [[ 1.15017829  0.5128515  -0.11998034]]

Accuracy of Logistic Regression: 0.80

probability p,1-p of all test data:

[[ 0.81112387  0.18887613]

[ 0.87276998  0.12723002]

[ 0.39654165  0.60345835]

[ 0.55972308  0.44027692]

[ 0.64804161  0.35195839]]

probability p of the 1st test datum [ 0.18887613]