geek開発日誌

超一流のプログラマーを目指してます!作成したポートフォリオを記録していく開発日誌です。

教師なし機械学習(クラスタリング)

こんにちは。nakatatsuです。

今回も機械学習について学んでいきます。基本的にはPythonとscikit-learnライブラリを使用します。
参考書は
Pythonではじめる機械学習オライリー・ジャパン
を使用します。

今回は教師なし学習の「クラスタリングアルゴリズムです。クラスタリングはデータセットクラスタと呼ばれるグループに分割するアルゴリズムです。同じクラスタ内のデータは類似していて、異なるクラスタのデータは異なるようにデータを分割します。


種類

今回は以下の3つの方法を用います。

・k-means
k-meansは最も単純で最も広く用いられているクラスタリングアルゴリズムです。このアルゴリズムはデータのある領域を代表するようなクラスタセンタを見つけようとします。
3つのクラスタセンタを探す過程を下図を用いて説明します。まず、ランダムで3つのデータポイントをクラスタセンタとして選びます。次に、個々のデータポイントが最も近いクラスタセンタに割り当てられます。そして、割り当てられた点の平均にクラスタセンタを更新します。この過程を繰り返し、クラスタセンタが動かなくなったら終了します。
新しいデータポイントが与えられると、クラスタセンタのうち、最も近いものに割り当てられます。

f:id:nakatatsu_com:20190916222601p:plain

・凝集型クラスタリング
凝集型クラスタリングは個々のデータポイントをそれぞれ個別のクラスタとして開始し、最も類似した2つのクラスタを併合していきます。
3つのクラスタを生成する過程を下図を用いて説明します。最初は個々のデータポイントがクラスタになっており、ステップごとに最も近い2つのクラスタが併合されていきます。指定されたクラスタが生成されるまでこれを繰り返し、クラスタの併合を続けていきます。
アルゴリズムの動作の関係上、凝集型クラスタリングは新しいデータに対して予測をすることができません。

f:id:nakatatsu_com:20190916224753p:plain

・DBSCAN
DBSCANは、クラスタはデータ中で高密度領域を構成していて、比較的空虚な領域で区切られているという考えに基づいています。
DBSCANにはmin_samples、epsという2つのパラメータがあります。あるデータポイントから距離eps以内にmin_samples以上のデータポイントがある場合、そのデータポイントをコアサンプルと呼びます。
まず、適当に1つのデータポイントを選び、そのデータポイントから距離eps以内にあるすべてのデータポイントを見つけ、その数がmin_samples以上であればその点をコアサンプル、min_samples以下であればその点をノイズと呼びます。次に、eps以内にあるすべての点をテストし、それらの点がまだクラスタに割り当てられていなければクラスタを割り当てます。これを境界ポイントと呼びます。その点がコアサンプルであればその近傍をさらにテストします。これをすべての点で繰り返し、最終的にコアポイント、境界ポイント、ノイズの3種類のデータポイントに分類されます。
DBSCANは凝集型クラスタリングと同じように新しいデータに対して予測をすることができません。


準備

Anacondaを使用するための準備についてはこのブログの過去記事「教師あり機械学習(k-最近傍法)」を参照してください。

nakatatsu-com.hatenablog.com


結果

AnacondaのJupyterLabを使用して作成しました。


まずはk-meansです。


流れとしては、

○k-meansの適用
○k-meansが機能しない例
○ベクトル量子化による成分分解

となっています。


import mglearn
import matplotlib.pyplot as plt
import numpy as np


○k-meansの適用
合成データセットに対してk-meansを適用してみます。

# k-meansを適用
from sklearn.datasets import make_blobs
from sklearn.cluster import KMeans

X, y = make_blobs(random_state=1)

# クラスタセンタ数3
kmeans = KMeans(n_clusters=3)
kmeans.fit(X)

mglearn.discrete_scatter(X[:, 0], X[:, 1], kmeans.labels_, markers='o')
mglearn.discrete_scatter(
    kmeans.cluster_centers_[:, 0], kmeans.cluster_centers_[:, 1], [0, 1, 2],
    markers='^', markeredgewidth=2);

f:id:nakatatsu_com:20190916231736p:plain

fig, axes = plt.subplots(1, 2, figsize=(10, 5))

# クラスタセンタ数2
kmeans = KMeans(n_clusters=2)
kmeans.fit(X)
assignments = kmeans.labels_
mglearn.discrete_scatter(X[:, 0], X[:, 1], assignments, ax=axes[0])

# クラスタセンタ数5
kmeans = KMeans(n_clusters=5)
kmeans.fit(X)
assignments = kmeans.labels_
mglearn.discrete_scatter(X[:, 0], X[:, 1], assignments, ax=axes[1]);

f:id:nakatatsu_com:20190916231750p:plain

クラスタリングにはクラスタ数を指定する必要があります。しかし、実世界ではクラスタ数はわからないことも多いです。
クラスタリングとクラス分類は両方ともラベル付けをするという意味では同じです。しかし、クラスタリングには真のラベルというものがないので、クラスタリングからわかることは同じラベルのものは相互に似ているということだけです。


○k-meansが機能しない例
k-meansは比較的単純なデータセットにしか機能しません。ここでは、k-meansが機能しない例を見てみます。

# k-meansがうまく機能しない場合
X_varied, y_varied = make_blobs(n_samples=200,
                                cluster_std=[1.0, 2.5, 0.5],
                                random_state=170)
y_pred = KMeans(n_clusters=3, random_state=0).fit_predict(X_varied)

mglearn.discrete_scatter(X_varied[:, 0], X_varied[:, 1], y_pred)
plt.legend(['cluster 0', 'cluster 1', 'cluster 2'], loc='best')
plt.xlabel('Feature 0')
plt.ylabel('Feature 1');

f:id:nakatatsu_com:20190916233128p:plain

X, y = make_blobs(random_state=170, n_samples=600)
rng = np.random.RandomState(74)

transformation = rng.normal(size=(2, 2))
X = np.dot(X, transformation)

kmeans = KMeans(n_clusters=3)
kmeans.fit(X)
y_pred = kmeans.predict(X)

plt.scatter(X[:, 0], X[:, 1], c=y_pred, cmap=mglearn.cm3)
plt.scatter(kmeans.cluster_centers_[:, 0], kmeans.cluster_centers_[:, 1],
            marker='^', c=[0, 1, 2], s=100, linewidth=2, cmap=mglearn.cm3)
plt.xlabel('Feature 0')
plt.ylabel('Feature 1');

f:id:nakatatsu_com:20190916233141p:plain

from sklearn.datasets import make_moons

X, y = make_moons(n_samples=200, noise=0.05, random_state=0)

kmeans = KMeans(n_clusters=2)
kmeans.fit(X)
y_pred = kmeans.predict(X)

plt.scatter(X[:, 0], X[:, 1], c=y_pred, cmap=mglearn.cm2, s=60)
plt.scatter(kmeans.cluster_centers_[:, 0], kmeans.cluster_centers_[:, 1],
            marker='^', c=[mglearn.cm2(0), mglearn.cm2(1)], s=100, linewidth=2)
plt.xlabel('Feature 0')
plt.ylabel('Feature 1');

f:id:nakatatsu_com:20190916233153p:plain

k-meansは、異なる濃度を持ったクラスタ、丸くないクラスタ、two_moonsデータセットのような複雑な形状に対してはうまく機能しないことがわかります。


○ベクトル量子化による成分分解
PCAやNMFはデータポイントを複数の成分の和として表現しますが、これに対してk-meansはクラスタセンタで個々のデータポイントを表現します。これをベクトル量子化と呼びます。
「Labeled Faces in the Wild」データセットに対して、PCA、NMF、k-meansを適用し、抽出された100成分で比較してみます。さらに、100成分で再構築した画像についても比較してみます。k-meansについてはクラスタセンタのうち最も近いものを用います。

# ベクトル量子化
from sklearn.datasets import fetch_lfw_people
from sklearn.model_selection import train_test_split
from sklearn.decomposition import PCA
from sklearn.decomposition import NMF

people = fetch_lfw_people(min_faces_per_person=20, resize=0.7)
image_shape = people.images[0].shape
mask = np.zeros(people.target.shape, dtype=np.bool)
for target in np.unique(people.target):
    mask[np.where(people.target == target)[0][:50]] = 1
    
X_people = people.data[mask]
y_people = people.target[mask]

X_train, X_test, y_train, y_test = train_test_split(
    X_people, y_people, stratify=y_people, random_state=0)
nmf = NMF(n_components=100, random_state=0)
nmf.fit(X_train)
pca = PCA(n_components=100, random_state=0)
pca.fit(X_train)
kmeans = KMeans(n_clusters=100, random_state=0)
kmeans.fit(X_train)

X_reconstructed_pca = pca.inverse_transform(pca.transform(X_test))
X_reconstructed_kmeans = kmeans.cluster_centers_[kmeans.predict(X_test)]
X_reconstructed_nmf = np.dot(nmf.transform(X_test), nmf.components_)
# 抽出された成分
fig, axes = plt.subplots(3, 5, figsize=(8, 8),
                         subplot_kw={'xticks': (), 'yticks': ()})
fig.suptitle('Extracted Components')
for ax, comp_kmeans, comp_pca, comp_nmf in zip(
    axes.T, kmeans.cluster_centers_, pca.components_, nmf.components_):
    ax[0].imshow(comp_kmeans.reshape(image_shape))
    ax[1].imshow(comp_pca.reshape(image_shape), cmap='viridis')
    ax[2].imshow(comp_nmf.reshape(image_shape))
    
axes[0, 0].set_ylabel('kmeans')
axes[1, 0].set_ylabel('pca')
axes[2, 0].set_ylabel('nmf')

# 再構築
fig, axes = plt.subplots(4, 5, figsize=(8, 8),
                         subplot_kw={'xticks': (), 'yticks': ()})
fig.suptitle('Reconstructions')
for ax, orig, rec_kmeans, rec_pca, rec_nmf in zip(
    axes.T, X_test, X_reconstructed_kmeans, X_reconstructed_pca, X_reconstructed_nmf):
    ax[0].imshow(orig.reshape(image_shape))
    ax[1].imshow(rec_kmeans.reshape(image_shape))
    ax[2].imshow(rec_pca.reshape(image_shape))
    ax[3].imshow(rec_nmf.reshape(image_shape))
    
axes[0, 0].set_ylabel('original')
axes[1, 0].set_ylabel('kmeans')
axes[2, 0].set_ylabel('pca')
axes[3, 0].set_ylabel('nmf');

f:id:nakatatsu_com:20190916234805p:plain

f:id:nakatatsu_com:20190916234817p:plain

k-meansを用いたベクトル量子化では、入力次元よりもはるかに多くのクラスタを使うことができます。例えば、2次元しかないデータセットに対して、PCAやNMFにはできることはほとんどありませんが、多数のクラスタセンタでk-meansを用いれば、より強力な表現を見つけることができます。



次は凝集型クラスタリングです。


流れとしては、

○凝集型クラスタリングの適用
○階層型クラスタリング、デンドログラム

となっています。


○凝集型クラスタリングの適用
合成データセットに対して凝集型クラスタリングを適用してみます。

# 凝縮型クラスタリングを適用
from sklearn.cluster import AgglomerativeClustering

X, y = make_blobs(random_state=1)

agg = AgglomerativeClustering(n_clusters=3)
assignment = agg.fit_predict(X)

mglearn.discrete_scatter(X[:, 0], X[:, 1], assignment)
plt.xlabel('Feature 0')
plt.ylabel('Feature 1');

f:id:nakatatsu_com:20190917000255p:plain

クラスタは完全に再現されていることがわかります。
しかし、凝集型クラスタリングもk-meansと同じようにクラスタ数を指定しなければなりません。


○階層型クラスタリング、デンドログラム
ここでは、正しいクラスタ数を選択する手助けをしてくれる方法を見てみます。
凝集型クラスタリングは、いわゆる階層型クラスタリングを行い、これを可視化することで、階層化の詳細を理解することができます。
また、別の方法としてデンドログラムと呼ばれる方法があります。

# 階層型クラスタリング
mglearn.plots.plot_agglomerative()

f:id:nakatatsu_com:20190917001349p:plain

# デンドログラム
from scipy.cluster.hierarchy import dendrogram, ward

X, y = make_blobs(random_state=0, n_samples=12)
linkage_array = ward(X)
dendrogram(linkage_array)

ax = plt.gca()
bounds = ax.get_xbound()
ax.plot(bounds, [7.25, 7.25], '--', c='k')
ax.plot(bounds, [4, 4], '--', c='k')

ax.text(bounds[1], 7.25, ' two clusters', va='center', fontdict={'size': 15})
ax.text(bounds[1], 4, ' three clusters', va='center', fontdict={'size': 15})
plt.xlabel('Sample index')
plt.ylabel('Cluster distance');

f:id:nakatatsu_com:20190917001332p:plain

階層型クラスタリングでは、それぞれのクラスタがより小さいクラスタに分割されていることがわかります。この可視化によって階層化の詳細を理解することができます。しかし、この方法はデータが2次元であるということに依存しているので、多次元のデータセットを扱うことができません。
これに対して、デンドログラムは多次元のデータセットを扱うことができます。y軸はいつ2つのクラスタを併合したかを示しているだけではなく、枝の長さが2つのクラスタがどれだけ離れていたかを示しています。このことから、このデンドログラムでは3クラスタから2クラスタに併合する際に、かなり離れた点を併合していることがわかります。



最後はDBSCANです。


流れとしては、

○異なるeps、min_samplesによるクラスタリング
○two_moonsデータセットに対してDBSCANを適用

となっています。


○異なるeps、min_samplesによるクラスタリング
DBSCANの主な利点はクラスタ数を指定する必要がない点です。ここでは、合成データセットに対してDBSCANを適用し、epsとmin_samplesをさまざまな値に設定してクラスタリングを行います。

# DBSCANを適用
from sklearn.cluster import DBSCAN

dbscan = DBSCAN()
clusters = dbscan.fit_predict(X)
print('Cluster memberships:\n{}'.format(clusters))
Cluster memberships:
[-1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1]
mglearn.plots.plot_dbscan()
min_samples: 2 eps: 1.000000  cluster: [-1  0  0 -1  0 -1  1  1  0  1 -1 -1]
min_samples: 2 eps: 1.500000  cluster: [0 1 1 1 1 0 2 2 1 2 2 0]
min_samples: 2 eps: 2.000000  cluster: [0 1 1 1 1 0 0 0 1 0 0 0]
min_samples: 2 eps: 3.000000  cluster: [0 0 0 0 0 0 0 0 0 0 0 0]
min_samples: 3 eps: 1.000000  cluster: [-1  0  0 -1  0 -1  1  1  0  1 -1 -1]
min_samples: 3 eps: 1.500000  cluster: [0 1 1 1 1 0 2 2 1 2 2 0]
min_samples: 3 eps: 2.000000  cluster: [0 1 1 1 1 0 0 0 1 0 0 0]
min_samples: 3 eps: 3.000000  cluster: [0 0 0 0 0 0 0 0 0 0 0 0]
min_samples: 5 eps: 1.000000  cluster: [-1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1]
min_samples: 5 eps: 1.500000  cluster: [-1  0  0  0  0 -1 -1 -1  0 -1 -1 -1]
min_samples: 5 eps: 2.000000  cluster: [-1  0  0  0  0 -1 -1 -1  0 -1 -1 -1]
min_samples: 5 eps: 3.000000  cluster: [0 0 0 0 0 0 0 0 0 0 0 0]

f:id:nakatatsu_com:20190917003308p:plain

epsを増やすと、より多くの点がクラスタに含まれるようになり、クラスタが大きくなりますが、複数のクラスタが併合されることにもなります。
min_samplesを増やすと、コアポイントになるデータポイントが少なくなり、より多くのデータポイントがノイズとなります。


○two_moonsデータセットに対してDBSCANを適用
k-means、凝集型クラスタリングではtwo_moonsデータセットに対してうまく機能しませんでした。ここでは、two_moonsデータセットに対してDBSCANを適用してみます。

# two_moonsデータセットに適用
from sklearn.preprocessing import StandardScaler

X, y = make_moons(n_samples=200, noise=0.05, random_state=0)

scaler = StandardScaler()
scaler.fit(X)
X_scaled = scaler.transform(X)

dbscan = DBSCAN()
clusters = dbscan.fit_predict(X_scaled)
plt.scatter(X_scaled[:, 0], X_scaled[:, 1], c=clusters, cmap=mglearn.cm2, s=60)
plt.xlabel('Feature 0')
plt.xlabel('Feature 1');

f:id:nakatatsu_com:20190917004503p:plain

DBSCANではtwo_moonsデータセットに対してうまく機能していることがわかります。うまく機能するためには、適したepsを設定する必要がありますが、StandardScalerやMinMaxScalerでスケール変換してからの方がうまくいく場合が多いです。


まとめ

k-mean、凝集型クラスタリングクラスタ数を指定しなければいけませんが、DBSCANはクラスタ数を指定する必要はなく、自動的にクラスタ数を決めることができます。
k-measは、クラスタセンタによってクラスタの特徴を表すことができます。
凝集型クラスタリングは、階層的クラスタリングやデンドログラムによって階層化の詳細を理解できるので、正しいクラスタ数を選択することができます。
DBSCANは、どのクラスタにも属さないノイズを検出することができます。k-means、凝集型クラスタリングと異なり、複雑な形状のクラスタを見つけることができます。


参考文献

参考文献です。

Pythonではじめる機械学習 ―scikit-learnで学ぶ特徴量エンジニアリングと機械学習の基礎

Pythonではじめる機械学習 ―scikit-learnで学ぶ特徴量エンジニアリングと機械学習の基礎

教師なし機械学習(教師なし変換)

こんにちは。nakatatsuです。

今回も機械学習について学んでいきます。基本的にはPythonとscikit-learnライブラリを使用します。
参考書は
Pythonではじめる機械学習オライリー・ジャパン
を使用します。

今回は教師なし学習の「教師なし変換」アルゴリズムです。元のデータ表現を変換して、人間や他の機械学習アルゴリズムにとって、よりわかりやすい新しいデータ表現を作るアルゴリズムです。

最も一般的な利用法は次元削減です。
次元削減によって、可視化や特徴量の抽出を行います。

・可視化
次元数を2次元に減らし、可視化を行います。

・特徴量の抽出
たくさんの特徴量で構成されるデータの高次元表現を入力として、少量の本質的な特徴を表す特徴量を抽出します。抽出した特徴量を教師あり学習に用いることで精度が上がったり、メモリ使用量や計算時間が削減される場合があります。


種類

今回は可視化、特徴量の抽出を行うために、以下の3つの方法を用います。

・主成分分析(PCA)
PCAは可視化、特徴量の抽出に広く用いられています。
左上のプロットは元のデータを表しています。このアルゴリズムは、まず最も分散が大きい方向を見つけ、それに第1成分というラベルを付けます。データはこの方向に対して最も情報を持ちます。次に第1成分と直交する方向の中から最も情報を持っている方向を探します。2次元だと直交する方向は1つしかありませんが、高次元空間ではたくさんの直交する方向があります。このようにして見つけていく方向を主成分と呼びます。
右上のプロットはそれぞれの主成分がx軸、y軸に沿うように回転させたものです。さらに、原点の周辺にデータが来るように回転させる前にデータから平均値を引いています。
左下のプロットは第1主成分だけを残したものを表しています。このように、主成分のうちいくつかだけを残すことで、PCAを次元削減に使うこともできます。
右下のプロットは逆回転して平均を足し、データを元に戻したものです。第1主成分に含まれている情報しか維持されていないことがわかります。

f:id:nakatatsu_com:20190915163741p:plain

・非負値行列因子分解(NMF)
NMFは特徴量の抽出に用いられています。
PCAの有用な解釈の1つにデータを主成分の重み付き和として表現するというものがあります。NMFでも個々のデータをいくつかの成分の重み付き和として表現します。
PCAでは個々の成分は互いに直交するものである必要がありますが、NMFでは係数と成分が非負である必要があります。
抽出された非負の成分は、原点からデータへの方向だと考えることができます。左のプロットはNMFで見つかった2つの成分を表しています。十分な数の成分がある場合にはデータの極端な部分の方向を向くことがわかります。右のプロットはNMFで見つかった1つの成分を表しています。1つしか成分を使わない場合にはデータの平均値へ向かう成分を作ることがわかります。

f:id:nakatatsu_com:20190915171031p:plain

非負の重み付き和に分解する方法は、いくつもの独立した発生源から得られたデータを重ね合わせて作られるようなデータに対して特に有効です。例えば、複数の人が話している音声データや多数の楽器からなる音楽です。NMFを用いると、組み合わされたデータを作り上げている元の成分を特定することができます。

・t-SNE
t-SNEは可視化に用いられています。
可視化したい場合、PCAを用いる方法が一般的ですが、この手法の性質上、その有用性は限られます。t-SNEアルゴリズムははるかに複雑なマッピングを行い、より良い可視化を実現できます。


準備

Anacondaを使用するための準備についてはこのブログの過去記事「教師あり機械学習(k-最近傍法)」を参照してください。

nakatatsu-com.hatenablog.com


結果

AnacondaのJupyterLabを使用して作成しました。


まずはPCAです。


流れとしては、

○可視化
○特徴量の抽出

となっています。


import matplotlib.pyplot as plt
import mglearn
import numpy as np


○可視化
特徴量30のcancerデータセットをPCAを用いて次元削減し、可視化してみます。

# スケール変換
from sklearn.datasets import load_breast_cancer
from sklearn.preprocessing import StandardScaler

cancer = load_breast_cancer()

scaler = StandardScaler()
scaler.fit(cancer.data)
X_scaled = scaler.transform(cancer.data)
# PCAによる次元削減
from sklearn.decomposition import PCA

pca = PCA(n_components=2)
pca.fit(X_scaled)

X_pca = pca.transform(X_scaled)
print('Original shape: {}'.format(str(X_scaled.shape)))
print('Reduced shape: {}'.format(str(X_pca.shape)))
Original shape: (569, 30)
Reduced shape: (569, 2)
# PCAによる可視化
plt.figure(figsize=(8, 8))
mglearn.discrete_scatter(X_pca[:, 0], X_pca[:, 1], cancer.target)
plt.legend(cancer.target_names, loc='best')
plt.gca().set_aspect('equal')
plt.xlabel('First principal component')
plt.ylabel('Second principal component');

f:id:nakatatsu_com:20190915172947p:plain

print('PCA component shape: {}'.format(pca.components_.shape))
PCA component shape: (2, 30)
# PCAの成分
print('PCA components:\n{}'.format(pca.components_))
PCA components:
[[ 0.21890244  0.10372458  0.22753729  0.22099499  0.14258969  0.23928535
   0.25840048  0.26085376  0.13816696  0.06436335  0.20597878  0.01742803
   0.21132592  0.20286964  0.01453145  0.17039345  0.15358979  0.1834174
   0.04249842  0.10256832  0.22799663  0.10446933  0.23663968  0.22487053
   0.12795256  0.21009588  0.22876753  0.25088597  0.12290456  0.13178394]
 [-0.23385713 -0.05970609 -0.21518136 -0.23107671  0.18611302  0.15189161
   0.06016536 -0.0347675   0.19034877  0.36657547 -0.10555215  0.08997968
  -0.08945723 -0.15229263  0.20443045  0.2327159   0.19720728  0.13032156
   0.183848    0.28009203 -0.21986638 -0.0454673  -0.19987843 -0.21935186
   0.17230435  0.14359317  0.09796411 -0.00825724  0.14188335  0.27533947]]
# PCAの成分(ヒートマップ)
plt.matshow(pca.components_, cmap='viridis')
plt.yticks([0, 1], ['First component', 'Second component'])
plt.colorbar()
plt.xticks(range(len(cancer.feature_names)), cancer.feature_names, rotation=60, ha='left')
plt.xlabel('Feature')
plt.ylabel('Principal components');

f:id:nakatatsu_com:20190915173017p:plain

2つの主成分により、プロットした結果を見てみると、2つのクラスをきれいに分離できていることがわかります。このデータに対してクラス分類を行う際、この結果からどのアルゴリズムを使うかを決めるのに役立ちます。例えば今回の場合、線形モデルでもそれなりに分類できそうだとわかります。
PCAで抽出した主成分の組み合わせを見てみると、主成分には元の特徴量すべてが含まれていることがわかります。


○特徴量の抽出
画像は通常数千ものピクセルで構成されており、それらが集まって初めて意味を持ちます。
しかし、元々の表現よりも解析に適した表現があると予想できるので、「Labeled Faces in the Wild」データセットをPCAに用いて次元削減し、特徴量を抽出してみます。

# 特徴量抽出(「labeled in the Wild」データセット)
from sklearn.datasets import fetch_lfw_people

people = fetch_lfw_people(min_faces_per_person=20, resize=0.7)
image_shape = people.images[0].shape

fix, axes = plt.subplots(2, 5, figsize=(15, 8), subplot_kw={'xticks': (), 'yticks': ()})
for target, image, ax in zip(people.target, people.images, axes.ravel()):
    ax.imshow(image)
    ax.set_title(people.target_names[target])

f:id:nakatatsu_com:20190915173213p:plain

print('people.images.shape: {}'.format(people.images.shape))
print('Number of classes: {}'.format(len(people.target_names)))
people.images.shape: (3023, 87, 65)
Number of classes: 62
mask = np.zeros(people.target.shape, dtype=np.bool)
for target in np.unique(people.target):
    mask[np.where(people.target == target)[0][:50]] = 1
    
X_people = people.data[mask]
y_people = people.target[mask]

X_people = X_people / 255.0
# k-最近傍法を使用
from sklearn.neighbors import KNeighborsClassifier
from sklearn.model_selection import train_test_split

X_train, X_test, y_train, y_test = train_test_split(
    X_people, y_people, stratify=y_people, random_state=0)
knn = KNeighborsClassifier(n_neighbors=1)
knn.fit(X_train, y_train)
print('Test set score of 1-nn: {:.2f}'.format(knn.score(X_test, y_test)))
Test set score of 1-nn: 0.23
# PCAで特徴量を抽出し、k-最近傍法を使用
pca = PCA(n_components=100, whiten=True, random_state=0).fit(X_train)
X_train_pca = pca.transform(X_train)
X_test_pca = pca.transform(X_test)

print('X_train_pca.shape: {}'.format(X_train_pca.shape))
X_train_pca.shape: (1547, 100)
knn = KNeighborsClassifier(n_neighbors=1)
knn.fit(X_train_pca, y_train)
print('Test set score of 1-nn: {:.2f}'.format(knn.score(X_test_pca, y_test)))
Test set score of 1-nn: 0.31
print('pca.components_.shape: {}'.format(pca.components_.shape))
pca.components_.shape: (100, 5655)
# 最初の15主成分
fix, axes = plt.subplots(3, 5, figsize=(15, 12), subplot_kw={'xticks': (), 'yticks': ()})
for i, (component, ax) in enumerate(zip(pca.components_, axes.ravel())):
    ax.imshow(component.reshape(image_shape), cmap='viridis')
    ax.set_title('{}.component'.format(i + 1))

f:id:nakatatsu_com:20190915173257p:plain

mglearn.plots.plot_pca_faces(X_train, X_test, image_shape)

f:id:nakatatsu_com:20190915173315p:plain

元のデータに1ー最近傍法を用いて精度を導出してみると、23%という結果になりました。
次に、PCAにより5655次元(87×65)から100次元に削減してから、1ー最近傍法を用いて精度を導出してみると、31%まで上昇しました。このことから、主成分がデータのより良い表現となっていることがわかります。
抽出された最初の15成分を見てみると、最初の主成分は顔と背景のコンストラストをコーディングしていて、2つ目の主成分は光の当たり方による顔の左右の明るさの差をコーディングしていると予想できます。このことから、人間が解釈する方法とはまったく異なった方法でデータを解釈していることがわかります。
抽出した成分によって元画像を再構築してみると、主成分を増やしていくにつれて、画像の詳細が再構築されていることがわかります。



次はNMFです。


流れとしては、

○顔画像への適用
○信号データへの適用

となっています。


○顔画像への適用
NMFを「Labeled Faces in the Wild」データセットに適用して、PCAとの違いを見てみます。

mglearn.plots.plot_nmf_faces(X_train, X_test, image_shape)

f:id:nakatatsu_com:20190915181806p:plain

# 最初の15主成分
from sklearn.decomposition import NMF

nmf = NMF(n_components=15, random_state=0)
nmf.fit(X_train)
X_train_nmf = nmf.transform(X_train)
X_test_nmf = nmf.transform(X_test)

fix, axes = plt.subplots(3, 5, figsize=(15, 12), subplot_kw={'xticks': (), 'yticks': ()})
for i, (component, ax) in enumerate(zip(nmf.components_, axes.ravel())):
    ax.imshow(component.reshape(image_shape))
    ax.set_title('{}.component'.format(i))

f:id:nakatatsu_com:20190915181833p:plain

# 成分3が大きい画像を表示
compn = 3
inds = np.argsort(X_train_nmf[:, compn])[::-1]
fig, axes = plt.subplots(2, 5, figsize=(15, 8), subplot_kw={'xticks': (), 'yticks': ()})
for i, (ind, ax) in enumerate(zip(inds, axes.ravel())):
    ax.imshow(X_train[ind].reshape(image_shape))

f:id:nakatatsu_com:20190915181845p:plain

# 成分7が大きい画像を表示
compn = 7
inds = np.argsort(X_train_nmf[:, compn])[::-1]
fig, axes = plt.subplots(2, 5, figsize=(15, 8), subplot_kw={'xticks': (), 'yticks': ()})
for i, (ind, ax) in enumerate(zip(inds, axes.ravel())):
    ax.imshow(X_train[ind].reshape(image_shape))

f:id:nakatatsu_com:20190915181900p:plain

抽出した成分によって元画像を再構築した結果は、PCAを用いたときと傾向は似ていますが、少し質が悪くなっています。
抽出された最初の15成分を見てみると、PCAよりもはるかに顔のプロトタイプを捉えていることがわかります。
これらは、PCAは再構築に最適な方向を見つける、NFMはデータから興味深いパターンを見つけるという特徴を裏付ける結果となっています。
また、成分3は少し右を向いた顔を、成分7は少し左を向いた画像を表していると予想できますので、これらの成分が特に強い画像をプロットしてみると、予想通りの結果が得られました。


○信号データへの適用
NMFはいくつもの独立した発生源から得られたデータを重ね合わせて作られるようなデータに対して特に有効です。ここでは、NMFを使用して3つの信号源からの信号が組み合わされた信号から元の成分を取り出してみます。

# NMFとPCAを使用した信号の復元
A = np.random.RandomState(0).uniform(size=(100, 3))
X = np.dot(S, A.T)
print('Shape of measurements: {}'.format(X.shape))
Shape of measurements: (2000, 100)
nmf = NMF(n_components=3, random_state=42)
S_ = nmf.fit_transform(X)
print('Recovered signal shape: {}'.format(S_.shape))
Recovered signal shape: (2000, 3)
pca = PCA(n_components=3)
H = pca.fit_transform(X)
models = [X, S, S_, H]
names = ['Observations (first three measurements)',
         'True sources',
         'NMF recovered signals',
         'PCA recovered signals']

fig, axes = plt.subplots(4, figsize=(8, 4), gridspec_kw={'hspace': 0.5},
                         subplot_kw={'xticks': (), 'yticks': ()})

for model, name, ax in zip(models, names, axes):
    ax.set_title(name)
    ax.plot(model[:, :3], '-')

f:id:nakatatsu_com:20190915183249p:plain

NMFでは元の信号源をうまく特定できていることがわかります。



最後はt-SNEです。


○数字データセットを可視化
t-SNEを手書きデータセットに適用して、PCAとの違いを見てみます。

# 手書き数字データセット
from sklearn.datasets import load_digits

digits = load_digits()

fig, axes = plt.subplots(2, 5, figsize=(10, 5), subplot_kw={'xticks': (), 'yticks': ()})
for ax, img in zip(axes.ravel(), digits.images):
    ax.imshow(img)

f:id:nakatatsu_com:20190915183954p:plain

# PCAを構築し、可視化
pca = PCA(n_components=2)
pca.fit(digits.data)
digits_pca = pca.transform(digits.data)
colors = ['#476A2A', '#7851B8', '#BD3430', '#4A2D4E', '#875525',
          '#A83683', '#4E655E', '#853541', '#3A3120', '#535D8E']
plt.figure(figsize=(10, 10))
plt.xlim(digits_pca[:, 0].min(), digits_pca[:, 0].max())
plt.ylim(digits_pca[:, 1].min(), digits_pca[:, 1].max())
for i in range(len(digits.data)):
    plt.text(digits_pca[i, 0], digits_pca[i, 1], str(digits.target[i]),
             color=colors[digits.target[i]],
             fontdict={'weight': 'bold', 'size': 9})
plt.xlabel('First principal component')
plt.ylabel('Second principal component');

f:id:nakatatsu_com:20190915183938p:plain

# t-SNEを構築し、可視化
from sklearn.manifold import TSNE

tsne = TSNE(random_state=42)
digits_tsne = tsne.fit_transform(digits.data)
plt.figure(figsize=(10, 10))
plt.xlim(digits_tsne[:, 0].min(), digits_tsne[:, 0].max() + 1)
plt.ylim(digits_tsne[:, 1].min(), digits_tsne[:, 1].max() + 1)
for i in range(len(digits.data)):
    plt.text(digits_tsne[i, 0], digits_tsne[i, 1], str(digits.target[i]),
             color=colors[digits.target[i]],
             fontdict={'weight': 'bold', 'size': 9})
plt.xlabel('t-SNE feature 0')
plt.ylabel('t-SNE feature 1');

f:id:nakatatsu_com:20190915183924p:plain

PCAではうまく分離できていないが、t-SNEではかなり明確に分類されていることがわかります。


まとめ

PCAは可視化や特徴量の抽出に広く用いられます。
NMFは主に特徴量の抽出に、t-SNEは主に可視化に用いられます。
特徴量の抽出について、PCAは再構築やデータのエンコードに、NFMはデータから興味深いパターンを見つけるのに用いられます。
可視化について、t-SNEはPCAよりはるかに複雑なマッピングを行い、より良い可視化を実現することができます。


参考文献

参考文献です。

Pythonではじめる機械学習 ―scikit-learnで学ぶ特徴量エンジニアリングと機械学習の基礎

Pythonではじめる機械学習 ―scikit-learnで学ぶ特徴量エンジニアリングと機械学習の基礎

教師あり機械学習(ニューラルネットワーク)

こんにちは。nakatatsuです。

今回も機械学習について学んでいきます。基本的にはPythonとscikit-learnライブラリを使用します。
参考書は
Pythonではじめる機械学習オライリー・ジャパン
を使用します。

今回は教師あり学習の「ニューラルネットワークアルゴリズムです。ここでは比較的簡単な多層パーセプトロン(以下、MLP)について見ていきます。

MLPは線形モデルを一般化し、決定までに複数のステージで計算するものと見ることができます。

下図のように、左側のノード群は入力特徴量、接続している線が学習された係数、右側のノードが出力を表しています。MLPでは中間処理ステップを表す隠れユニットの計算で重み付き和が行われ、次にこの隠れユニットの値に対して重み付き和が行われて最後の結果が出力されます。

f:id:nakatatsu_com:20190913132050p:plain

さらにこのモデルを強力にするために、個々の隠れユニットの重み付き和を計算したら、その結果に対して非線形関数を適用します。非線形関数としては、reluやtanhが用いられます。この非線形関数によって、より複雑なモデルを構築することができます。

f:id:nakatatsu_com:20190913134613p:plain

yを計算する式は以下のようになります。
wは入力xと隠れ層hの重み、vは隠れ層hと出力yの重みです。
この式により、クラス分類や回帰を行います。

h[0] = tanh(w[0, 0] × x[0] + w[1, 0] × x[1] + w[2, 0] × x[2] + w[3, 0] × x[3] + b[0])
h[1] = tanh(w[0, 1] × x[0] + w[1, 1] × x[1] + w[2, 1] × x[2] + w[3, 1] × x[3] + b[1])
h[2] = tanh(w[0, 2] × x[0] + w[1, 2] × x[1] + w[2, 2] × x[2] + w[3, 2] × x[3] + b[2])

y = v[0] × h[0] +v[1] × h[1] + v[2] × h[2] + b


種類

ニューラルネットワークによるクラス分類(クラスラベルの予測)

ニューラルネットワークによる回帰(連続値の予測)

回帰については、今回取り扱いません。


準備

Anacondaを使用するための準備についてはこのブログの過去記事「教師あり機械学習(k-最近傍法)」を参照してください。

nakatatsu-com.hatenablog.com


結果

AnacondaのJupyterLabを使用して作成しました。


流れとしては、

ニューラルネットワークの構築
○隠れユニットを変更してモデル構築
○隠れ層および非線形関数を変更してモデル構築
○alphaを変更してモデル構築
○乱数を変更してモデル構築
○特徴量のスケール変換を行い、精度を導出

となっています。


import mglearn
import matplotlib.pyplot as plt


ニューラルネットワークの構築
非線形関数relu、隠れ層1、隠れユニット100の条件でニューラルネットワークを構築してみます。

# ニューラルネットワークの構築(relu、隠れ層1、隠れユニット100)
from sklearn.neural_network import MLPClassifier
from sklearn.datasets import make_moons # two_moonsデータセット
from sklearn.model_selection import train_test_split

X, y = make_moons(n_samples=100, noise=0.25, random_state=3)

X_train, X_test, y_train, y_test = train_test_split(X, y, stratify=y, random_state=42)

mlp = MLPClassifier(solver='lbfgs', random_state=0).fit(X_train, y_train)
mglearn.plots.plot_2d_separator(mlp, X_train, fill=True, alpha=0.3)
mglearn.discrete_scatter(X_train[:, 0], X_train[:, 1], y_train)
plt.xlabel('Feature 0')
plt.ylabel('Feature 1');

f:id:nakatatsu_com:20190913141353p:plain

比較的滑らかな決定境界を学習していることがわかります。隠れユニット100はこの小さいデータセットに対しては明らかに大きすぎると予想されます。


○隠れユニットを変更してモデル構築
非線形関数relu、隠れ層1、隠れユニット10の条件でニューラルネットワークを構築してみます。

# ニューラルネットワークの構築(relu、隠れ層1、隠れユニット10)
mlp = MLPClassifier(solver='lbfgs', random_state=0, hidden_layer_sizes=[10])
mlp.fit(X_train, y_train)
mglearn.plots.plot_2d_separator(mlp, X_train, fill=True, alpha=0.3)
mglearn.discrete_scatter(X_train[:, 0], X_train[:, 1], y_train)
plt.xlabel('Feature 0')
plt.ylabel('Feature 1');

f:id:nakatatsu_com:20190913141835p:plain

隠れユニットを10にすると、決定境界は少しギザギザになっているので、隠れユニットを減らすと複雑さは低減することがわかります。


○隠れ層および非線形関数を変更してモデル構築
隠れ層を増やしニューラルネットワークを構築してみます。次に、非線形関数にtanhを用いてニューラルネットワークを構築してみます。

# ニューラルネットワークの構築(relu、隠れ層2、隠れユニット10)
mlp = MLPClassifier(solver='lbfgs', random_state=0, hidden_layer_sizes=[10, 10])
mlp.fit(X_train, y_train)
mglearn.plots.plot_2d_separator(mlp, X_train, fill=True, alpha=0.3)
mglearn.discrete_scatter(X_train[:, 0], X_train[:, 1], y_train)
plt.xlabel('Feature 0')
plt.ylabel('Feature 1');

f:id:nakatatsu_com:20190913142615p:plain

# ニューラルネットワークの構築(tanh、隠れ層2、隠れユニット10)
mlp = MLPClassifier(solver='lbfgs', activation='tanh',
                    random_state=0, hidden_layer_sizes=[10, 10])
mlp.fit(X_train, y_train)
mglearn.plots.plot_2d_separator(mlp, X_train, fill=True, alpha=0.3)
mglearn.discrete_scatter(X_train[:, 0], X_train[:, 1], y_train)
plt.xlabel('Feature 0')
plt.ylabel('Feature 1');

f:id:nakatatsu_com:20190913142628p:plain

隠れ層を2にすると、決定境界は少し滑らかになっているので、 隠れ層を増やすとより複雑なモデルになることがわかります。
非線形関数にtanhを用いると、決定境界は少し滑らかになっているので、非線形関数によって複雑さを制御できることがわかります。


○alphaを変更してモデル構築
複雑さを制御できるパラメータalphaを変更して、ニューラルネットワークを構築してみます。

# 異なる隠れユニットとalphaに対する決定境界
fig, axes = plt.subplots(2, 4, figsize=(20, 8))
for axx, n_hidden_nodes in zip(axes, [10, 100]):
    for ax, alpha in zip(axx, [0.0001, 0.01, 0.1, 1]):
        mlp = MLPClassifier(solver='lbfgs', random_state=0,
                            hidden_layer_sizes=[n_hidden_nodes, n_hidden_nodes],
                            alpha=alpha)
        mlp.fit(X_train, y_train)
        mglearn.plots.plot_2d_separator(mlp, X_train, fill=True, alpha=0.3, ax=ax)
        mglearn.discrete_scatter(X_train[:, 0], X_train[:, 1], y_train, ax=ax)
        ax.set_title('n_hidden=[{}, {}]\nalpha={:.4f}'.format(
            n_hidden_nodes, n_hidden_nodes, alpha))

f:id:nakatatsu_com:20190913143859p:plain

alphaを増やすと、決定境界は少しギザギザになっているので、複雑さは低減することがわかります。


○乱数を変更してモデル構築
乱数による初期化がモデルに与える影響を確認するために、乱数を変更してニューラルネットワークを構築してみます。

# 異なる乱数で初期化された決定境界
fig, axes = plt.subplots(2, 4, figsize=(20, 8))
for i, ax in enumerate(axes.ravel()):
    mlp = MLPClassifier(solver='lbfgs', random_state=i,
                        hidden_layer_sizes=[100, 100])
    mlp.fit(X_train, y_train)
    mglearn.plots.plot_2d_separator(mlp, X_train, fill=True, alpha=0.3, ax=ax)
    mglearn.discrete_scatter(X_train[:, 0], X_train[:, 1], y_train, ax=ax)

f:id:nakatatsu_com:20190913145658p:plain

同じパラメータを用いても、異なる乱数を用いると異なったモデルが得られることがわかります。これもニューラルネットワークの重要な性質の1つです。


○特徴量のスケール変換を行い、精度を導出
サポートベクタマシンと同様にニューラルネットワークもデータのスケールの影響を強く受けるので、スケール変換を行い、その後精度を導出してみます。

# 各特徴量の最大値
from sklearn.datasets import load_breast_cancer # cancerデータセット

cancer = load_breast_cancer()
print('Cancer data per-feature maxima:\n{}'.format(cancer.data.max(axis=0)))
Cancer data per-feature maxima:
[2.811e+01 3.928e+01 1.885e+02 2.501e+03 1.634e-01 3.454e-01 4.268e-01
 2.012e-01 3.040e-01 9.744e-02 2.873e+00 4.885e+00 2.198e+01 5.422e+02
 3.113e-02 1.354e-01 3.960e-01 5.279e-02 7.895e-02 2.984e-02 3.604e+01
 4.954e+01 2.512e+02 4.254e+03 2.226e-01 1.058e+00 1.252e+00 2.910e-01
 6.638e-01 2.075e-01]
# 精度
X_train, X_test, y_train, y_test = train_test_split(
    cancer.data, cancer.target, random_state=0)

mlp = MLPClassifier(random_state=42)
mlp.fit(X_train, y_train)

print('Accuracy on training set: {:.2f}'.format(mlp.score(X_train, y_train)))
print('Accuracy on test set: {:.2f}'.format(mlp.score(X_test, y_test)))
Accuracy on training set: 0.94
Accuracy on test set: 0.92
# スケール変換(平均0、分散1)し、精度を導出
mean_on_train = X_train.mean(axis=0)
std_on_train = X_train.std(axis=0)

X_train_scaled = (X_train - mean_on_train) / std_on_train
X_test_scaled = (X_test - mean_on_train) / std_on_train

mlp = MLPClassifier(random_state=0)
mlp.fit(X_train_scaled, y_train)

print('Accuracy on training set: {:.3f}'.format(mlp.score(X_train_scaled, y_train)))
print('Accuracy on test set: {:.3f}'.format(mlp.score(X_test_scaled, y_test)))
Accuracy on training set: 0.991
Accuracy on test set: 0.965
# 学習繰り返し回数を増やし、精度を導出
mlp = MLPClassifier(max_iter=1000, random_state=0)
mlp.fit(X_train_scaled, y_train)

print('Accuracy on training set: {:.3f}'.format(mlp.score(X_train_scaled, y_train)))
print('Accuracy on test set: {:.3f}'.format(mlp.score(X_test_scaled, y_test)))
Accuracy on training set: 1.000
Accuracy on test set: 0.972
# alphaを大きくし、精度を導出
mlp = MLPClassifier(max_iter=1000, alpha=1, random_state=0)
mlp.fit(X_train_scaled, y_train)

print('Accuracy on training set: {:.3f}'.format(mlp.score(X_train_scaled, y_train)))
print('Accuracy on test set: {:.3f}'.format(mlp.score(X_test_scaled, y_test)))
Accuracy on training set: 0.988
Accuracy on test set: 0.972

データのスケール変換前は、訓練セット精度が94%、テストセット精度が92%でした。
データのスケール変換を行い、再度精度を導出すると、訓練セット精度が99.1%、テストセット精度が96.5%と大きく改善したことがわかります。
さらに学習繰り返し回数を増やし、alphaを調整すると、訓練セット精度が98・8%、テストセット精度が97.2%となり、より汎化性能が上がり、良いモデルを構築することができました。


まとめ

以上より、ニューラルネットワークにはモデルの複雑さを制御する方法が、隠れ層の数、隠れユニットの数、alphaとたくさんあることがわかりました。
他にも、モデルは乱数による初期化やデータのスケールの影響を受けることがわかりました。


長所・短所

長所
・非常に複雑なモデルを構築できる
・特に大きなデータセットに有効
短所
・データのスケールを調整する必要がある
・パラメータを調整する必要がある
・大きいモデルは訓練に時間がかかる


参考文献

参考文献です。

Pythonではじめる機械学習 ―scikit-learnで学ぶ特徴量エンジニアリングと機械学習の基礎

Pythonではじめる機械学習 ―scikit-learnで学ぶ特徴量エンジニアリングと機械学習の基礎

教師あり機械学習(サポートベクタマシン)

こんにちは。nakatatsuです。

今回も機械学習について学んでいきます。基本的にはPythonとscikit-learnライブラリを使用します。
参考書は
Pythonではじめる機械学習オライリー・ジャパン
を使用します。

今回は教師あり学習の「サポートベクタマシン」アルゴリズムです。サポートベクタマシンは入力空間の超平面のような簡単なモデルではなく、より複雑なモデルの構築を可能にするために線形サポートベクタマシンを拡張したものです。

線形モデルは直線や超平面が柔軟性を制限するため、非常に制約が強くなっています。線形モデルを柔軟にする1つの方法は、新たな特徴量を追加することです。例えば、入力特徴量の交互作用や多項式項です。

サポートベクタマシンはクラス分類にも回帰にも利用できます。


種類

・サポートベクタマシンによるクラス分類(クラスラベルの予測)

・サポートベクタマシンによる回帰(連続値の予測)

回帰については、今回取り扱いません。


準備

Anacondaを使用するための準備についてはこのブログの過去記事「教師あり機械学習(k-最近傍法)」を参照してください。

nakatatsu-com.hatenablog.com


結果

AnacondaのJupyterLabを使用して作成しました。


流れとしては、

○線形モデルに新たな特徴量を追加
○サポートベクタマシンの構築
○パラメータごとに決定境界を描画し、精度を導出
○特徴量のスケール変換を行い、精度を導出

となっています。


import mglearn
import matplotlib.pyplot as plt
import numpy as np


○線形モデルに新たな特徴量を追加
線形モデルに新たな特徴量を追加してみます。ここでは、2番目の特徴量の2乗を加えてみます。

# ガウス分布でサンプリングした2次元データセット
from sklearn.datasets import make_blobs

X, y = make_blobs(centers=4, random_state=8)
y = y % 2
# 決定境界を描画
from sklearn.svm import LinearSVC

linear_svm = LinearSVC().fit(X, y)

mglearn.plots.plot_2d_separator(linear_svm, X)
mglearn.discrete_scatter(X[:, 0], X[:, 1], y)
plt.xlabel('Feature 0')
plt.ylabel('Feature 1');

f:id:nakatatsu_com:20190912150325p:plain

# 入力特徴量を拡張し、3Dで可視化
X_new = np.hstack([X, X[:, 1:] ** 2]) # 2番目の特徴量の2乗を追加

from mpl_toolkits.mplot3d import Axes3D, axes3d

figure = plt.figure()
ax = Axes3D(figure, elev=-152, azim=-26)
# 拡張された3次元空間で決定境界を描画
linear_svm_3d = LinearSVC().fit(X_new, y)
coef, intercept = linear_svm_3d.coef_.ravel(), linear_svm_3d.intercept_

figure = plt.figure()
ax = Axes3D(figure, elev=-152, azim=-26)
xx = np.linspace(X_new[:, 0].min() - 2, X_new[:, 0].max() + 2, 50)
yy = np.linspace(X_new[:, 1].min() - 2, X_new[:, 1].max() + 2, 50)

XX, YY = np.meshgrid(xx, yy)
ZZ = (coef[0] * XX + coef[1] * YY + intercept) / -coef[2]
ax.plot_surface(XX, YY, ZZ, rstride=8, cstride=8, alpha=0.3)

mask = y == 0
# y==0の点をプロット
ax.scatter(X_new[mask, 0], X_new[mask, 1], X_new[mask, 2], c='b', cmap=mglearn.cm2, s=60)
# y==1の点をプロット
ax.scatter(X_new[~mask, 0], X_new[~mask, 1], X_new[~mask, 2], c='r', marker='^', cmap=mglearn.cm2, s=60)

ax.set_xlabel('feature0')
ax.set_ylabel('feature1')
ax.set_zlabel('feature1 ** 2');

f:id:nakatatsu_com:20190912150858p:plain

# 元の特徴量で決定境界を描画
ZZ = YY ** 2
dec = linear_svm_3d.decision_function(np.c_[XX.ravel(), YY.ravel(), ZZ.ravel()])
plt.contourf(XX, YY, dec.reshape(XX.shape), levels=[dec.min(), 0, dec.max()], cmap=mglearn.cm2, alpha=0.5)
mglearn.discrete_scatter(X[:, 0], X[:, 1], y)
plt.xlabel('Feature 0')
plt.ylabel('Feature 1');

f:id:nakatatsu_com:20190912150930p:plain

線形モデルによるクラス分類では、直線で分離することしかできないのでうまく分離できていないことがわかります。
新たな特徴量を追加し3Dで可視化してみると、うまく分離できていることがわかります。
元の特徴量で決定境界を描画してみると、直線ではなく楕円に近くなっており線形ではなくなっていることがわかります。
このことから、非線形の特徴量を加えることで線形モデルが柔軟になり、より複雑なモデルを構築することが可能になると言えます。


○サポートベクタマシンの構築
入力特徴量を拡張することができるRBFカーネル法を用いて、サポートベクタマシンを構築してみます。

# RBFカーネル法による決定境界とサポートベクタ
from sklearn.svm import SVC

X, y = mglearn.tools.make_handcrafted_dataset()
svm = SVC(kernel='rbf', C=10, gamma=0.1).fit(X, y)
mglearn.plots.plot_2d_separator(svm, X, eps=0.5)
mglearn.discrete_scatter(X[:, 0], X[:, 1], y)

# サポートベクタをプロット
sv = svm.support_vectors_
sv_labels = svm.dual_coef_.ravel() > 0
mglearn.discrete_scatter(sv[:, 0], sv[:, 1], sv_labels, s=15, markeredgewidth=3)
plt.xlabel('Feature 0')
plt.ylabel('Feature 1');

f:id:nakatatsu_com:20190912152735p:plain

サポートベクタマシンによる境界は、滑らかで非線形になっていることがわかります。
縁取りされた点はサポートベクタと呼ばれており、決定境界を決定するために用いられます。


○パラメータごとに決定境界を描画し、精度を導出
パラメータを変化させてそれぞれの決定境界を描画してみます。その後、精度を導出してみます。
主なパラメータとしては、gammaとCがあります。

・gamma
点が近いということを意味するスケールを決定します。gammaが大きくなると、モデルは複雑になります。
・C
個々のデータポイントの重要度を制限します。Cが大きくなると、モデルは複雑になります。

# パラメータの調整
fig, axes = plt.subplots(3, 3, figsize=(15, 10))

for ax, C in zip(axes, [-1, 0, 3]):
    for a, gamma in zip(ax, range(-1, 2)):
        mglearn.plots.plot_svm(log_C=C, log_gamma=gamma, ax=a)
        
axes[0, 0].legend(['class 0', 'class 1', 'sv class 0', 'sv class 1'], ncol=4, loc=(0.9, 1.2));

f:id:nakatatsu_com:20190912153957p:plain

# 精度
from sklearn.datasets import load_breast_cancer # cancerデータセット
from sklearn.model_selection import train_test_split

cancer = load_breast_cancer()
X_train, X_test, y_train, y_test = train_test_split(
    cancer.data, cancer.target, random_state=0)

svc = SVC()
svc.fit(X_train, y_train)

print('Accuracy on training set: {:.2f}'.format(svc.score(X_train, y_train)))
print('Accuracy on test set: {:.2f}'.format(svc.score(X_test, y_test)))
Accuracy on training set: 1.00
Accuracy on test set: 0.63

gammaを0.1から10に変化させています。gammaが小さい場合、決定境界はゆっくりとしか変化せず、モデルの複雑さは小さくなります。gammaが大きい場合、個々のデータポイントをより重視するようになり、モデルはより複雑になります。
Cを0.1から1000に変化させています。Cに関してもgammaと同様のことが言えます。
デフォルトパラメータ(C=1、gamma=1/n_features)で精度を導出してみると、訓練セット精度100%、テストセット精度63%となり、強く過剰適合していることがわかります。これは特徴量のスケールが大きく異なることが主な原因であると考えられます。


○特徴量のスケール変換を行い、精度を導出
それぞれのスケールが大体同じになるように変換します。その後、精度を導出してみます。

# 特徴量のレンジ
plt.plot(X_train.min(axis=0), 'o', label='min')
plt.plot(X_train.max(axis=0), '^', label='max')
plt.legend(loc=4)
plt.xlabel('Feature index')
plt.ylabel('Feature magnitude')
plt.yscale('log')

f:id:nakatatsu_com:20190912155431p:plain

# データの前処理(訓練セット)
min_on_training = X_train.min(axis=0)
range_on_training = (X_train - min_on_training).max(axis=0)

X_train_scaled = (X_train - min_on_training) / range_on_training
print('Minimum for each feature\n{}'.format(X_train_scaled.min(axis=0)))
print('Maximum for each feature\n{}'.format(X_train_scaled.max(axis=0)))
Minimum for each feature
[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
 0. 0. 0. 0. 0. 0.]
Maximum for each feature
[1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
 1. 1. 1. 1. 1. 1.]
# データの前処理(テストセット)
X_test_scaled = (X_test - min_on_training) / range_on_training

svc = SVC()
svc.fit(X_train_scaled, y_train)

print('Accuracy on training set: {:.3f}'.format(svc.score(X_train_scaled, y_train)))
print('Accuracy on test set: {:.3f}'.format(svc.score(X_test_scaled, y_test)))
Accuracy on training set: 0.948
Accuracy on test set: 0.951
# パラメータの調整
svc = SVC(C=1000)
svc.fit(X_train_scaled, y_train)

print('Accuracy on training set: {:.3f}'.format(svc.score(X_train_scaled, y_train)))
print('Accuracy on test set: {:.3f}'.format(svc.score(X_test_scaled, y_test)))
Accuracy on training set: 0.988
Accuracy on test set: 0.972

特徴量のスケールを確認してみると、桁違いにスケールが異なることが確認できました。
すべての特徴量を0から1の間になるように変換し精度を導出してみると、訓練セット精度94.8%、テストセット精度95.1%と結果が一変しました。この結果から適合不足と考えられるので、パラメータを調整します。
パラメータを調整し、再度精度を導出してみると、訓練セット精度98.8%、テストセット精度97.2%と大きく改良されました。


長所・短所

長所
・同じような意味を持つ特徴量からなる中規模なデータセットに対しては強力
短所
・データのスケールを調整する必要がある
・パラメータを調整する必要がある


参考文献

参考文献です。

Pythonではじめる機械学習 ―scikit-learnで学ぶ特徴量エンジニアリングと機械学習の基礎

Pythonではじめる機械学習 ―scikit-learnで学ぶ特徴量エンジニアリングと機械学習の基礎

教師あり機械学習(勾配ブースティング回帰木)

こんにちは。nakatatsuです。

今回も機械学習について学んでいきます。基本的にはPythonとscikit-learnライブラリを使用します。
参考書は
Pythonではじめる機械学習オライリー・ジャパン
を使用します。

今回は教師あり学習の「勾配ブースティング回帰木」アルゴリズムです。勾配ブースティング回帰木は決定木のアンサンブル法(複数の機械学習モデルを組み合わせることでより強力なモデルを構築する方法)です。決定木の利点の多くを残したまま、決定木の欠点の一部を補うことができます。

前々回の記事で述べたように、決定木の最大の問題点は訓練データに対して過剰適合してしまうことです。
勾配ブースティング回帰木は、1つ前の決定木の誤りを次の決定木が修正するようにして、決定木を順番に作成していきます。深さ1から5ぐらいの非常に浅い決定木が用いられ、これを多数組み合わせていきます。それぞれの決定木はデータの一部に対してしか良い予測を行えないので、繰り返し追加していくことで性能を向上させます。逆に言うと、これらを制御することで過剰適合を低減させることができます。
クラス分類でも回帰でも広く使われている機械学習手法です。


種類

・勾配ブースティング回帰木によるクラス分類(クラスラベルの予測)

・勾配ブースティング回帰木による回帰(連続値の予測)

回帰については、今回取り扱いません。


準備

Anacondaを使用するための準備についてはこのブログの過去記事「教師あり機械学習(k-最近傍法)」を参照してください。

nakatatsu-com.hatenablog.com


結果

AnacondaのJupyterLabを使用して作成しました。


流れとしては、

○勾配ブースティング回帰木の構築および精度の導出
○パラメータの調整
○特徴量の重要度の可視化

となっています。


import matplotlib.pyplot as plt
import numpy as np


○勾配ブースティング回帰木の構築および精度の導出
cancerデータセットを用いて、勾配ブースティング回帰木を構築し、精度を導出してみます。
デフォルト値の、深さの最大値3、決定木の数は100個、学習率は0.1のパラメータで行います。これらはそれぞれ、max_depth、n_estimators、learning_rateのパラメータで指定できます。

・n_estimators
n_estimatorsを増やすことで、訓練セットに対する過ちを補正する機会が増えるので、モデルは複雑になります。
・learning_rate
learning_rateを増やすことで、それまでの決定木の過ちをより強く補正しようとするので、モデルは複雑になります。

# max_depth=3, n_estimators=100, learning_rate=0.1
from sklearn.datasets import load_breast_cancer # cancerデータセット
from sklearn.model_selection import train_test_split
from sklearn.ensemble import GradientBoostingClassifier

cancer = load_breast_cancer()
X_train, X_test, y_train, y_test = train_test_split(
    cancer.data, cancer.target, random_state=0)

gbrt = GradientBoostingClassifier(random_state=0)
gbrt.fit(X_train, y_train)

print('Accuracy on training set: {:.3f}'.format(gbrt.score(X_train, y_train)))
print('Accuracy on test set: {:.3f}'.format(gbrt.score(X_test, y_test)))
Accuracy on training set: 1.000
Accuracy on test set: 0.958

訓練セットに対する精度が100%になっているので、おそらくは過剰適合しています。
深さの最大値を制限するか、学習率を下げることで、過剰適合を低減することができます。


○パラメータの調整
深さの最大値を制限した場合と学習率を下げた場合の両方で精度がどのように変化するか見てみます。

# 深さの最大値を制限した場合
# max_depth=1, n_estimators=100, learning_rate=0.1
gbrt = GradientBoostingClassifier(random_state=0, max_depth=1)
gbrt.fit(X_train, y_train)

print('Accuracy on training set: {:.3f}'.format(gbrt.score(X_train, y_train)))
print('Accuracy on test set: {:.3f}'.format(gbrt.score(X_test, y_test)))
Accuracy on training set: 0.991
Accuracy on test set: 0.972
# 学習率を下げた場合
# max_depth=3, n_estimators=100, learning_rate=0.01
gbrt = GradientBoostingClassifier(random_state=0, learning_rate=0.01)
gbrt.fit(X_train, y_train)

print('Accuracy on training set: {:.3f}'.format(gbrt.score(X_train, y_train)))
print('Accuracy on test set: {:.3f}'.format(gbrt.score(X_test, y_test)))
Accuracy on training set: 0.988
Accuracy on test set: 0.965

深さの最大値を制限した場合、学習率を下げた場合のどちらも訓練セットに対する精度は下がり、テストセットに対する精度は向上していることがわかります。すなわち、過剰適合を低減することができています。
また、学習率を下げた場合より深さの最大値を制限した場合の方がモデル性能は大きく向上していることがわかります。


○特徴量の重要度の可視化
特徴量の重要度を可視化して、モデルの詳細を見てみます。

# 特徴量の重要度
def plot_feature_importances_cancer(model):
    n_features = cancer.data.shape[1]
    plt.barh(range(n_features), model.feature_importances_, align='center') # 横棒グラフ
    plt.yticks(np.arange(n_features), cancer.feature_names)
    plt.xlabel('Feature importance')
    plt.ylabel('Feature')

gbrt = GradientBoostingClassifier(random_state=0, max_depth=1)
gbrt.fit(X_train, y_train)

plot_feature_importances_cancer(gbrt)

f:id:nakatatsu_com:20190911125643p:plain

ランダムフォレストの場合に似ていますが、こちらはいくつかの特徴量が完全に無視されていることがわかります。
ランダムフォレストを先に試し、予測時間が非常に重要な場合や機械学習モデルから最後の1%まで性能を絞り出したい場合には、勾配ブースティング回帰木を試してみるのが一般的です。


長所・短所

長所
・ランダムフォレストより精度が高い
・ランダムフォレストより予測が早い
・ランダムフォレストよりメモリ使用量が小さい
・データのスケールを考慮する必要がない
短所
・ランダムフォレストより訓練に時間がかかる
・パラメータのチューニングに細心の注意が必要である
・高次元の疎なデータには適さない


参考文献

参考文献です。

Pythonではじめる機械学習 ―scikit-learnで学ぶ特徴量エンジニアリングと機械学習の基礎

Pythonではじめる機械学習 ―scikit-learnで学ぶ特徴量エンジニアリングと機械学習の基礎

教師あり機械学習(ランダムフォレスト)

こんにちは。nakatatsuです。

今回も機械学習について学んでいきます。基本的にはPythonとscikit-learnライブラリを使用します。
参考書は
Pythonではじめる機械学習オライリー・ジャパン
を使用します。

今回は教師あり学習の「ランダムフォレスト」アルゴリズムです。ランダムフォレストは決定木のアンサンブル法(複数の機械学習モデルを組み合わせることでより強力なモデルを構築する方法)です。決定木の利点の多くを残したまま、決定木の欠点の一部を補うことができます。

前回の記事で述べたように、決定木の最大の問題点は訓練データに対して過剰適合してしまうことです。
ランダムフォレストは、少しずつ異なる決定木をたくさん集め平均を取ることで、過剰適合の度合いを減らすことができます。クラス分類でも回帰でも、現在最も広く使われている機械学習手法です。


種類

・ランダムフォレストによるクラス分類(クラスラベルの予測)

・ランダムフォレストによる回帰(連続値の予測)

回帰については、今回取り扱いません。


準備

Anacondaを使用するための準備についてはこのブログの過去記事「教師あり機械学習(k-最近傍法)」を参照してください。

nakatatsu-com.hatenablog.com


結果

AnacondaのJupyterLabを使用して作成しました。


流れとしては、

○ランダムフォレストの構築
○精度の導出
○特徴量の重要度の可視化

となっています。


import matplotlib.pyplot as plt
import mglearn
import numpy as np


○ランダムフォレストの構築
two_moonsデータセットを用いて、ランダムフォレストを構築してみます。決定木の数は5個とします。
ブートストラップサンプリングによって、ランダムフォレストの中の個々の決定木が少しずつ違うデータセットに対して構築されることになります。今回はそのパラメータであるn_samplesを100個としました。
さらに、個々のノードでの特徴量の選択によってそれぞれの決定木は異なる特徴量のサブセットに対して分割を行うことになります。今回はそのパラメータであるmax_featuresをデフォルト値としました。
これらの機構が組み合わされることでランダルフォレスト中の個々の決定木が異なるものとなります。

# ランダムフォレストの構築
from sklearn.datasets import make_moons # two_moonsデータセット
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestClassifier

X, y = make_moons(n_samples=100, noise=0.25, random_state=3) # n_samples:ブートストラップサンプリング数
X_train, X_test, y_train, y_test = train_test_split(X, y, stratify=y, random_state=42)

forest = RandomForestClassifier(n_estimators=5, random_state=2) # n_estimators:決定木数
forest.fit(X_train, y_train);
# プロット
fig, axes = plt.subplots(2, 3, figsize=(20, 10))
for i, (ax, tree) in enumerate(zip(axes.ravel(), forest.estimators_)): # axes.raverl():配列を1次元化
    ax.set_title('Tree {}'.format(i))
    mglearn.plots.plot_tree_partition(X_train, y_train, tree, ax=ax)
    
mglearn.plots.plot_2d_separator(forest, X_train, fill=True, ax=axes[-1, -1], alpha=0.4)
axes[-1, -1].set_title('Random Forest')
mglearn.discrete_scatter(X_train[:, 0], X_train[:, 1], y_train);

f:id:nakatatsu_com:20190910124855p:plain

5つの決定木が学習した決定境界は明らかに相互に異なっていることがわかります。ランダムフォレストは個々のどの決定木よりも過剰適合が少なく、直観に合致した決定境界を描いています。実際のアプリケーションは数百から数千にもなる決定木を使うので、決定境界はさらに滑らかになります。


○精度の導出
cancersデータセットを用いて、100個の決定木を用いたランダムフォレストを構築し、精度を見てみます。

# 精度
from sklearn.datasets import load_breast_cancer # cancerデータセット

cancer = load_breast_cancer()
X_train, X_test, y_train, y_test = train_test_split(
    cancer.data, cancer.target, random_state=0)
forest = RandomForestClassifier(n_estimators=100, random_state=0)
forest.fit(X_train, y_train)

print('Accuracy on training set: {:.3f}'.format(forest.score(X_train, y_train)))
print('Accuracy on test set: {:.3f}'.format(forest.score(X_test, y_test)))
Accuracy on training set: 1.000
Accuracy on test set: 0.972

パラメータを全く調整していないにも関わらず、97%の精度を示しています。これは線形モデルや個別の決定木よりも高い数字となっています。


○特徴量の重要度の可視化
ランダムフォレストでも特徴量の重要度を見ることができます。これは個々の決定木の特徴量の重要度を平均したものとなっています。多くの場合、ランダムフォレストによる特徴量の重要度は個々の決定木のそれよりも信頼できます。

# 特徴量の重要度
def plot_feature_importances_cancer(model):
    n_features = cancer.data.shape[1]
    plt.barh(range(n_features), model.feature_importances_, align='center') # 横棒グラフ
    plt.yticks(np.arange(n_features), cancer.feature_names)
    plt.xlabel('Feature importance')
    plt.ylabel('Feature')

plot_feature_importances_cancer(forest)

f:id:nakatatsu_com:20190910125944p:plain

ランダムフォレストでは決定木の場合よりもはるかに多くの特徴量に対して0以上の重要度を与えています。1つの決定木の場合と同じように、「worst radius」特徴量に高い重要度を与えているが、全体としては「worst perimeter」特徴量が最も重要となっています。
ランダムフォレストの結果は個々の決定木の結果よりも広い視野で見た全体像を捉えることができます。


長所・短所

長所
・単一の決定木よりも高速
・頑健(イレギュラーに強い) ・それほどパラメータチューニングをせずに使える
・データのスケールを考慮する必要がない
短所
・高次元の疎なデータには適さない


参考文献

参考文献です。

Pythonではじめる機械学習 ―scikit-learnで学ぶ特徴量エンジニアリングと機械学習の基礎

Pythonではじめる機械学習 ―scikit-learnで学ぶ特徴量エンジニアリングと機械学習の基礎

教師あり機械学習(決定木)

こんにちは。nakatatsuです。

今回も機械学習について学んでいきます。基本的にはPythonとscikit-learnライブラリを使用します。
参考書は
Pythonではじめる機械学習オライリー・ジャパン
を使用します。

今回は教師あり学習の「決定木」アルゴリズムです。決定木ではYes/Noで答えられる質問で構成された階層的な木構造を学習します。クラス分類と回帰に広く用いられています。

例えば、4種類の動物をなるべく少ない質問で区別したいとして、以下の例では、4クラスの動物(Hawk、Penguin、Dolphin、Bear)を3つの特徴量(Has feathers?、Can fly?、Has fins?)で識別するモデルを作成したことになります。このようなモデルを手で作成するのではなく、データから教師あり学習によって作成することができます。
f:id:nakatatsu_com:20190909090214p:plain


種類

・決定木によるクラス分類(クラスラベルの予測)

・決定木による回帰(連続値の予測)


準備

Anacondaを使用するための準備についてはこのブログの過去記事「教師あり機械学習(k-最近傍法)」を参照してください。

nakatatsu-com.hatenablog.com


結果

AnacondaのJupyterLabを使用して作成しました。


まずはクラス分類です。

流れとしては、

○決定木の構築
○複雑さの制御
○決定木の可視化
○特徴量の重要度の可視化

となっています。


import matplotlib.pyplot as plt
import numpy as np


○決定木の構築
two_moonsデータセットを用いて、深さ1、2、9の決定木を構築してみます。

# 深さ1、2、9の決定境界(左)と決定木(右)  
# two_moonsデータセット
mglearn.plots.plot_tree_progressive()

f:id:nakatatsu_com:20190909091909p:plain

f:id:nakatatsu_com:20190909091926p:plain

f:id:nakatatsu_com:20190909091939p:plain

f:id:nakatatsu_com:20190909091951p:plain

最初のレベルで2つのクラスがかなりきれいに分離できていることがわかります。
データの再帰分割は対象の領域に1つの対象値しか含まれなくなるまで繰り返します。すると、モデルは複雑になり過ぎて、訓練データに対して過剰適合してしまいます。深さ9では過剰適合が起きていることがわかります。


○複雑さの制御
過剰適合を防ぐための戦略として、構築過程で木の生成を早めに止める「事前枝刈り」があります。ここでは、事前枝刈りの1つである「木の深さを制御する方法」をcancerデータセットに適用して複雑さの制御を行います。

# 過剰適合
from sklearn.tree import DecisionTreeClassifier
from sklearn.datasets import load_breast_cancer # cancerデータセット
from sklearn.model_selection import train_test_split

cancer = load_breast_cancer()
X_train, X_test, y_train, y_test = train_test_split(
    cancer.data, cancer.target, stratify=cancer.target, random_state=42)
tree = DecisionTreeClassifier(random_state=0)
tree.fit(X_train, y_train)
print('Accuracy on training set: {:.3f}'.format(tree.score(X_train, y_train)))
print('Accuracy on test set: {:.3f}'.format(tree.score(X_test, y_test)))
Accuracy on training set: 1.000
Accuracy on test set: 0.937
# 事前枝刈り(max_depth=4)
tree = DecisionTreeClassifier(max_depth=4, random_state=0)
tree.fit(X_train, y_train)

print('Accuracy on training set: {:.3f}'.format(tree.score(X_train, y_train)))
print('Accuracy on test set: {:.3f}'.format(tree.score(X_test, y_test)))
Accuracy on training set: 0.988
Accuracy on test set: 0.951

木の深さを制限することで過剰適合が抑制され、訓練セットに対する精度は下がりますが、テストセットに対する制度は向上しました。


○決定木の可視化
決定木を可視化してみます。

# dotファイル形式で書き出す
from sklearn.tree import export_graphviz
export_graphviz(tree, out_file='tree.dot', class_names=['malignant', 'benign'],
                feature_names=cancer.feature_names, impurity=False, filled=True)
# graphvizモジュールで可視化
import graphviz

with open('tree.dot') as f:
    dot_graph = f.read()
graphviz.Source(dot_graph)

f:id:nakatatsu_com:20190909100519p:plain

最初の分岐には「worst radius」という特徴量が使われており、このレベルで大多数は分離できていることがわかる。


○特徴量の重要度の可視化
決定木を構築する際、アルゴリズムは重要な特徴量を選択し、分岐に使用します。「特徴量の重要度」と呼ばれる個々の特徴量がどの程度重要かを示す割合を見てみます。0から1の間の数で、0は「まったく使われていない」、1は「完全にターゲットを予想できる」を意味し、その和は1になります。

# 特徴量の重要度
print('Feature importances:\n{}'.format(tree.feature_importances_))
Feature importances:
[0.         0.         0.         0.         0.         0.
 0.         0.         0.         0.         0.01019737 0.04839825
 0.         0.         0.0024156  0.         0.         0.
 0.         0.         0.72682851 0.0458159  0.         0.
 0.0141577  0.         0.018188   0.1221132  0.01188548 0.        ]
# 特徴量の重要度を可視化
def plot_feature_importances_cancer(model):
    n_features = cancer.data.shape[1]
    plt.barh(range(n_features), model.feature_importances_, align='center') # 横棒グラフ
    plt.yticks(np.arange(n_features), cancer.feature_names)
    plt.xlabel('Feature importance')
    plt.ylabel('Feature')
    
plot_feature_importances_cancer(tree)

f:id:nakatatsu_com:20190909101415p:plain

最初の分岐に用いた特徴量「worst radius」が群を抜いて重要だということがわかります。これは最初のレベルで大多数は分離できているという観察結果と一致します。




次は回帰です。


import pandas as pd
import mglearn
import matplotlib.pyplot as plt
import numpy as np


○訓練データレンジ外側の予測
メモリ(RAM)価格の履歴データセットを使用して、訓練データレンジの外側に対して予測を行います。決定木と線形モデルを比較してみます。

# メモリ(RAM)価格の歴史的推移
import os
ram_prices = pd.read_csv(os.path.join(mglearn.datasets.DATA_PATH, 'ram_price.csv'))

plt.semilogy(ram_prices.date, ram_prices.price) # 片対数プロット
plt.xlabel('Year')
plt.ylabel('Price in $/Mbyte');

f:id:nakatatsu_com:20190909105347p:plain

# 2000年以降を予測
from sklearn.tree import DecisionTreeRegressor
from sklearn.linear_model import LinearRegression
data_train = ram_prices[ram_prices.date < 2000]
data_test = ram_prices[ram_prices.date >= 2000]

X_train = data_train.date[:, np.newaxis]
y_train = np.log(data_train.price)

tree = DecisionTreeRegressor().fit(X_train, y_train)
linear_reg = LinearRegression().fit(X_train, y_train)

X_all = ram_prices.date[:, np.newaxis]

pred_tree = tree.predict(X_all)
pred_lr = linear_reg.predict(X_all)

price_tree = np.exp(pred_tree)
price_lr = np.exp(pred_lr)
# 片対数プロット
plt.semilogy(data_train.date, data_train.price, label='Training data')
plt.semilogy(data_test.date, data_test.price, label='Test data')
plt.semilogy(ram_prices.date, price_tree, label='Tree prediction')
plt.semilogy(ram_prices.date, price_lr, label='Linear prediction')
plt.legend();

f:id:nakatatsu_com:20190909110007p:plain

訓練データレンジの外側に対して、以下の結果となりました。

・線形モデル
データを直線で近似し、訓練データとテストデータの双方において細かい変異は取りこぼしているものの、テストデータに対して良い予測を与えています。

・決定木
訓練データに対しては、複雑さを制御していないので完全な予測(グラフが重なっている)を行なっています。しかし、訓練データレンジの外側に対しては、最後の点を返すだけになっており、予測することができていません。

このように決定木は訓練データにない領域に関しては新しい答えを生成することができません。


長所・短所

長所
・非常に高速
・結果のモデルが容易に可視化可能
・データのスケールを考慮する必要がない
短所
・過剰適合しやすく、汎化性能が低い傾向がある


参考文献

参考文献です。

Pythonではじめる機械学習 ―scikit-learnで学ぶ特徴量エンジニアリングと機械学習の基礎

Pythonではじめる機械学習 ―scikit-learnで学ぶ特徴量エンジニアリングと機械学習の基礎