SVM～カーネル法

2020-06-24 / tau / コメントする

概要

書籍”Pythonではじめる機械学習”の2.3.7 カーネル法を用いた”サポートベクタマシン”の写経

線形特徴量の非線形化

線形モデルでは分離不可能なデータ

以下は、scikit-learnのmake_blobs()により生成した2特徴量、2クラス分類のデータに線形サポートベクターマシンを適用した例。このとき、決定境界は以下のように得られる。

(1) $\begin{gather*} b + w_0 f_0 + w_1 f_1 = 0\\ b \approx -0.2817,\; w_0 \approx 0.1261,\; w_1 \approx -0.0918 \end{gather*}$

決定境界より上側では多項式の値は負となり、下側では正となるが、この境界は明らかに2つのクラスを分割していない。このように単純な例でも、線形モデルでは的確なクラス分類はできない。

以下のコードでは、原典と以下が異なっている。

収束しないという警告を受けて、LinearSVCのmax_iterをデフォルトの1000より大きな値としている

import matplotlib.pyplot as plt
from sklearn.datasets import make_blobs
from sklearn.svm import LinearSVC

X, y = make_blobs(centers=4, random_state=8)
y = y % 2

linsvm = LinearSVC(max_iter=5500).fit(X, y)

y_min, y_max = -15, 15
b = linsvm.intercept_
w0 = linsvm.coef_[0][0]
w1 = linsvm.coef_[0][1]
x_lower = -(b + w1 * y_min) / w0
x_upper = -(b + w1 * y_max) / w0

fig, ax = plt.subplots()
X0 = X[y==0]
X1 = X[y==1]
ax.scatter(X0[:, 0], X0[:, 1], marker='o', s=60, ec='k')
ax.scatter(X1[:, 0], X1[:, 1], marker='^', s=60, ec='k')
ax.plot([x_lower, x_upper], [y_min, y_max], linewidth=2, c='tab:green')
ax.set_ylim(y_min, y_max)
ax.set_xlabel("Feature-0")
ax.set_ylabel("Feature-1")
plt.show()

import matplotlib.pyplot as plt

from sklearn.datasets import make_blobs

from sklearn.svm import LinearSVC

X, y = make_blobs(centers=4, random_state=8)

y = y % 2

linsvm = LinearSVC(max_iter=5500).fit(X, y)

y_min, y_max = -15, 15

b = linsvm.intercept_

w0 = linsvm.coef_[0][0]

w1 = linsvm.coef_[0][1]

x_lower = -(b + w1 * y_min) / w0

x_upper = -(b + w1 * y_max) / w0

fig, ax = plt.subplots()

X0 = X[y==0]

X1 = X[y==1]

ax.scatter(X0[:, 0], X0[:, 1], marker='o', s=60, ec='k')

ax.scatter(X1[:, 0], X1[:, 1], marker='^', s=60, ec='k')

ax.plot([x_lower, x_upper], [y_min, y_max], linewidth=2, c='tab:green')

ax.set_ylim(y_min, y_max)

ax.set_xlabel("Feature-0")

ax.set_ylabel("Feature-1")

plt.show()

非線形特徴量の追加

ここで、特徴量1の2乗を新たな特徴量として加える。この場合、3つの特徴量に対して3次元空間内に各点が位置し、それぞれがクラス0/1に属している。新たな特徴量の追加によって、その軸の方向に各点が立ち上がり、真ん中の三角形の点群と両側の丸印の点群が平面でうまく分割できそうである。

このデータセットに対して、線形SVMを適用し、決定境界を描いたのが以下の画像。特徴量が2つの場合の決定境界は直線だったが、特徴が3つになると決定境界は平面となる。予想通り、単純な平面で2つのクラスが分けられている。この決定境界の式は以下のようになる。

(2) $\begin{gather*} b + w_0 f_0 + w_1 f_1 + w_2 {f_1}^2 = 0\\ b \approx 1.1734,\; w_0 \approx 0.1301,\; w_1 \approx -0.2203,\; w_2 = -0.0597 \end{gather*}$

この平面に対して上側（f₁²が小さい側）では多項式の値は正となり、その反対側では負となる。

以下のコードは、原典と以下が異なっている。

収束しないという警告を受けて、LinearSVCのmax_iterをデフォルトの1000より大きな値としている
Axes3Dの生成の仕方を最新のバージョンに合ったものとしている
- 原典ではFigureオブジェクトを生成し、それをビューに関する引数とともにAxes3Dコンストラクターに渡している
- 本コードでは、subplotsの引数でprojectionを指定してFigureとAxes3Dを同時に生成し、veiw_init()でビューに関する引数を指定

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import make_blobs
from sklearn.svm import LinearSVC
from mpl_toolkits.mplot3d import Axes3D

X, y = make_blobs(centers=4, random_state=8)
y = y % 2
X_new = np.hstack((X, X[:, 1].reshape(-1, 1)**2))
X0, X1 = X_new[y==0], X_new[y==1]

linsvc = LinearSVC(max_iter=3700).fit(X_new, y)
intercept = linsvc.intercept_
coef = linsvc.coef_.ravel()

fig, ax = plt.subplots(subplot_kw=dict(projection='3d'))
ax.view_init(elev=-152, azim=-23)

u = np.linspace(X_new[:, 0].min() - 2, X_new[:, 1].max() + 2)
v = np.linspace(X_new[:, 0].min() - 2, X_new[:, 1].max() + 2)
u, v = np.meshgrid(u, v)
w = -(coef[0] * u + coef[1] * v + intercept) / coef[2]

ax.scatter(X0[:, 0], X0[:, 1], X0[:, 2], marker='o', s=40, ec='k')
ax.scatter(X1[:, 0], X1[:, 1], X1[:, 2], marker='^', s=40, ec='k')
ax.plot_wireframe(u, v, w, rstride=8, cstride=8, color='tab:green', alpha=0.5)

ax.set_xlabel("Feature-0")
ax.set_ylabel("Feature-1")
ax.set_zlabel("Feature-1**2")

plt.show()

import numpy as np

import matplotlib.pyplot as plt

from sklearn.datasets import make_blobs

from sklearn.svm import LinearSVC

from mpl_toolkits.mplot3d import Axes3D

X, y = make_blobs(centers=4, random_state=8)

y = y % 2

X_new = np.hstack((X, X[:, 1].reshape(-1, 1)**2))

X0, X1 = X_new[y==0], X_new[y==1]

linsvc = LinearSVC(max_iter=3700).fit(X_new, y)

intercept = linsvc.intercept_

coef = linsvc.coef_.ravel()

fig, ax = plt.subplots(subplot_kw=dict(projection='3d'))

ax.view_init(elev=-152, azim=-23)

u = np.linspace(X_new[:, 0].min() - 2, X_new[:, 1].max() + 2)

v = np.linspace(X_new[:, 0].min() - 2, X_new[:, 1].max() + 2)

u, v = np.meshgrid(u, v)

w = -(coef[0] * u + coef[1] * v + intercept) / coef[2]

ax.scatter(X0[:, 0], X0[:, 1], X0[:, 2], marker='o', s=40, ec='k')

ax.scatter(X1[:, 0], X1[:, 1], X1[:, 2], marker='^', s=40, ec='k')

ax.plot_wireframe(u, v, w, rstride=8, cstride=8, color='tab:green', alpha=0.5)

ax.set_xlabel("Feature-0")

ax.set_ylabel("Feature-1")

ax.set_zlabel("Feature-1**2")

plt.show()

元の特徴量に対する決定境界

上の例では特徴量は3つだが、最後の特徴量は2つ目の特徴量f₁から計算される量であり、実質は2つの特徴量が決まれば決定境界が決まる。3次元空間内の平面の決定境界を以下のように書きなおすと、このことが確認できる。

(3) $\begin{align*} f_0 &= \frac{-b - w_1 f_1 - w_2 {f_1}^2}{w_0} \\ &= -9.02 +1.69 f_1 + 0.46 {f_1}^2 \\ &= 0.46(f_1 +1.83)^2 - 10.56 \end{align*}$

これを2つの特徴量に対する決定境界として描画したのが以下の図で、境界が2次関数となっているのが確認できる。

SVMそのものは線形の決定境界しか得られないが、非線形化した特徴量を追加することによって、より複雑な決定境界とすることができる。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import make_blobs
from sklearn.svm import LinearSVC
from mpl_toolkits.mplot3d import Axes3D

X, y = make_blobs(centers=4, random_state=8)
y = y % 2
X_new = np.hstack((X, X[:, 1].reshape(-1, 1)**2))

linsvc = LinearSVC(max_iter=3700).fit(X_new, y)
intercept = linsvc.intercept_

x_min, x_max = -12, 12
y_min, y_max = -15, 15
u = np.linspace(x_min, x_max, 500)
v = np.linspace(y_min, y_max, 500)
u, v = np.meshgrid(u, v)
w = v**2
decision = linsvc.predict(np.c_[u.ravel(), v.ravel(), w.ravel()])

fig, ax = plt.subplots()

ax.scatter(X_new[y==0, 0], X_new[y==0, 1], marker='o', s=60, ec='k')
ax.scatter(X_new[y==1, 0], X_new[y==1, 1], marker='^', s=60, ec='k')
ax.contourf(u, v, decision.reshape(u.shape),
    levels=1, colors=['tab:blue', 'tab:orange'], alpha=0.4)
ax.contour(u, v, decision.reshape(u.shape), levels=1, colors='k')

ax.set_xlim(x_min, x_max)
ax.set_ylim(y_min, y_max)
ax.set_xlabel("Feature-0")
ax.set_ylabel("Feature-1")

plt.show()

import numpy as np

import matplotlib.pyplot as plt

from sklearn.datasets import make_blobs

from sklearn.svm import LinearSVC

from mpl_toolkits.mplot3d import Axes3D

X, y = make_blobs(centers=4, random_state=8)

y = y % 2

X_new = np.hstack((X, X[:, 1].reshape(-1, 1)**2))

linsvc = LinearSVC(max_iter=3700).fit(X_new, y)

intercept = linsvc.intercept_

x_min, x_max = -12, 12

y_min, y_max = -15, 15

u = np.linspace(x_min, x_max, 500)

v = np.linspace(y_min, y_max, 500)

u, v = np.meshgrid(u, v)

w = v**2

decision = linsvc.predict(np.c_[u.ravel(), v.ravel(), w.ravel()])

fig, ax = plt.subplots()

ax.scatter(X_new[y==0, 0], X_new[y==0, 1], marker='o', s=60, ec='k')

ax.scatter(X_new[y==1, 0], X_new[y==1, 1], marker='^', s=60, ec='k')

ax.contourf(u, v, decision.reshape(u.shape),

levels=1, colors=['tab:blue', 'tab:orange'], alpha=0.4)

ax.contour(u, v, decision.reshape(u.shape), levels=1, colors='k')

ax.set_xlim(x_min, x_max)

ax.set_ylim(y_min, y_max)

ax.set_xlabel("Feature-0")

ax.set_ylabel("Feature-1")

plt.show()

カーネルトリック

概要

上記の例では特徴量の1つを2次として新たな特徴量とした。特徴量の非線形化としては、このように特徴量の累乗とするほか、異なる特徴量同士の積を交互作用として導入することが考えられる。ただし、特徴量の数が多くなった時に、それらの全ての組み合わせに対する積を考えると、計算量が膨れ上がる。カーネルトリック(kernel trick)とは、拡張された特徴量空間でのデータ間の距離を、実際の拡張計算をせずに行う方法らしい。

受け売りをそのまま書いておくと、SVMで広く用いられているカーネルトリックのマッピング方法は以下の2つとのこと。

多項式カーネル(polynomial kernel)：もとの特徴量の特定の次数までの全ての多項式を計算
放射既定関数(radial basis function: RBF)カーネルとも呼ばれるガウシアンカーネル：直感的には全次数の全ての多項式を考えるが、次数が高くなるにつれて特徴量の重要性を小さくする

以下はforgeデータセットに対して、カーネルトリックを用いたSVCを適用した例。直線はLinearSVCによる決定境界で、曲線はガウシアンカーネル(RBF)によるSVCの決定境界で、カーネル関数は以下のような形。

(4) $\begin{equation*} k_{\rm rbf}(x_1, x_2) = \exp \left( -\gamma || x_1 - x_2 ||^2 \right) \end{equation*}$

scikit-learnのSVCの引数で、kernel='rbf'、C=10、gamma=0.1と指定している。

線形モデルの決定境界が直線なのに対して、カーネルトリックによる決定境界は、非線形化した特徴量を導入していることから曲線となっている。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.svm import SVC, LinearSVC

X = np.array( \
[[ 9.96346605,  4.59676542],
 [11.0329545,  -0.16816717],
 [11.54155807,  5.21116083],
 .....
 [ 9.50169345,  1.93824624],
 [ 9.15072323,  5.49832246],
 [11.563957,    1.3389402 ]]
)
y = np.array( \
[1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0]
)

f0min, f0max = 7.5, 12.5
f1min, f1max = -1, 6

linsvc = LinearSVC(max_iter=4000).fit(X, y)
intercept = linsvc.intercept_
coef = linsvc.coef_.ravel()

svc = SVC(kernel='rbf', C=10, gamma=0.1).fit(X, y)
u = np.linspace(f0min, f0max, 400)
v = np.linspace(f1min, f1max, 400)
u, v = np.meshgrid(u, v)
pred = svc.predict(np.c_[u.ravel(), v.ravel()]).reshape(u.shape)

fig, ax = plt.subplots()

ax.scatter(X[y==0][:, 0], X[y==0][:, 1],
    marker='o', s=60, fc='tab:blue', ec='k')
ax.scatter(X[y==1][:, 0], X[y==1][:, 1],
    marker='^', s=60, fc='tab:orange', ec='k')

sv_class = y[svc.support_]
ax.scatter(svc.support_vectors_[sv_class==0][:, 0],
           svc.support_vectors_[sv_class==0][:, 1],
           marker='o', s=150, fc='tab:blue', ec='blue', linewidth=3)
ax.scatter(svc.support_vectors_[sv_class==1][:, 0],
           svc.support_vectors_[sv_class==1][:, 1],
           marker='^', s=150, fc='tab:orange', ec='red', linewidth=3)

f1 = lambda f0: -(intercept + coef[0]*f0) / coef[1]
ax.plot([f0min, f0max], [f1(f0min), f1(f0max)])

ax.contour(u, v, pred, levels=[0.5])

ax.set_xlim(f0min, f0max)
ax.set_ylim(f1min, f1max)
ax.tick_params(bottom=False, left=False, labelbottom=False, labelleft=False)
ax.set_xlabel("Feature-0")
ax.set_ylabel("Feature-1")

plt.show()

import numpy as np

import matplotlib.pyplot as plt

from sklearn.svm import SVC, LinearSVC

X = np.array( \

[[ 9.96346605, 4.59676542],

[11.0329545, -0.16816717],

[11.54155807, 5.21116083],

.....

[ 9.50169345, 1.93824624],

[ 9.15072323, 5.49832246],

[11.563957, 1.3389402 ]]

)

y = np.array( \

[1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0]

)

f0min, f0max = 7.5, 12.5

f1min, f1max = -1, 6

linsvc = LinearSVC(max_iter=4000).fit(X, y)

intercept = linsvc.intercept_

coef = linsvc.coef_.ravel()

svc = SVC(kernel='rbf', C=10, gamma=0.1).fit(X, y)

u = np.linspace(f0min, f0max, 400)

v = np.linspace(f1min, f1max, 400)

u, v = np.meshgrid(u, v)

pred = svc.predict(np.c_[u.ravel(), v.ravel()]).reshape(u.shape)

fig, ax = plt.subplots()

ax.scatter(X[y==0][:, 0], X[y==0][:, 1],

marker='o', s=60, fc='tab:blue', ec='k')

ax.scatter(X[y==1][:, 0], X[y==1][:, 1],

marker='^', s=60, fc='tab:orange', ec='k')

sv_class = y[svc.support_]

ax.scatter(svc.support_vectors_[sv_class==0][:, 0],

svc.support_vectors_[sv_class==0][:, 1],

marker='o', s=150, fc='tab:blue', ec='blue', linewidth=3)

ax.scatter(svc.support_vectors_[sv_class==1][:, 0],

svc.support_vectors_[sv_class==1][:, 1],

marker='^', s=150, fc='tab:orange', ec='red', linewidth=3)

f1 = lambda f0: -(intercept + coef[0]*f0) / coef[1]

ax.plot([f0min, f0max], [f1(f0min), f1(f0max)])

ax.contour(u, v, pred, levels=[0.5])

ax.set_xlim(f0min, f0max)

ax.set_ylim(f1min, f1max)

ax.tick_params(bottom=False, left=False, labelbottom=False, labelleft=False)

ax.set_xlabel("Feature-0")

ax.set_ylabel("Feature-1")

plt.show()

scikit-learnのSVCクラスには、サポートベクターに関する以下のパラメーターがある。

support_：データセットにおけるサポートベクターのインデックス（1次元配列）
support_vector_：サポートベクターの配列（2次元配列）

38～44行目で、これらのパラメーターを使ってサポートベクターを強調表示している。

パラメータ調整

SVCモデルでパラメーターCとgammaの値を変化させたときの決定境界は以下の通り。

gammaはガウシアンカーネルの直径（σ2に相当）の逆数で、この値が小さいと直径が大きくなり、より多くの点を近いと判断するようになる。左の方はgammaが小さく広域のデータをまとめようとするため、決定境界は大まかとなり、右の方はgammaが大きく近いもの同士をまとめようとする傾向となる。

Cは正則化の強さの逆数で、上の方ほどCの値が小さく正則化が強く効くため、決定境界はよりまっすぐとなり、下の方ほど正則化が弱く個々のデータの影響を受ける。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.svm import SVC
from itertools import product

X = np.array( \
[[ 9.96346605,  4.59676542],
 [11.0329545,  -0.16816717],
 [11.54155807,  5.21116083],
 .....
 [ 9.50169345,  1.93824624],
 [ 9.15072323,  5.49832246],
 [11.563957,    1.3389402 ]]
)
y = np.array( \
[1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0]
)

f0min, f0max = 7.5, 12.5
f1min, f1max = -1, 6
u = np.linspace(f0min, f0max, 400)
v = np.linspace(f1min, f1max, 400)
u, v = np.meshgrid(u, v)

C_list = [0.1, 1, 1000]
gamma_list = [0.1, 1, 10]
params = product(C_list, gamma_list)

plt.rcParams['font.size'] = 6
fig, axs = plt.subplots(3, 3, figsize=(6.4, 4.8))
axs_1d = axs.ravel()
fig.subplots_adjust(hspace=0.3)

for ax, param in zip(axs_1d, params):
    svc = SVC(kernel='rbf', C=param[0], gamma=param[1]).fit(X, y)
    pred = svc.predict(
            np.hstack([u.ravel().reshape(-1, 1), v.ravel().reshape(-1, 1)])
        ).reshape(u.shape)
    ax.scatter(X[y==0][:, 0], X[y==0][:, 1], marker='o')
    ax.scatter(X[y==1][:, 0], X[y==1][:, 1], marker='^')
    ax.contour(u, v, pred, levels=[0.5])
    ax.tick_params(left=False, bottom=False, labelleft=False, labelbottom=False)
    ax.set_title("C={:.1f}, gamma={:.1f}".format(param[0], param[1]))

plt.show()

import numpy as np

import matplotlib.pyplot as plt

from sklearn.svm import SVC

from itertools import product

X = np.array( \

[[ 9.96346605, 4.59676542],

[11.0329545, -0.16816717],

[11.54155807, 5.21116083],

.....

[ 9.50169345, 1.93824624],

[ 9.15072323, 5.49832246],

[11.563957, 1.3389402 ]]

)

y = np.array( \

[1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0]

)

f0min, f0max = 7.5, 12.5

f1min, f1max = -1, 6

u = np.linspace(f0min, f0max, 400)

v = np.linspace(f1min, f1max, 400)

u, v = np.meshgrid(u, v)

C_list = [0.1, 1, 1000]

gamma_list = [0.1, 1, 10]

params = product(C_list, gamma_list)

plt.rcParams['font.size'] = 6

fig, axs = plt.subplots(3, 3, figsize=(6.4, 4.8))

axs_1d = axs.ravel()

fig.subplots_adjust(hspace=0.3)

for ax, param in zip(axs_1d, params):

svc = SVC(kernel='rbf', C=param[0], gamma=param[1]).fit(X, y)

pred = svc.predict(

np.hstack([u.ravel().reshape(-1, 1), v.ravel().reshape(-1, 1)])

).reshape(u.shape)

ax.scatter(X[y==0][:, 0], X[y==0][:, 1], marker='o')

ax.scatter(X[y==1][:, 0], X[y==1][:, 1], marker='^')

ax.contour(u, v, pred, levels=[0.5])

ax.tick_params(left=False, bottom=False, labelleft=False, labelbottom=False)

ax.set_title("C={:.1f}, gamma={:.1f}".format(param[0], param[1]))

plt.show()

Breast cancerデータへの適用例

Breast cancerデータへの適用例で、特徴量データの大きさやレンジが大きくばらついていること、特徴量データをそのまま使った場合に過学習となること、特徴量データに前処理を施して正規化(normalize)した場合に精度が向上することを示している。

SVMの特徴

SVMの特徴量を受け売りのまままとめておく。

データにわずかな特徴量しかない場合も複雑な決定境界を生成可能（低次元でも高次元でもうまく機能）
サンプルの個数が大きくなるとうまく機能しない（10万サンプルくらいになると、実行時間やメモリ使用量の面で難しくなる
注意深いデータの前処理とパラメーター調整が必要
検証が難しい（予測に対する理由を理解することが難しい）
RBFの場合、Cやgammaを大きくするとより複雑なモデルになる（2つのパラメーターは強く相関するため、同時に調整する必要がある）

今後の課題～覚え書き

カーネル関数

(5) $\begin{equation*} K(\boldsymbol{x}_1, \boldsymbol{x}_2) = \sum \phi(\boldsymbol{x}_1) \phi(\boldsymbol{x}_2) \end{equation*}$

多項式カーネル

(6) $\begin{equation*} K(\boldsymbol{x}_1, \boldsymbol{x}_2) = (\boldsymbol{x}_1 \cdot \boldsymbol{x}_2 + 1 )^d \end{equation*}$

ガウシアンカーネル

(7) $\begin{equation*} K(\boldsymbol{x}_1, \boldsymbol{x}_2) = \exp \left(- \frac{||(\boldsymbol{x}_1 - \boldsymbol{x}_2 ||^2}{2\sigma^2} \right) \end{equation*}$

SVMの定式化

2020-06-24 / tau / コメントする

SVMの定式化

SVMのクラス分類の条件

2つの特徴量を持つデータが2つのクラスに分かれているとする。ここで下図のように、1つの直線によって、2つのクラスを完全に分離できるとする。

このとき、直線lによって分割したとして、以下の符号によってクラスを分離する。

(1) $\begin{equation*} \left\{ \begin{align} a x_1 + b x_2 + c > 0 &\rightarrow \rm{Class1} \\ a x_1 + b x_2 + c < 0 & \rightarrow \rm{Class2} \end{align} \right. \end{equation*}$

ここでラベル変数t_iを導入する。t_iはデータiがClass1/2のいずれに属するかを示す変数で、Class1ならt_i > 0、Class2ならt_i < 0と定義する。

(2) $\begin{equation*} \left\{ \begin{array}{lll} t_i = 1 & x_i \in \rm{Class1} & (a x_{i1} + b x_{i2} + c > 0) \\ t_i = -1 & x_i \in \rm{Class2} & (a x_{i1} + b x_{i2} + c < 0) \\ \end{array} \right. \end{equation*}$

このラベル変数を用いて、クラスの条件式は以下のように統一される。

(3) $\begin{equation*} t_i (a x_{i1} + b x_{i2} + c) > 0 \end{equation*}$

SVMにおいては、すべてのデータについてこの式が満足されるようにa, b, cを決定する。これらはすべてa, b, cに対する制約条件だが、どのようにこれらの値を求めるべきか、その目的関数が必要になる。SVMでは、これをマージン最大化により行う。

マージン最大化

ある直線l₁によって、下図のようにデータセットがClass1/2に分類できるとする。このときl₁に対してClass1/2の最も直線に近いデータを”サポートベクター”と呼ぶ。また、これらのサポートベクターに対応するl₁と平行な直線間の距離を”マージン”と呼ぶ。

ところで、l₁とは異なる別の直線l₂を選ぶと、異なるサポートベクターに対してより大きなマージンを得ることができる。SVMでは、式(3)のもとでこのマージンを最大化するような直線lを探すこととなる。

直線lに対するサポートベクターの対を(x⁺, x⁻)とすると、それぞれからlへの距離dは以下のように表現される。

(4) $\begin{equation*} d = \frac{|a x^+_1 + b x^+_2 + c|}{a^2 + b^2} = \frac{|a x^-_1 + b x^-_2 + c|}{a^2 + b^2} \end{equation*}$

ここで直線lはマージンの端にある平行な2つの直線の中央にあることから、上式の分子は同じ値となる。この値でdを除したものを改めて $\tilde{d}$ と置くと、dの最大化問題は $1/\tilde{d}=\sqrt{a^2+b^2}$ の最小化問題となる。これに式(3)の制約条件を加味して、問題は以下の制約条件付き最小化問題となる。

(5) $\begin{align*} \min a^2 + b^2 \quad {\rm s.t.} \; t_i (a x_{i1} + b x_{i2} + c) > 0 \; (i=1~n) \end{align*}$

今後の課題

ここから先の定式化
ソフトマージンの導出

勾配ブースティング

2020-06-21 / tau / コメントする

概要

勾配ブースティング(gradient boosthing)は、ランダムフォレストと同じく複数の決定木を組み合わせてモデルを強化する手法。ランダムフォレストと異なる点は、最初から複数の決定木を使うのではなく、1つずつ順番に決定木を増やしていく。その際に追加される決定木はそれぞれ深さ1～5くらいの浅い木（弱学習機：weak learner）で、直前の適合不足を補うように学習する。

勾配ブースティングの主なパラメーターは弱学習機の数(n_estimators)と学習率(learning_rate)で、学習率を大きくすると個々の弱学習機の補正を強化しモデルは複雑になる。

cancerデータへの適用

Pythonのscikit-learnにあるGradienBoostingClassifierをbreast_cancerデータに適用する例。”Pythonではじめる機械学習”の”2.3.6.2 勾配ブースティング回帰木”掲載のコードに沿って確認するが、バージョンの違いのためか、結果が異なる。いくつかのデフォルトのパラメーターを明示的に設定／変更してみたが、書籍に掲載されている結果には至っていない。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_breast_cancer
from sklearn.model_selection import train_test_split
from sklearn.ensemble import GradientBoostingClassifier

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

gbcf = GradientBoostingClassifier(random_state=1)
gbcf.fit(X_train, y_train)
print("Training score: {:.3f}".format(gbcf.score(X_train, y_train)))
print("Test score    : {:.3f}".format(gbcf.score(X_test, y_test)))

gbcf = GradientBoostingClassifier(max_depth=1, random_state=0)
gbcf.fit(X_train, y_train)
print("Training score: {:.3f}".format(gbcf.score(X_train, y_train)))
print("Test score    : {:.3f}".format(gbcf.score(X_test, y_test)))

fig, ax = plt.subplots(figsize=(8, 4.8))
fig.subplots_adjust(left=0.3)
ax.barh(ds.feature_names, gbcf.feature_importances_)
ax.set_xlabel("feature importance")
plt.show()

gbcf = GradientBoostingClassifier(learning_rate=0.01, random_state=0)
gbcf.fit(X_train, y_train)
print("Training score: {:.3f}".format(gbcf.score(X_train, y_train)))
print("Test score    : {:.3f}".format(gbcf.score(X_test, y_test)))

import numpy as np

import matplotlib.pyplot as plt

from sklearn.datasets import load_breast_cancer

from sklearn.model_selection import train_test_split

from sklearn.ensemble import GradientBoostingClassifier

ds = load_breast_cancer()

X_train, X_test, y_train, y_test =\

train_test_split(ds.data, ds.target, random_state=0)

gbcf = GradientBoostingClassifier(random_state=1)

gbcf.fit(X_train, y_train)

print("Training score: {:.3f}".format(gbcf.score(X_train, y_train)))

print("Test score : {:.3f}".format(gbcf.score(X_test, y_test)))

gbcf = GradientBoostingClassifier(max_depth=1, random_state=0)

gbcf.fit(X_train, y_train)

print("Training score: {:.3f}".format(gbcf.score(X_train, y_train)))

print("Test score : {:.3f}".format(gbcf.score(X_test, y_test)))

fig, ax = plt.subplots(figsize=(8, 4.8))

fig.subplots_adjust(left=0.3)

ax.barh(ds.feature_names, gbcf.feature_importances_)

ax.set_xlabel("feature importance")

plt.show()

gbcf = GradientBoostingClassifier(learning_rate=0.01, random_state=0)

gbcf.fit(X_train, y_train)

print("Training score: {:.3f}".format(gbcf.score(X_train, y_train)))

print("Test score : {:.3f}".format(gbcf.score(X_test, y_test)))

最初に試したのが以下のコード。ここでテストスコアが書籍にある0.958にならない。min_samples_split=5とすると書籍と同じ結果になるが、以降の特徴量重要度やlearning_rateの変更結果は再現されない。

gbcf = GradientBoostingClassifier(random_state=0)
gbcf.fit(X_train, y_train)
print("Training score: {:.3f}".format(gbcf.score(X_train, y_train)))
print("Test score    : {:.3f}".format(gbcf.score(X_test, y_test)))

# Training score: 1.000
# Test score    : 0.965

gbcf = GradientBoostingClassifier(random_state=0)

gbcf.fit(X_train, y_train)

print("Training score: {:.3f}".format(gbcf.score(X_train, y_train)))

print("Test score : {:.3f}".format(gbcf.score(X_test, y_test)))

# Training score: 1.000

# Test score : 0.965

過剰適合に対してmax_depth=1と強力な枝刈りをした場合。この結果は小数点以下3桁の表示で書籍と一致している。

gbcf = GradientBoostingClassifier(max_depth=1, random_state=0)
gbcf.fit(X_train, y_train)
print("Training score: {:.3f}".format(gbcf.score(X_train, y_train)))
print("Test score    : {:.3f}".format(gbcf.score(X_test, y_test)))

# Training score: 0.991
# Test score    : 0.972

gbcf = GradientBoostingClassifier(max_depth=1, random_state=0)

gbcf.fit(X_train, y_train)

print("Training score: {:.3f}".format(gbcf.score(X_train, y_train)))

print("Test score : {:.3f}".format(gbcf.score(X_test, y_test)))

# Training score: 0.991

# Test score : 0.972

learning_rateをデフォルトの0.1から0.01に変更した場合の結果も書籍と一致する。今回の再現結果では、デフォルト状態からテストスコアは改善されていない。

gbcf = GradientBoostingClassifier(learning_rate=0.01, random_state=0)
gbcf.fit(X_train, y_train)
print("Training score: {:.3f}".format(gbcf.score(X_train, y_train)))
print("Test score    : {:.3f}".format(gbcf.score(X_test, y_test)))

# Training score: 0.988
# Test score    : 0.965

gbcf = GradientBoostingClassifier(learning_rate=0.01, random_state=0)

gbcf.fit(X_train, y_train)

print("Training score: {:.3f}".format(gbcf.score(X_train, y_train)))

print("Test score : {:.3f}".format(gbcf.score(X_test, y_test)))

# Training score: 0.988

# Test score : 0.965

なお、事前剪定を強化したケースのグラフが、書籍と大きく異なる。横軸の値が倍ほどになっており、worst concave points、worst perimeter、mean concave pointsが重要度の大半を占めている。書籍では他の多くの特徴量も重要度がある程度高い点と異なっている。

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

fig, ax = plt.subplots(figsize=(8, 4.8))
fig.subplots_adjust(left=0.3)
ax.barh(ds.feature_names, gbcf.feature_importances_)
ax.set_xlabel("feature importance")
plt.show()

gbcf = GradientBoostingClassifier(max_depth=1, random_state=0)

gbcf.fit(X_train, y_train)

fig, ax = plt.subplots(figsize=(8, 4.8))

fig.subplots_adjust(left=0.3)

ax.barh(ds.feature_names, gbcf.feature_importances_)

ax.set_xlabel("feature importance")

plt.show()

今後確認したい点

勾配ブースティングの基本的な考え方の整理
簡単な事例での勾配ブースティングの挙動確認
回帰への適用
異なるモデルの組み合わせ

ランダムフォレストの概要

2020-06-20 / tau / コメントする

概要

“Pythonではじめる機械学習”のランダムフォレストの写経。

ランダムフォレストは決定木のアンサンブル法の1つ。異なる複数の決定木をランダムに発生させて平均をとることで、個々の決定木の過剰適合を打ち消すという考え方。

実行例

以下は、scikit-learnのランダムフォレストのモデルRandomForestClassifierでmoonsデータセットをクラス分類した例で、ランダムに生成された5つの木とそれらを平均したランダムフォレストの決定領域を示している。

以下は実装コードで、draw_decision_field()とdraw_tree_boundary()はコードを省略。

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patch
from sklearn.datasets import make_moons
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestClassifier


def draw_decision_field(clf, ax, x0s, x1s, n_areas=2,
        colors=['tab:blue', 'tab:orange'], alpha=0.5, fill=True):
# コード省略

def draw_tree_boundary(tree, ax, left, right, bottom, top,
        i_node=0, stop_level=None, n_level=0):
# コード省略 

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)

x0_min, x0_max = -1.8, 2.5
x1_min, x1_max = -1.0, 1.8
x0s = np.linspace(x0_min, x0_max, 50)
x1s = np.linspace(x1_min, x1_max, 50)

forest = RandomForestClassifier(n_estimators=5, random_state=2)
forest.fit(X_train, y_train)

fig, axs = plt.subplots(2, 3, figsize=(8, 4.8))
fig.subplots_adjust(hspace=0.3)
axs_1d = axs.reshape(-1)

for i, ax in enumerate(axs_1d[:-1]):
    draw_tree_boundary(tree=forest.estimators_[i].tree_, ax=ax,
        left=x0_min, right=x0_max, bottom=x1_min, top=x1_max)
    ax.set_title("Tree-{}".format(i))

draw_decision_field(forest, axs_1d[-1], x0s, x1s)
axs_1d[-1].set_title("Random Forest")

for ax in axs_1d:
    ax.scatter(X_train[y_train==0][:, 0], X_train[y_train==0][:, 1],
        marker='o', s=15, fc='tab:blue', ec='k')
    ax.scatter(X_train[y_train==1][:, 0], X_train[y_train==1][:, 1],
        marker='^', s=15, fc='tab:orange', ec='k')

    ax.set_xlim(x0_min, x0_max)
    ax.set_ylim(x1_min, x1_max)
    ax.tick_params(bottom=False, left=False, labelbottom=False, labelleft=False)

plt.show()

import numpy as np

import matplotlib.pyplot as plt

import matplotlib.patches as patch

from sklearn.datasets import make_moons

from sklearn.model_selection import train_test_split

from sklearn.ensemble import RandomForestClassifier

def draw_decision_field(clf, ax, x0s, x1s, n_areas=2,

colors=['tab:blue', 'tab:orange'], alpha=0.5, fill=True):

# コード省略

def draw_tree_boundary(tree, ax, left, right, bottom, top,

i_node=0, stop_level=None, n_level=0):

# コード省略

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)

x0_min, x0_max = -1.8, 2.5

x1_min, x1_max = -1.0, 1.8

x0s = np.linspace(x0_min, x0_max, 50)

x1s = np.linspace(x1_min, x1_max, 50)

forest = RandomForestClassifier(n_estimators=5, random_state=2)

forest.fit(X_train, y_train)

fig, axs = plt.subplots(2, 3, figsize=(8, 4.8))

fig.subplots_adjust(hspace=0.3)

axs_1d = axs.reshape(-1)

for i, ax in enumerate(axs_1d[:-1]):

draw_tree_boundary(tree=forest.estimators_[i].tree_, ax=ax,

left=x0_min, right=x0_max, bottom=x1_min, top=x1_max)

ax.set_title("Tree-{}".format(i))

draw_decision_field(forest, axs_1d[-1], x0s, x1s)

axs_1d[-1].set_title("Random Forest")

for ax in axs_1d:

ax.scatter(X_train[y_train==0][:, 0], X_train[y_train==0][:, 1],

marker='o', s=15, fc='tab:blue', ec='k')

ax.scatter(X_train[y_train==1][:, 0], X_train[y_train==1][:, 1],

marker='^', s=15, fc='tab:orange', ec='k')

ax.set_xlim(x0_min, x0_max)

ax.set_ylim(x1_min, x1_max)

ax.tick_params(bottom=False, left=False, labelbottom=False, labelleft=False)

plt.show()

上の例でn_estimatorsの数を1から9まで増やしていったときの様子を示す。徐々に細かい枝が取り払われて、境界が滑らかになっていく様子がわかる。

ランダムフォレストの考え方

概要

ランダムフォレストの大まかな手順は以下の通り。

訓練／ランダムフォレストの構築
1. 決定木の数を指定する
2. 決定木の数の分だけ、訓練データからランダムにデータを生成する
3. 各決定木を学習させる
予測段階
1. 予測したいデータの特徴量を与える
2. その特徴量に対して各決定木がクラス分類
3. 各決定木の分類結果の多数決でクラスを決定

ランダムフォレストの構築

決定木の数

決定木の数をRandomForestClassifierのn_estimatorsパラメーターで、構築する決定木の数を指定する。

決定木に与えるデータの生成

各決定木に与えるデータセットをランダムに設定。特徴量の一部または全部を選び、各特徴量のデータをブートストラップサンプリングによってランダムに選び出す。

RandomForestClassifierでは、選ぶ特徴量の数をmax_featuresで指定し、各特徴量のデータをブートストラップサンプリングで選ぶかどうかをbootstrap(=True/False)で設定する。これらの乱数系列をrandom_stateパラメーターで指定する。

冒頭の事例の場合、特徴量は2つなので常にいずれも選び、bootstrapはデフォルトのTrue、決定木の数を3つ、乱数系列を3と指定している。

forest = RandomForestClassifier(n_estimators=3, random_state=2)
forest.fit(X_train, y_train)

1 2	forest = RandomForestClassifier(n_estimators=3, random_state=2) forest.fit(X_train, y_train)

ランダムフォレストによる予測

ある特徴量セットを持つデータのクラスを予測する。与えられたデータに対して、各決定木がクラスを判定し、その結果の多数決でそのデータのクラスを決定する。

クラス分類でランダムフォレストが多数決で任意の点のクラスを決定する様子を確認する。以下は10個のmoonsデータセットに対して3つの分類木を適用した例。

訓練後のランダムフォレストにおいて、各点のクラスが多数決で決められていることが、以下の図で確認できる。

cancerデータによる確認

精度

breast_cancerデータに対してランダムフォレストを適用する。100個の決定木を準備し、他のパラメーターはデフォルトのままで実行する。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_breast_cancer
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestClassifier

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

print("Training score: {:.3f}".format(forest.score(X_train, y_train)))
print("Test score    : {:.3f}".format(forest.score(X_test, y_test)))

fig, ax = plt.subplots()
fig.subplots_adjust(left=0.3)
ax.barh(ds.feature_names, forest.feature_importances_)
plt.show()

import numpy as np

import matplotlib.pyplot as plt

from sklearn.datasets import load_breast_cancer

from sklearn.model_selection import train_test_split

from sklearn.ensemble import RandomForestClassifier

ds = load_breast_cancer()

X_train, X_test, y_train, y_test = \

train_test_split(ds.data, ds.target, random_state=0)

forest = RandomForestClassifier(n_estimators=100, random_state=0)

forest.fit(X_train, y_train)

print("Training score: {:.3f}".format(forest.score(X_train, y_train)))

print("Test score : {:.3f}".format(forest.score(X_test, y_test)))

fig, ax = plt.subplots()

fig.subplots_adjust(left=0.3)

ax.barh(ds.feature_names, forest.feature_importances_)

plt.show()

このコードの前半の実行結果は以下の通りで、訓練データに対して完全適合し、テストデータに対しても97.2%の精度を示している。

Training score: 1.000
Test score    : 0.972

1 2	Training score: 1.000 Test score : 0.972

【注】決定木の数による違い

n_estimatorの数によって、テストスコアが違ってくるが、このケースでは100以上決定木の数を多くしてもテストスコアは向上しない。

n_estimators = 10 → 0.951
n_estimators = 50 → 0.965
n_estimators ≥ 100 → 0.972

特徴量重要度

後半の実行結果は以下のように表示される。単独の決定木による特徴量重要度と比べて全ての特徴量が0以上の重要度となっている。決定木の時に最重要であったworst radiusも重要度が高いが、worst perimeterが最も重要度が高く、worst concave pointsやconcave pointsも重要度が高い。

今後の課題

ランダムフォレストにおける特徴量重要度の意義
n_jobsの効果
RandomForestClassifierのrandom_stateによる違い
max_features、max_depth、min_samples_leafなどの影響
- デフォルトはmax_features=sqrt(n_features)、max_depth=None(all)、min_samples_leaf=1

決定境界／クラス分類の分布を描く関数

2020-06-11 / tau / コメントする

概要

2つの特徴量を持つデータセットを学習したモデルに対し、2次元の特徴量空間における決定境界やクラス分類の分布を描く関数の例。

draw_decision_boundary()で決定境界の線を描き、draw_decision_area()で領域のクラス分布を色分けで表示する。

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patch
from sklearn.datasets import make_moons
from sklearn.neighbors import KNeighborsClassifier


def draw_decision_boundary(clf, ax, x0s, x1s,
        threshold=0, color='k', alpha=1.0):

    y_predicted = np.empty((len(x1s), len(x0s)))

    for row, x1 in enumerate(x1s):
        for col, x0 in enumerate(x0s):
            y_predicted[row, col] = clf.predict(np.array([[x0, x1]]))

    ax.contour(x0s, x1s, y_predicted,
        colors=color, levels=[threshold] , alpha=alpha)


def draw_decision_field(clf, ax, x0s, x1s, n_areas=2,
        colors=['tab:blue', 'tab:orange'], alpha=0.5, fill=True):

    y_predicted = np.empty((len(x1s), len(x0s)))

    for row, x1 in enumerate(x1s):
        for col, x0 in enumerate(x0s):
            y_predicted[row, col] = clf.predict(np.array([[x0, x1]]))
    ax.contourf(x0s, x1s, y_predicted,
        colors=colors, levels=n_areas-1, alpha=alpha)


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

x0_min, x0_max = -2.0, 2.5
x1_min, x1_max = -1.0, 1.5
x0s = np.linspace(x0_min, x0_max, 50)
x1s = np.linspace(x1_min, x1_max, 50)

fig, axs = plt.subplots(1, 2, figsize=(11, 4.8))
fig.subplots_adjust(wspace=0.3)
axs_1d = axs.reshape(-1)

knn = KNeighborsClassifier(n_neighbors=3)
knn.fit(X, y)

for ax in axs_1d:
    ax.scatter(X[y==0][:, 0], X[y==0][:, 1], marker='o')
    ax.scatter(X[y==1][:, 0], X[y==1][:, 1], marker='^')

    ax.set_xlim(x0_min, x0_max)
    ax.set_ylim(x1_min, x1_max)
    ax.set_xlabel("Feature-0")
    ax.set_ylabel("Feature-1")

draw_decision_boundary(knn, axs[0], x0s, x1s, threshold=0.5)
draw_decision_field(knn, axs[1], x0s, x1s, alpha=0.2)

axs[0].set_title("decision boundary")
axs[1].set_title("decision field")

plt.show()

import numpy as np

import matplotlib.pyplot as plt

import matplotlib.patches as patch

from sklearn.datasets import make_moons

from sklearn.neighbors import KNeighborsClassifier

def draw_decision_boundary(clf, ax, x0s, x1s,

threshold=0, color='k', alpha=1.0):

y_predicted = np.empty((len(x1s), len(x0s)))

for row, x1 in enumerate(x1s):

for col, x0 in enumerate(x0s):

y_predicted[row, col] = clf.predict(np.array([[x0, x1]]))

ax.contour(x0s, x1s, y_predicted,

colors=color, levels=[threshold] , alpha=alpha)

def draw_decision_field(clf, ax, x0s, x1s, n_areas=2,

colors=['tab:blue', 'tab:orange'], alpha=0.5, fill=True):

y_predicted = np.empty((len(x1s), len(x0s)))

for row, x1 in enumerate(x1s):

for col, x0 in enumerate(x0s):

y_predicted[row, col] = clf.predict(np.array([[x0, x1]]))

ax.contourf(x0s, x1s, y_predicted,

colors=colors, levels=n_areas-1, alpha=alpha)

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

x0_min, x0_max = -2.0, 2.5

x1_min, x1_max = -1.0, 1.5

x0s = np.linspace(x0_min, x0_max, 50)

x1s = np.linspace(x1_min, x1_max, 50)

fig, axs = plt.subplots(1, 2, figsize=(11, 4.8))

fig.subplots_adjust(wspace=0.3)

axs_1d = axs.reshape(-1)

knn = KNeighborsClassifier(n_neighbors=3)

knn.fit(X, y)

for ax in axs_1d:

ax.scatter(X[y==0][:, 0], X[y==0][:, 1], marker='o')

ax.scatter(X[y==1][:, 0], X[y==1][:, 1], marker='^')

ax.set_xlim(x0_min, x0_max)

ax.set_ylim(x1_min, x1_max)

ax.set_xlabel("Feature-0")

ax.set_ylabel("Feature-1")

draw_decision_boundary(knn, axs[0], x0s, x1s, threshold=0.5)

draw_decision_field(knn, axs[1], x0s, x1s, alpha=0.2)

axs[0].set_title("decision boundary")

axs[1].set_title("decision field")

plt.show()

関数の使い方

それぞれの関数単体では特にパッケージは必要ないが、いくつかのパラメーターは一定のクラスを想定している。

`draw_decision_boundary()`

draw_decision_boundary(clf, ax, x0s, x1s, threshold, color, alpha)

clf: 学習済みのクラス分類モデルのインスタンスを指定する。predict()メソッドを持つこと（引数は2次元配列を想定）。
ax: 決定境界を描くAxesオブジェクト。
x0s, x1s: クラスを計算する領域の計算点の座標を1次元配列で指定。
threshold: 決定境界の値を整数で与える。デフォルトは0。決定値がクラスラベル（例えば0と1）で与えられる場合はその平均（たとえば0.5）を与える。
color: 決定境界の場合はカラーコード。デフォルトは’k’（黒）
alpha: 分布図の場合の塗りつぶしの透明度を実数で指定。デフォルトは1（不透明）

`draw_decision_field()`

draw_decision_field(clf, ax, x0s, x1s, n_areas=2, colors, alpha)

clf: 学習済みのクラス分類モデルのインスタンスを指定する。predict()メソッドを持つこと（引数は2次元配列を想定）。
ax: 決定境界を描くAxesオブジェクト。
x0s, x1s: クラスを計算する領域の計算点の座標を1次元配列で指定。
n_areas: 分割される領域の数を整数で指定。デフォルトは2（2つの領域）
colors: 分割される領域を塗りつぶす色をカラーコードの配列で与える。デフォルトは['tab:blue', 'tab:oranbe']。
alpha: 分布図の場合の塗りつぶしの透明度を実数で指定。デフォルトは0.5（半透明）。

関数の内容

`draw_decision_boundary()`

pyplotのcontourを利用している。

def draw_decision_boundary(clf, ax, x0s, x1s,
        threshold=0, color='k', alpha=1.0):

    y_predicted = np.empty((len(x1s), len(x0s)))

    for row, x1 in enumerate(x1s):
        for col, x0 in enumerate(x0s):
            y_predicted[row, col] = clf.predict(np.array([[x0, x1]]))

    ax.contour(x0s, x1s, y_predicted,
        colors=color, levels=[threshold] , alpha=alpha)

def draw_decision_boundary(clf, ax, x0s, x1s,

threshold=0, color='k', alpha=1.0):

y_predicted = np.empty((len(x1s), len(x0s)))

for row, x1 in enumerate(x1s):

for col, x0 in enumerate(x0s):

y_predicted[row, col] = clf.predict(np.array([[x0, x1]]))

ax.contour(x0s, x1s, y_predicted,

colors=color, levels=[threshold] , alpha=alpha)

`draw_decision_field()`

pyplotのcontourfを利用している。

def draw_decision_field(clf, ax, x0s, x1s, n_areas=2,
        colors=['tab:blue', 'tab:orange'], alpha=0.5, fill=True):

    y_predicted = np.empty((len(x1s), len(x0s)))

    for row, x1 in enumerate(x1s):
        for col, x0 in enumerate(x0s):
            y_predicted[row, col] = clf.predict(np.array([[x0, x1]]))
    ax.contourf(x0s, x1s, y_predicted,
        colors=colors, levels=n_areas-1, alpha=alpha)

def draw_decision_field(clf, ax, x0s, x1s, n_areas=2,

colors=['tab:blue', 'tab:orange'], alpha=0.5, fill=True):

y_predicted = np.empty((len(x1s), len(x0s)))

for row, x1 in enumerate(x1s):

for col, x0 in enumerate(x0s):

y_predicted[row, col] = clf.predict(np.array([[x0, x1]]))

ax.contourf(x0s, x1s, y_predicted,

colors=colors, levels=n_areas-1, alpha=alpha)

決定木の境界描画

2020-06-10 / tau / コメントする

概要

書籍”Pythonではじめる機械学習”の決定木のところで、ノードの分割をするごとの境界を描いている。

書籍ではmglearnパッケージを使っているが、これを自前の関数で再現した例。

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patch
from sklearn.datasets import make_moons
from sklearn.tree import DecisionTreeClassifier


def draw_tree_boundary(tree, ax, left, right, bottom, top,
        i_node=0, stop_level=None, n_level=0):

    if tree.children_left[i_node] == -1 or stop_level == n_level:
        fc =\
            'tab:orange' if np.argmax(tree.value[i_node][0])==0 else 'tab:blue'
        rect = patch.Rectangle(xy=(left, bottom),
            width=right-left, height=top-bottom, fc=fc, alpha=0.2)
        ax.add_patch(rect)
        return

    if tree.feature[i_node] == 0:
        f0 = tree.threshold[i_node]
        ax.plot([f0, f0], [top, bottom])
        draw_tree_boundary(tree=tree, ax=ax,
            left=left, right=f0, top=top, bottom=bottom,
            i_node=tree.children_left[i_node],
            stop_level=stop_level, n_level=n_level+1,)
        draw_tree_boundary(tree=tree, ax=ax,
            left=f0, right=right, top=top, bottom=bottom,
            i_node=tree.children_right[i_node],
            stop_level=stop_level, n_level=n_level+1)
    else:
        f1 = tree.threshold[i_node]
        ax.plot([left, right], [f1, f1])
        draw_tree_boundary(tree=tree, ax=ax,
            left=left, right=right, top=f1, bottom=bottom,
            i_node=tree.children_left[i_node],
            stop_level=stop_level, n_level=n_level+1)
        draw_tree_boundary(tree=tree, ax=ax,
            left=left, right=right, top=top, bottom=f1,
            i_node=tree.children_right[i_node],
            stop_level=stop_level, n_level=n_level+1)


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

treeclf = \
    DecisionTreeClassifier(max_depth=None, min_samples_leaf=1, random_state=0)
treeclf.fit(X, y)

fig, ax = plt.subplots()

ax.scatter(X[y==0][:, 0], X[y==0][:, 1],
    ec='k', s=60, marker='o', fc='tab:orange', label="Class 0")
ax.scatter(X[y==1][:, 0], X[y==1][:, 1],
    ec='k', s=60, marker='^', fc='tab:blue', label="Class 1")

x0_min, x0_max = -2, 2.5
x1_min, x1_max = -1, 1.5

draw_tree_boundary(tree=treeclf.tree_, i_node=0, ax=ax,
    left=x0_min, right=x0_max, bottom=x1_min, top=x1_max)

ax.set_xlim(x0_min, x0_max)
ax.set_ylim(x1_min, x1_max)
ax.legend()

plt.show()

import numpy as np

import matplotlib.pyplot as plt

import matplotlib.patches as patch

from sklearn.datasets import make_moons

from sklearn.tree import DecisionTreeClassifier

def draw_tree_boundary(tree, ax, left, right, bottom, top,

i_node=0, stop_level=None, n_level=0):

if tree.children_left[i_node] == -1 or stop_level == n_level:

fc =\

'tab:orange' if np.argmax(tree.value[i_node][0])==0 else 'tab:blue'

rect = patch.Rectangle(xy=(left, bottom),

width=right-left, height=top-bottom, fc=fc, alpha=0.2)

ax.add_patch(rect)

return

if tree.feature[i_node] == 0:

f0 = tree.threshold[i_node]

ax.plot([f0, f0], [top, bottom])

draw_tree_boundary(tree=tree, ax=ax,

left=left, right=f0, top=top, bottom=bottom,

i_node=tree.children_left[i_node],

stop_level=stop_level, n_level=n_level+1,)

draw_tree_boundary(tree=tree, ax=ax,

left=f0, right=right, top=top, bottom=bottom,

i_node=tree.children_right[i_node],

stop_level=stop_level, n_level=n_level+1)

else:

f1 = tree.threshold[i_node]

ax.plot([left, right], [f1, f1])

draw_tree_boundary(tree=tree, ax=ax,

left=left, right=right, top=f1, bottom=bottom,

i_node=tree.children_left[i_node],

stop_level=stop_level, n_level=n_level+1)

draw_tree_boundary(tree=tree, ax=ax,

left=left, right=right, top=top, bottom=f1,

i_node=tree.children_right[i_node],

stop_level=stop_level, n_level=n_level+1)

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

treeclf = \

DecisionTreeClassifier(max_depth=None, min_samples_leaf=1, random_state=0)

treeclf.fit(X, y)

fig, ax = plt.subplots()

ax.scatter(X[y==0][:, 0], X[y==0][:, 1],

ec='k', s=60, marker='o', fc='tab:orange', label="Class 0")

ax.scatter(X[y==1][:, 0], X[y==1][:, 1],

ec='k', s=60, marker='^', fc='tab:blue', label="Class 1")

x0_min, x0_max = -2, 2.5

x1_min, x1_max = -1, 1.5

draw_tree_boundary(tree=treeclf.tree_, i_node=0, ax=ax,

left=x0_min, right=x0_max, bottom=x1_min, top=x1_max)

ax.set_xlim(x0_min, x0_max)

ax.set_ylim(x1_min, x1_max)

ax.legend()

plt.show()

関数の仕様

描画用の関数draw_tree_boundary()の引数は以下の通り。

draw_tree_boundary(tree, ax, left, right, bottom, top, i_node=0, stop_level=None, n_level=0)

tree: 描きたい決定木モデルのtree_オブジェクトを渡す。
ax: 境界図を描くターゲットのAxesオブジェクトを渡す。
left, right, bottom, top: その時点でのノードの描画範囲をaxに即した座標で指定する。
i_node: エリアを描画するノード。省略した場合のデフォルトは0で、ルートノード（全域）以下を描画。
stop_level: 描画する木の深さを指定。デフォルトはNoneで、この場合は最深部まで描く。
n_level: この関数の再帰呼び出しの際に内部的に使われる。

この関数を呼び出し方の例は以下の通りで、stop_levelを省略しているので、リーフノードまで含めた木全体を描いている。

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

treeclf = \
    DecisionTreeClassifier(max_depth=None, min_samples_leaf=1, random_state=0)
treeclf.fit(X, y)

fig, ax = plt.subplots()

ax.scatter(X[y==0][:, 0], X[y==0][:, 1],
    ec='k', s=60, marker='o', fc='tab:orange', label="Class 0")
ax.scatter(X[y==1][:, 0], X[y==1][:, 1],
    ec='k', s=60, marker='^', fc='tab:blue', label="Class 1")

x0_min, x0_max = -2, 2.5
x1_min, x1_max = -1, 1.5

draw_tree_boundary(tree=treeclf.tree_, i_node=0, ax=ax,
    left=x0_min, right=x0_max, bottom=x1_min, top=x1_max)

ax.set_xlim(x0_min, x0_max)
ax.set_ylim(x1_min, x1_max)
ax.legend()

plt.show()

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

treeclf = \

DecisionTreeClassifier(max_depth=None, min_samples_leaf=1, random_state=0)

treeclf.fit(X, y)

fig, ax = plt.subplots()

ax.scatter(X[y==0][:, 0], X[y==0][:, 1],

ec='k', s=60, marker='o', fc='tab:orange', label="Class 0")

ax.scatter(X[y==1][:, 0], X[y==1][:, 1],

ec='k', s=60, marker='^', fc='tab:blue', label="Class 1")

x0_min, x0_max = -2, 2.5

x1_min, x1_max = -1, 1.5

draw_tree_boundary(tree=treeclf.tree_, i_node=0, ax=ax,

left=x0_min, right=x0_max, bottom=x1_min, top=x1_max)

ax.set_xlim(x0_min, x0_max)

ax.set_ylim(x1_min, x1_max)

ax.legend()

plt.show()

関数の処理内容

この関数の大まかな処理の流れは、以下の通り。

ルートノードの分割から初めて、リーフノードに行きつくまで分割と下の階層の探索を再帰的に進める
リーフノードであればそのノードのクラスで色を塗り、親のノードに戻る
あるノードの左の子ノードの下のリーフノードの処理が全部終わったら、右の子ノードの処理に移り、それも終わったら親のノードに戻る

関数の処理内容を最初の呼び出しから追うと以下の通り。

i_nodeとstop_levelを省略して呼び出し→ルートノードから木全体を描く
現在のノードがリーフノード（子ノードのインデックスが–1）あるいは現在の深さがstop_levelに達したなら、以下を実行してreturn（親ノードに戻る）
1. 現在のノードの卓越クラスに応じてtab:orangeかtab:blueでフェイスカラーを設定
2. 引数で得られた矩形領域をフェイスカラーで塗りつぶす
3. 塗りつぶした矩形をaxに追加
現在のノードがリーフノードでなく、終了深さにも達していない場合は、現在のノードを分割する特徴量によって以下を実行してreturn（親ノードに戻る）
1. ノードの分割基準が特徴量0の場合
  1. 分割基準の特徴量0の値で領域の上から下まで境界線を引く
  2. 左側のエリアを指定して左子ノードを処理する
  3. 戻ってきたら右側のエリアを指定して右子ノードを処理する
2. ノードの分割基準が特徴量1の場合
  1. 分割基準の特徴量1の値で領域の左から右まで境界線を引く
  2. 下側のエリアを指定して左子ノードを処理する
  3. 戻ってきたら上側のエリアを指定して右子ノードを処理する

def draw_tree_boundary(tree, ax, left, right, bottom, top,
        i_node=0, stop_level=None, n_level=0):

    if tree.children_left[i_node] == -1 or stop_level == n_level:
        fc =\
            'tab:orange' if np.argmax(tree.value[i_node][0])==0 else 'tab:blue'
        rect = patch.Rectangle(xy=(left, bottom),
            width=right-left, height=top-bottom, fc=fc, alpha=0.2)
        ax.add_patch(rect)
        return

    if tree.feature[i_node] == 0:
        f0 = tree.threshold[i_node]
        ax.plot([f0, f0], [top, bottom])
        draw_tree_boundary(tree=tree, ax=ax,
            left=left, right=f0, top=top, bottom=bottom,
            i_node=tree.children_left[i_node],
            stop_level=stop_level, n_level=n_level+1,)
        draw_tree_boundary(tree=tree, ax=ax,
            left=f0, right=right, top=top, bottom=bottom,
            i_node=tree.children_right[i_node],
            stop_level=stop_level, n_level=n_level+1)
    else:
        f1 = tree.threshold[i_node]
        ax.plot([left, right], [f1, f1])
        draw_tree_boundary(tree=tree, ax=ax,
            left=left, right=right, top=f1, bottom=bottom,
            i_node=tree.children_left[i_node],
            stop_level=stop_level, n_level=n_level+1)
        draw_tree_boundary(tree=tree, ax=ax,
            left=left, right=right, top=top, bottom=f1,
            i_node=tree.children_right[i_node],
            stop_level=stop_level, n_level=n_level+1)

def draw_tree_boundary(tree, ax, left, right, bottom, top,

i_node=0, stop_level=None, n_level=0):

if tree.children_left[i_node] == -1 or stop_level == n_level:

fc =\

'tab:orange' if np.argmax(tree.value[i_node][0])==0 else 'tab:blue'

rect = patch.Rectangle(xy=(left, bottom),

width=right-left, height=top-bottom, fc=fc, alpha=0.2)

ax.add_patch(rect)

return

if tree.feature[i_node] == 0:

f0 = tree.threshold[i_node]

ax.plot([f0, f0], [top, bottom])

draw_tree_boundary(tree=tree, ax=ax,

left=left, right=f0, top=top, bottom=bottom,

i_node=tree.children_left[i_node],

stop_level=stop_level, n_level=n_level+1,)

draw_tree_boundary(tree=tree, ax=ax,

left=f0, right=right, top=top, bottom=bottom,

i_node=tree.children_right[i_node],

stop_level=stop_level, n_level=n_level+1)

else:

f1 = tree.threshold[i_node]

ax.plot([left, right], [f1, f1])

draw_tree_boundary(tree=tree, ax=ax,

left=left, right=right, top=f1, bottom=bottom,

i_node=tree.children_left[i_node],

stop_level=stop_level, n_level=n_level+1)

draw_tree_boundary(tree=tree, ax=ax,

left=left, right=right, top=top, bottom=f1,

i_node=tree.children_right[i_node],

stop_level=stop_level, n_level=n_level+1)

決定木による回帰

2020-06-09 / tau / コメントする

概要

決定木を回帰に用いる場合、回帰木(regression tree)とも呼ぶ。ここでは決定木の回帰における性質・挙動を確認する。

回帰木の学習過程

以下は、sin関数に対して回帰木を適用し、剪定の深さを深くしていった場合の推移。

剪定深さ1の場合、特徴量を2つに分割しそれぞれの領域のデータから学習し予測値を得ている。剪定深さ2の場合、さらに各領域を2分割して4つの領域で予測値を得ている。このようにして剪定深さnに対して2ⁿの領域のデータで学習する。この例の場合は訓練セットとして80個のデータを準備し、1000個のデータの予測をしている。

剪定深さ6で2⁶=64訓練セットの個数と近くなるが、サインカーブの山と谷のところで区間が長く、誤差が出ている。これは、回帰木のノード分割がy = sin xの値に基づいて行われるとき、その値がかなり近くなる山・谷のところでなかなか分離されないからと考えられる。

剪定深さ10で2¹⁰=1024のとき、分割数がテストセットと同じくらいの数になるので初めて値が近い点も区別され、全体がフィットする。

import numpy as np
import matplotlib.pyplot as plt
from pandas import DataFrame
from sklearn.tree import DecisionTreeRegressor

x_train = np.linspace(0, np.pi * 2, num=80)
y_train = np.sin(x_train).reshape(-1, 1)
X_train = x_train.reshape(-1, 1)

x_test = np.linspace(0, np.pi * 2, num=1000).reshape(-1, 1)
X_test = x_test.reshape(-1, 1)

depths = [1, 2, 3, 4, 6, 10]

fig, axs = plt.subplots(2, 3, figsize=(12, 6))
fig.subplots_adjust(hspace=0.3)
ax_1d = axs.reshape(-1)

for ax, depth in zip(ax_1d, depths):
    ax.scatter(x_train, y_train, s=10, ec='gray', fc="coral")
    treereg = DecisionTreeRegressor(max_depth=depth).fit(X_train, y_train)
    y_pred = treereg.predict(X_test)
    ax.plot(x_test, y_pred, c='blue', linewidth=1)
    ax.set_title("max_depth={}".format(depth))

plt.show()

import numpy as np

import matplotlib.pyplot as plt

from pandas import DataFrame

from sklearn.tree import DecisionTreeRegressor

x_train = np.linspace(0, np.pi * 2, num=80)

y_train = np.sin(x_train).reshape(-1, 1)

X_train = x_train.reshape(-1, 1)

x_test = np.linspace(0, np.pi * 2, num=1000).reshape(-1, 1)

X_test = x_test.reshape(-1, 1)

depths = [1, 2, 3, 4, 6, 10]

fig, axs = plt.subplots(2, 3, figsize=(12, 6))

fig.subplots_adjust(hspace=0.3)

ax_1d = axs.reshape(-1)

for ax, depth in zip(ax_1d, depths):

ax.scatter(x_train, y_train, s=10, ec='gray', fc="coral")

treereg = DecisionTreeRegressor(max_depth=depth).fit(X_train, y_train)

y_pred = treereg.predict(X_test)

ax.plot(x_test, y_pred, c='blue', linewidth=1)

ax.set_title("max_depth={}".format(depth))

plt.show()

ここで、学習途上の状況を、剪定深さ2(max_depth=2)の時の状態で確認してみる。

分割された4つの領域に対する境界(0.517, 3.142, 5.766)のうち、最初の境界3.142はπの値で、0～2πの領域においてπの両側で対称なことから自然な結果。4つの領域における予測値(value)はグラフ上でも確認でき、やはりπの両側で対称な値となっている。

import numpy as np
import matplotlib.pyplot as plt
import graphviz
from sklearn.tree import DecisionTreeRegressor, export_graphviz

x_train = np.linspace(0, np.pi * 2, num=80)
y_train = np.sin(x_train)
X_train = x_train.reshape(-1, 1)

treereg = DecisionTreeRegressor(max_depth=2).fit(X_train, y_train)

x_test = np.linspace(0, np.pi * 2, num=1000)
X_test = x_test.reshape(-1, 1)
y_pred = treereg.predict(X_test)

dot_data = export_graphviz(treereg, out_file=None,
    feature_names=["x"], class_names=["sin"])
graph = graphviz.Source(dot_data)
graph.render("image", view=True)

fig, ax = plt.subplots()
ax.scatter(x_train, y_train, s=20, ec='gray', fc="coral")
ax.plot(X_test, y_pred, c='blue')
ax.set_xlabel("x")
plt.show()

import numpy as np

import matplotlib.pyplot as plt

import graphviz

from sklearn.tree import DecisionTreeRegressor, export_graphviz

x_train = np.linspace(0, np.pi * 2, num=80)

y_train = np.sin(x_train)

X_train = x_train.reshape(-1, 1)

treereg = DecisionTreeRegressor(max_depth=2).fit(X_train, y_train)

x_test = np.linspace(0, np.pi * 2, num=1000)

X_test = x_test.reshape(-1, 1)

y_pred = treereg.predict(X_test)

dot_data = export_graphviz(treereg, out_file=None,

feature_names=["x"], class_names=["sin"])

graph = graphviz.Source(dot_data)

graph.render("image", view=True)

fig, ax = plt.subplots()

ax.scatter(x_train, y_train, s=20, ec='gray', fc="coral")

ax.plot(X_test, y_pred, c='blue')

ax.set_xlabel("x")

plt.show()

ノイズの影響と過学習

先のサインカーブにノイズが乗った場合の回帰木を見てみる。剪定深さ3、5としたときの回帰木による回帰線の形は以下の通りで、深さが深いと個別のデータに対して過学習となっている様子がわかる。

これらのモデルのスコアは以下の通りで、深さ5の場合には訓練スコアに対してテストスコアが低く、過学習となっている。

深さ3の場合に訓練スコアの方がテストスコアより低いが、これは訓練スコアにノイズが含まれるのに対してテストスコアのyの値をすべてノイズがないsin値としているためで、訓練セットにおいて乱数を加える程度を小さくするとこの逆転現象は解消される。

depth=3
 training score: 0.944
 test score    : 0.951
depth=5
 training score: 0.979
 test score    : 0.942

depth=3

training score: 0.944

test score : 0.951

depth=5

training score: 0.979

test score : 0.942

上記の実行コードは以下の通り。

import numpy as np
import numpy.random as rnd
import matplotlib.pyplot as plt
from pandas import DataFrame
from sklearn.tree import DecisionTreeRegressor

n_train = 80
n_rand = 10
depth1 = 3
depth2 = 5

rnd.seed(23)

x = np.linspace(0, np.pi * 2, num=n_train)
y = np.sin(x)
for i in range(0, len(x), n_train//n_rand):
    y[i] = y[i] + (rnd.rand() * 2 - 1)
df = DataFrame(x, columns=['x'])
df['y'] = y

x_train = np.array(df['x'])
y_train = np.array(df['y'])
X_train = x_train.reshape(-1, 1)

treereg1 = DecisionTreeRegressor(max_depth=depth1).fit(X_train, y_train)
treereg2 = DecisionTreeRegressor(max_depth=depth2).fit(X_train, y_train)

x_test = np.linspace(0, np.pi * 2, 1000)
y_test = np.sin(x_test)
X_test = x_test.reshape(-1, 1)
y_pred1 = treereg1.predict(X_test)
y_pred2 = treereg2.predict(X_test)

print("depth={}".format(depth1))
print(" training score:{:6.3f}".format(treereg1.score(X_train, y_train)))
print(" test score    :{:6.3f}".format(treereg1.score(X_test, y_test)))
print("depth={}".format(depth2))
print(" training score:{:6.3f}".format(treereg2.score(X_train, y_train)))
print(" test score    :{:6.3f}".format(treereg2.score(X_test, y_test)))

fig, ax = plt.subplots()

ax.scatter(x_train, y_train, s=20, ec='gray', fc="coral")
ax.plot(X_test, y_pred1,
    c='forestgreen', linestyle='dashed', label="depth={}".format(depth1))
ax.plot(X_test, y_pred2, c='blue', label="depth={}".format(depth2))
ax.legend()

plt.show()

import numpy as np

import numpy.random as rnd

import matplotlib.pyplot as plt

from pandas import DataFrame

from sklearn.tree import DecisionTreeRegressor

n_train = 80

n_rand = 10

depth1 = 3

depth2 = 5

rnd.seed(23)

x = np.linspace(0, np.pi * 2, num=n_train)

y = np.sin(x)

for i in range(0, len(x), n_train//n_rand):

y[i] = y[i] + (rnd.rand() * 2 - 1)

df = DataFrame(x, columns=['x'])

df['y'] = y

x_train = np.array(df['x'])

y_train = np.array(df['y'])

X_train = x_train.reshape(-1, 1)

treereg1 = DecisionTreeRegressor(max_depth=depth1).fit(X_train, y_train)

treereg2 = DecisionTreeRegressor(max_depth=depth2).fit(X_train, y_train)

x_test = np.linspace(0, np.pi * 2, 1000)

y_test = np.sin(x_test)

X_test = x_test.reshape(-1, 1)

y_pred1 = treereg1.predict(X_test)

y_pred2 = treereg2.predict(X_test)

print("depth={}".format(depth1))

print(" training score:{:6.3f}".format(treereg1.score(X_train, y_train)))

print(" test score :{:6.3f}".format(treereg1.score(X_test, y_test)))

print("depth={}".format(depth2))

print(" training score:{:6.3f}".format(treereg2.score(X_train, y_train)))

print(" test score :{:6.3f}".format(treereg2.score(X_test, y_test)))

fig, ax = plt.subplots()

ax.scatter(x_train, y_train, s=20, ec='gray', fc="coral")

ax.plot(X_test, y_pred1,

c='forestgreen', linestyle='dashed', label="depth={}".format(depth1))

ax.plot(X_test, y_pred2, c='blue', label="depth={}".format(depth2))

ax.legend()

plt.show()

同じデータで決定木の剪定深さを変えていったときの状況を如何に示す。訓練スコアとテストスコアの関係から、深さ3までは学習不足、深さ4以降は過学習となっていることが示され、過学習になるとノイズの影響を受けていることがわかる。

決定木の限界～外挿

決定木は、与えられた訓練データに対しては完全な予測も可能だが、訓練データの領域外のデータに対しては妥当な予測ができない。書籍”Pythonではじめる機械学習”で紹介されている、メモリー単価の推移によってこれを確認する（データについてはこちらのサイトのものを使わせてもらった）。

時間をx、メモリー単価をyとするとメモリー単価を対数で表したlog yはxに対して概ね線形関係になっている。以下は、縦軸を対数目盛とした場合のメモリー単価、xとlog yについて線形回帰と決定木による学習と予測の結果を示したもので、2000年より前のデータによって双方のモデルを学習させ、2000年以降の価格を予測している。

線形回帰はデータの細かい傾向は再現できないが、訓練セットの外側についてもその傾向をある程度予測できている。一方決定木については、訓練セットについては完全に予測しているが、その外側になった途端に、外側の直前のデータの値をそのまま予測値としている。

import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.tree import DecisionTreeRegressor

data_path = \
    r"C:\Users\tomo\GoogleDrive\IT_and_Mobile\dev\python\Machine_Learning"
memory_prices = \
    pd.read_csv(os.path.join(data_path, r"regression_tree\memory_prices.csv"))

data_train = memory_prices[memory_prices.date < 2000]
data_test = memory_prices[memory_prices.date >= 2000]

x_train = np.array(data_train.date)
x_test = np.array(data_test.date)
y_train = np.array(data_train.price)
y_test = np.array(data_test.price)
X_train = x_train.reshape(-1, 1)
X_test = x_test.reshape(-1, 1)

linreg = LinearRegression().fit(X_train, np.log(y_train))
y_pred_linreg_train = np.exp(linreg.predict(X_train))
y_pred_linreg_test = np.exp(linreg.predict(X_test))

treereg = DecisionTreeRegressor().fit(X_train, np.log(y_train))
y_pred_treereg_train = np.exp(treereg.predict(X_train))
y_pred_treereg_test = np.exp(treereg.predict(X_test))

fig, ax = plt.subplots()
ax.plot(data_train.date, data_train.price, c='tab:blue', label="memory prices")
ax.plot(data_test.date, data_test.price, c='tab:blue', linestyle='dashed')
ax.plot(x_train, y_pred_linreg_train, c='tab:green', label="Linear Regression")
ax.plot(x_test, y_pred_linreg_test, c='tab:green', linestyle='dashed')
ax.plot(x_train, y_pred_treereg_train,
    c='tab:red', linestyle='dotted', linewidth=3, label="Decision Tree")
ax.plot(x_test, y_pred_treereg_test, c='tab:red', linestyle='dashed')
ax.set_yscale('log')
ax.legend()
plt.show()

import os

import numpy as np

import pandas as pd

import matplotlib.pyplot as plt

from sklearn.linear_model import LinearRegression

from sklearn.tree import DecisionTreeRegressor

data_path = \

r"C:\Users\tomo\GoogleDrive\IT_and_Mobile\dev\python\Machine_Learning"

memory_prices = \

pd.read_csv(os.path.join(data_path, r"regression_tree\memory_prices.csv"))

data_train = memory_prices[memory_prices.date < 2000]

data_test = memory_prices[memory_prices.date >= 2000]

x_train = np.array(data_train.date)

x_test = np.array(data_test.date)

y_train = np.array(data_train.price)

y_test = np.array(data_test.price)

X_train = x_train.reshape(-1, 1)

X_test = x_test.reshape(-1, 1)

linreg = LinearRegression().fit(X_train, np.log(y_train))

y_pred_linreg_train = np.exp(linreg.predict(X_train))

y_pred_linreg_test = np.exp(linreg.predict(X_test))

treereg = DecisionTreeRegressor().fit(X_train, np.log(y_train))

y_pred_treereg_train = np.exp(treereg.predict(X_train))

y_pred_treereg_test = np.exp(treereg.predict(X_test))

fig, ax = plt.subplots()

ax.plot(data_train.date, data_train.price, c='tab:blue', label="memory prices")

ax.plot(data_test.date, data_test.price, c='tab:blue', linestyle='dashed')

ax.plot(x_train, y_pred_linreg_train, c='tab:green', label="Linear Regression")

ax.plot(x_test, y_pred_linreg_test, c='tab:green', linestyle='dashed')

ax.plot(x_train, y_pred_treereg_train,

c='tab:red', linestyle='dotted', linewidth=3, label="Decision Tree")

ax.plot(x_test, y_pred_treereg_test, c='tab:red', linestyle='dashed')

ax.set_yscale('log')

ax.legend()

plt.show()

決定木によるクラス分類

2020-06-06 / tau / コメントする

概要

決定木をクラス分類に用いる場合、分類木(classification tree)とも呼ぶ。決定木は、決定境界が単調ではなく混み入っていて線形モデルでは分類が難しい場合でも対応できる。ここでは決定木のクラス分類における性質・挙動を確認する。

決定木の構築過程

2つの特徴量を持ち、2つのクラスのいずれかに属するデータについて、決定木が作られていく過程を追っていく。

データとしては、scikit-learnのmake_moons()で得られる以下のデータを用いる。

import matplotlib.pyplot as plt
from sklearn.datasets import make_moons

X, y = make_moons(n_samples=30, noise=0.25, random_state=9)

fig, ax = plt.subplots()

ax.scatter(X[y==0][:, 0], X[y==0][:, 1],
    ec='k', s=60, marker='o', label="Class 0")
ax.scatter(X[y==1][:, 0], X[y==1][:, 1],
    ec='k', s=60, marker='^', label="Class 1")

ax.set_xlabel("Feature-0")
ax.set_ylabel("Feature-1")
ax.legend()

plt.show()

import matplotlib.pyplot as plt

from sklearn.datasets import make_moons

X, y = make_moons(n_samples=30, noise=0.25, random_state=9)

fig, ax = plt.subplots()

ax.scatter(X[y==0][:, 0], X[y==0][:, 1],

ec='k', s=60, marker='o', label="Class 0")

ax.scatter(X[y==1][:, 0], X[y==1][:, 1],

ec='k', s=60, marker='^', label="Class 1")

ax.set_xlabel("Feature-0")

ax.set_ylabel("Feature-1")

ax.legend()

plt.show()

このデータセットについて、2つの特徴量のいずれかを調節して順次領域を分割していく。このとき、どのように分割するのが最も妥当かということについては、決定木の分割の考え方を参照。また、以下の実行例のコードについてはDecisionTreeClassifierに関するテストプログラムを参照。

まず最初の分割は以下の通りで、特徴量1の値0.272が境界となり、それ以下がクラス0、それより大きいとクラス1が卓越していると判定される。

第2ステップは、第1ステップで分けられたそれぞれの領域を分割する。どちらの領域も特徴量0が境界となっていて、それぞれの領域／ノードの分布状態に応じた特性量によって左右に分けられている。

第3ステップでは、左上の領域は特徴量1で上下に、右下の領域は特徴量0で左右に分割されるが、いずれについても分割後の領域のクラスが同じになっている。利得が最も高くなるように分割しても、領域の中の擾乱クラスのデータが少ない場合はこのようになる。

なお、右上の領域はクラス1のデータが2個のみ、左下の領域クラス0のデータが2個のみと単一のクラスのデータしかないため、これ以上分割されない葉(leaf)となっている。

第4ステップでは、左上と右下の領域がいずれも特徴量0で左右に分割され、今度は分割後のノードが異なるクラスになっている。左上と右下の領域（ノード）がクラス1の葉となっている。

第5ステップでは下方右から2番目の領域を特徴量1で分割している。

分割はここで終了。全てのノードが単一のクラスで構成された純粋な状態となっている。

全体としてはクラス0が左側の上に凸な分布、クラス1が左側の下に凸な分布と判定されていて、make_moons()が意図した分布と合っている。ただし左上に一つクラス1のデータがあるために、本来意図しない領域がクラス1と分割されている。

過学習の抑制

過学習の視覚的な例

決定木の構築過程で、1つのデータの影響で予想と異なるノード区分が発生するのが見られた。別のパターンのデータの例を以下に示す。少ないデータの影響で領域分割が複雑になっていることがわかる。

このような場合、特定のパターンの教師データに対しては適合度が高く、純粋な葉のみで構成される場合には適合度が100%になるが、他のデータに対しての適合度を下げてしまう。いわゆる過学習となってしまう。

決定木の過学習を防ぐための方法として、決定木の階層をあるレベルまでに留める、いわゆる枝刈り／剪定(pruning)という考え方や、葉の切り分けを純粋なレベルより前に留める（ノード内のデータ数が2つ以上複数のデータ数になったら分割を止める）という考え方がある。

これらはscikit-learnのDecisionTreeClassifierでは、枝刈りのうち事前剪定(pre-pruning)のためのパラメーターとしてmax_depth、ノードの最小データ数のパラメーターとしてmin_samples_leafが設定できる。

max_depth～剪定

以下の例は、上記のデータに対してmax_depthを変化させたときの領域分割の様子。このパラメーターはデフォルトではNoneで可能な限り分割を行っていくが、正の整数値を指定すると、その深さまでで分割を止める。分割の深さが少なくなるにしたがってモデルが単純化されていく様子がわかる。

min_samples_leaf～葉の純度の制限

min_samples_leafはデフォルトでは1で、ノードが完全に純粋でない限りデータ数が1個になるまで分割を試みる。この値を変化させて、異なるクラスのデータを含んでいても分割を行わないようにした場合の領域分割の状況を見てみる。

枝刈りに比べてモデルの複雑さは回避しながら、それらしい領域分割になっている。これは、事前枝刈りが各ノードの不純度に関わらず同じレベルで計算を止めるのに対して、葉の純度を個別にコントロールしているために柔軟に分割が進められているためと考えられる。

cancerデータによる確認

剪定による過剰適合の抑制

breast_cancerデータセットに対してDecisionTreeClassifierを適用してクラス分類し、訓練セットとテストセットのスコアを計算する。リーフノードが純粋になるまで木を成長させた場合と、深さ4で事前剪定をした場合のスコアを比較してみる。

import pydotplus as pdp
from sklearn.model_selection import train_test_split
from sklearn.datasets import load_breast_cancer
from sklearn.tree import DecisionTreeClassifier, export_graphviz

ds = load_breast_cancer()
X_train, X_test, y_train, y_test = train_test_split(ds.data, ds.target, stratify=ds.target, random_state=42)

clf = DecisionTreeClassifier(random_state=0)
clf.fit(X_train, y_train)

print("Tree depth: {}".format(clf.get_depth()))
print("Training score:{:6.3f}".format(clf.score(X_train, y_train)))
print("Test score    :{:6.3f}".format(clf.score(X_test, y_test)))

print()

clf4 = DecisionTreeClassifier(max_depth=4, random_state=0)
clf4.fit(X_train, y_train)

print("Tree depth: {}".format(clf4.get_depth()))
print("Training score:{:6.3f}".format(clf4.score(X_train, y_train)))
print("Test score    :{:6.3f}".format(clf4.score(X_test, y_test)))

import pydotplus as pdp

from sklearn.model_selection import train_test_split

from sklearn.datasets import load_breast_cancer

from sklearn.tree import DecisionTreeClassifier, export_graphviz

ds = load_breast_cancer()

X_train, X_test, y_train, y_test = train_test_split(ds.data, ds.target, stratify=ds.target, random_state=42)

clf = DecisionTreeClassifier(random_state=0)

clf.fit(X_train, y_train)

print("Tree depth: {}".format(clf.get_depth()))

print("Training score:{:6.3f}".format(clf.score(X_train, y_train)))

print("Test score :{:6.3f}".format(clf.score(X_test, y_test)))

print()

clf4 = DecisionTreeClassifier(max_depth=4, random_state=0)

clf4.fit(X_train, y_train)

print("Tree depth: {}".format(clf4.get_depth()))

print("Training score:{:6.3f}".format(clf4.score(X_train, y_train)))

print("Test score :{:6.3f}".format(clf4.score(X_test, y_test)))

出力結果は以下の通り。完全な木は7層で、訓練セットに対しては全データに適合しており、テストセットに対しては93.7%の適合率。一方、深さ4で枝刈りをした場合は、訓練セットに対する適合率は下がるがテストセットの適合率は上がり、過剰適合が抑制されている。

Tree depth: 7
Training score: 1.000
Test score    : 0.937

Tree depth: 4
Training score: 0.988
Test score    : 0.951

Tree depth: 7

Training score: 1.000

Test score : 0.937

Tree depth: 4

Training score: 0.988

Test score : 0.951

剪定の深さに対するスコアの変化

事前剪定の深さレベルを変化させたときの、訓練セットとテストセットに対するスコアの変化を確認する。train_test_split()の乱数系列によって結果のパターンが異なるが、概ねmax_depth=4で適合不足と過剰適合のバランスが最もとれているようであり、スコアは0.95程度。

ただし、そもそも決定木の深さがそれほど深くなく、max_depthのバリエーションが数個となるので、線形モデルにおけるハイパーパラメーターのような連続的な曲線は描き難い。

max_depths = np.arange(8, 0, -1)
training_scores = np.empty(max_depths.size)
test_scores = np.empty(max_depths.size)

fig, axs = plt.subplots(2, 2, figsize=(8, 6))
fig.subplots_adjust(hspace=0.4)
ax_1d = axs.reshape(axs.size)

random_state_list = [0, 1, 2, 3]
for random_state, ax in zip(random_state_list, ax_1d):
    X_train, X_test, y_train, y_test =\
        train_test_split(ds.data, ds.target, random_state=random_state)
    for i, max_depth in enumerate(max_depths):
        clf = DecisionTreeClassifier(max_depth=max_depth, random_state=0)
        clf.fit(X_train, y_train)
        training_scores[i] = clf.score(X_train, y_train)
        test_scores[i] = clf.score(X_test, y_test)
    ax.plot(max_depths, training_scores, clip_on=False)
    ax.plot(max_depths, test_scores)
    ax.set_xlim(ax.get_xlim()[::-1])
    ax.set_ylim(0.85, 1)
    ax.set_xlabel("max_depth")
    ax.set_title("random_state for split:{}".format(random_state))

plt.show()

max_depths = np.arange(8, 0, -1)

training_scores = np.empty(max_depths.size)

test_scores = np.empty(max_depths.size)

fig, axs = plt.subplots(2, 2, figsize=(8, 6))

fig.subplots_adjust(hspace=0.4)

ax_1d = axs.reshape(axs.size)

random_state_list = [0, 1, 2, 3]

for random_state, ax in zip(random_state_list, ax_1d):

X_train, X_test, y_train, y_test =\

train_test_split(ds.data, ds.target, random_state=random_state)

for i, max_depth in enumerate(max_depths):

clf = DecisionTreeClassifier(max_depth=max_depth, random_state=0)

clf.fit(X_train, y_train)

training_scores[i] = clf.score(X_train, y_train)

test_scores[i] = clf.score(X_test, y_test)

ax.plot(max_depths, training_scores, clip_on=False)

ax.plot(max_depths, test_scores)

ax.set_xlim(ax.get_xlim()[::-1])

ax.set_ylim(0.85, 1)

ax.set_xlabel("max_depth")

ax.set_title("random_state for split:{}".format(random_state))

plt.show()

リーフの最小サンプル数に対するスコアの変化

葉ノードの最小サンプル数を一定値以上とするmin_samples_leafを変化させたときの、訓練セットとテストセットの変化を見てみたが、乱数系列によってけっこうパターンがばらついている。この傾向は、max_depthを変化させても変わらなかった。

特徴量重要度

特徴量重要度の特性

breast_cancerデータセットを深さ4で剪定した決定木によってクラス分類した場合の特徴量重要度は以下のようになる。このグラフを表示するコードや特徴量重要度の計算方法についてはこちらを参照。

                 feature  importance
10          radius error    0.010197
11         texture error    0.048398
14      smoothness error    0.002416
20          worst radius    0.726829
21         worst texture    0.045816
24      worst smoothness    0.014158
26       worst concavity    0.018188
27  worst concave points    0.122113
28        worst symmetry    0.011885

feature importance

10 radius error 0.010197

11 texture error 0.048398

14 smoothness error 0.002416

20 worst radius 0.726829

21 worst texture 0.045816

24 worst smoothness 0.014158

26 worst concavity 0.018188

27 worst concave points 0.122113

28 worst symmetry 0.011885

このグラフと以下の決定木を比べ、重要度が大きい順に調べてみると以下のようなことがわかる。

worst radius(0.726)、第0層の不純度が高いノードを2つのクラスに分割し、分割後のノードの純度が高い
worst concave points(0.122)、第1層の左側のノード、259個のクラス1のデータの大部分を左側の子ノードに切り出しつつ、クラス0のデータも25を4と21と切り分けている
texture error(0.048)、第1層の右側のノード、クラス0のデータを完全に右の子ノードに切り分け、左の子ノードは不純度はゼロ
worst textureは(0.045)、第2層の左から2番目のノード、データ数は少ないが、左の子ノードにクラス1のデータを9/11、右の子ノードにクラス0のデータを18/21と子ノードの純度が高い
radius error(0.010)、第2層の左のノード、クラス1のデータを完全に左の子ノードに切り分け、右の子ノードの不純度はゼロ

特徴量重要度の計算方法や上記の特徴から、その性質は以下のように整理できる。

重要度は0～1の間の値をとる
浅いノードの対象特徴量ほど重要度が高い傾向（ただし、データの分布による可能性あり）
分割後の子ノードの純度が高いほど重要度が高い傾向
重要度の値は、どのクラスの切り分けに効いているかとは無関係

特徴量重要度の大きさは、枝を分離するときに重要な特徴量を示唆するが、その特徴量の大小とクラス分類の関係までは知ることができない。

線形モデルの特徴量係数との対比

特徴量重要度をLogistic回帰における特徴量の係数と比較してみる。以下はL2正則化によるLogistic回帰モデルをbreast_cancerデータに適用した場合の特徴量係数。

worst radiusについては、決定木での重要度が最も高いが、Logistic回帰でも比較的特徴量の重みは大きい。Logistic回帰の場合はこの特徴量がターゲット1（malignant：良性）であることを示唆しているが、決定木の場合にはそのような情報は得られない。

単調でないクラス分類

以下の例では、2つのクラス分類の境界が必ずしも1つではない。

Logistic回帰モデルの場合、一つの直線で分離しようとした結果、境界は赤い線のようになり、スコアも0.66625とかなり低い。

決定木はこのような場合でも分類可能だが、そもそもこのようなケースでは（それが決定木の性質ではなく本質的に）決定境界に対する大小だけでクラス分類を論ずることができない。

import numpy as np
import numpy.random as rnd
import matplotlib.pyplot as plt
import graphviz
from sklearn.tree import DecisionTreeClassifier, export_graphviz
from sklearn.linear_model import LogisticRegression

n0 = 40
n1a = 20
n1b = 20

rnd.seed(0)
f0_0 = rnd.rand(n0) * 4 + 3
f1_0 = rnd.rand(n0) * 1 + 4
f0_1a = rnd.rand(n1a) * 2 + 2
f1_1a = rnd.rand(n1a) * 1 + 2
f0_1b = rnd.rand(n1b) * 2 + 6
f1_1b = rnd.rand(n1b) * 1 + 6

X0 = np.hstack((f0_0.reshape(-1, 1), f1_0.reshape(-1, 1)))
X1a = np.hstack((f0_1a.reshape(-1, 1), f1_1a.reshape(-1, 1)))
X1b = np.hstack((f0_1b.reshape(-1, 1), f1_1b.reshape(-1, 1)))
X = np.vstack((X0, X1a, X1b))

y0 = np.zeros(n0)
y1a = np.ones(n1a)
y1b = np.ones(n1b)
y = np.hstack((y0, y1a, y1b))

clf = DecisionTreeClassifier().fit(X, y)
logreg = LogisticRegression().fit(X, y)
print(logreg.score(X, y))

dot_data= export_graphviz(clf, out_file=None,
    feature_names=["f0", "f1"], class_names=["C0", "C1"], filled=True)
graph = graphviz.Source(dot_data)
graph.render("image", view=True)

f0_min, f0_max = 1, 10
f1_min, f1_max = 1, 8

b = logreg.intercept_
w = logreg.coef_.reshape(logreg.coef_.size)
border_left = (-b - w[0] * f0_min) / w[1]
border_right = (-b - w[0] * f0_max) / w[1]

fig, ax = plt.subplots()
ax.scatter(X[y==0][:, 0], X[y==0][:, 1])
ax.scatter(X[y==1][:, 0], X[y==1][:, 1])
th1 = clf.tree_.threshold[0]
th2 = clf.tree_.threshold[2]
ax.plot([f0_min, f0_max], [th1, th1], c='tab:green', linestyle='dashed')
ax.plot([f0_min, f0_max], [th2, th2], c='tab:green', linestyle='dashed')
ax.plot([f0_min, f0_max], [border_left, border_right],
    c='tab:red', linestyle='dashed')
ax.set_xlim(f0_min, f0_max)
ax.set_ylim(f1_min, f1_max)
plt.show()

import numpy as np

import numpy.random as rnd

import matplotlib.pyplot as plt

import graphviz

from sklearn.tree import DecisionTreeClassifier, export_graphviz

from sklearn.linear_model import LogisticRegression

n0 = 40

n1a = 20

n1b = 20

rnd.seed(0)

f0_0 = rnd.rand(n0) * 4 + 3

f1_0 = rnd.rand(n0) * 1 + 4

f0_1a = rnd.rand(n1a) * 2 + 2

f1_1a = rnd.rand(n1a) * 1 + 2

f0_1b = rnd.rand(n1b) * 2 + 6

f1_1b = rnd.rand(n1b) * 1 + 6

X0 = np.hstack((f0_0.reshape(-1, 1), f1_0.reshape(-1, 1)))

X1a = np.hstack((f0_1a.reshape(-1, 1), f1_1a.reshape(-1, 1)))

X1b = np.hstack((f0_1b.reshape(-1, 1), f1_1b.reshape(-1, 1)))

X = np.vstack((X0, X1a, X1b))

y0 = np.zeros(n0)

y1a = np.ones(n1a)

y1b = np.ones(n1b)

y = np.hstack((y0, y1a, y1b))

clf = DecisionTreeClassifier().fit(X, y)

logreg = LogisticRegression().fit(X, y)

print(logreg.score(X, y))

dot_data= export_graphviz(clf, out_file=None,

feature_names=["f0", "f1"], class_names=["C0", "C1"], filled=True)

graph = graphviz.Source(dot_data)

graph.render("image", view=True)

f0_min, f0_max = 1, 10

f1_min, f1_max = 1, 8

b = logreg.intercept_

w = logreg.coef_.reshape(logreg.coef_.size)

border_left = (-b - w[0] * f0_min) / w[1]

border_right = (-b - w[0] * f0_max) / w[1]

fig, ax = plt.subplots()

ax.scatter(X[y==0][:, 0], X[y==0][:, 1])

ax.scatter(X[y==1][:, 0], X[y==1][:, 1])

th1 = clf.tree_.threshold[0]

th2 = clf.tree_.threshold[2]

ax.plot([f0_min, f0_max], [th1, th1], c='tab:green', linestyle='dashed')

ax.plot([f0_min, f0_max], [th2, th2], c='tab:green', linestyle='dashed')

ax.plot([f0_min, f0_max], [border_left, border_right],

c='tab:red', linestyle='dashed')

ax.set_xlim(f0_min, f0_max)

ax.set_ylim(f1_min, f1_max)

plt.show()

max_features～特徴量選択

DecisionTreeClassifierのコンストラクターのパラメーターの1つ、max_featuresについて。このパラメーターがデフォルトのNoneの場合やautoを指定した場合、n_featuresすなわちすべての特徴量が比較され、最も分離の性能がいいものが選ばれる。一方、このパラメーターに整数を指定すると、分離の際にランダムにその数だけ特徴量が選ばれ、その中で分離の性能がいいものが選ばれる。max_features=1とすると、ランダムに選ばれた特徴量が、その分離性能に関わらず用いられる。

以下は、make_moonsで生成された特徴量数2のデータセットについて、max_featuresを1、2と変えた時の実行結果。max_features=1の場合は、各深さにおいて特徴量1、2がランダムに選ばれる。これは必ずしも最適な分割とならないため、ノイズのように分割が滑らかになっていない。

決定木の特徴量重要度の計算方法

2020-06-05 / tau / コメントする

概要

scikit-learnのDecisionTreeClassificationモデルにfeature_importances_というパラメーターがある。このパラメーターは1次元配列で、特徴量番号に対する重要度が実数で格納されている。

このfeature_importances_について、公式ドキュメントでは以下のように書かれている。

The importance of a feature is computed as the (normalized) total reduction of the criterion brought by that feature. It is also known as the Gini importance.

～特徴量の重要度は、対象とする特徴量から得られた基準値の減少分の（正規化された）合計値。ジニ重要度としても知られている。～

と書かれているが、ちょっと曖昧で定義がはっきりしない。ジニ重要度というのは日本語サイトではなかなかヒットしないが、英語では結構取り上げられている。たとえばこちらのサイトでは以下のように引用説明されている。

It is sometimes called “gini importance” or “mean decrease impurity” and is defined as the total decrease in node impurity (weighted by the probability of reaching that node (which is approximated by the proportion of samples reaching that node)) averaged over all trees of the ensemble.

これを読むと、それはジニ重要度／平均不純度減少量と呼ばれ、ノードの不純度の減少分の重み付き和（重みはそのノードにたどり着いたサンプル数の比率）を決定木全体にわたって平均した値、となる。

定式化

あるノードの不純度をI(t_P)、その左右の子ノードの不純度をI(t_L), (t_R)とし、それぞれのノードのサンプル数をn_P, n_L, n_Rとする(n_P = n_L + n_R)。このとき、ノードt_Pの不純度の減少分の重み付き和は以下のようになる。

(1) $\begin{equation*} \Delta I(t_P) = \frac{n_p}{N} I(t_P) - \frac{n_L}{N} I(t_L) - \frac{n_R}{N} I(t_R) \end{equation*}$

ここでNは全サンプル数。この値を決定木全体にわたって平均したものが特徴量重要度となるので、これをM(t_P)とすると、以下のようになる。

(2) $\begin{equation*} M(t_P) = \frac{\Delta I(t_P)}{\displaystyle \sum_P^{all~nodes} \Delta I(t_P)} \end{equation*}$

なお分母分子でNが共通なので、式(1)においてNで割らずに計算しても結果は同じになる。

以上から、特徴量重要度の計算は以下の手順となる。

決定木の葉ノードを除く各ノードについて以下を計算
1. ノードの不純度とサンプル数を掛けた値(wI)を計算
2. ノードのwIから、左右の子ノードのwIを減じた値gを計算
決定木全体のgの合計でこれを除した値を、そのノードの分割基準となった特徴量の特徴量重要度とする

feature_importances_の内容

この流れを、breast_cancerデータセットに対して以下のコードで確認してみる。パラメーターの設定はO’Reillyの”Pythonではじめる機械学習”の例に合わせていて、深さ4で事前刈込をしている。

import matplotlib.pyplot as plt
from pandas import DataFrame
from sklearn.model_selection import train_test_split
from sklearn.datasets import load_breast_cancer
from sklearn.tree import DecisionTreeClassifier

ds = load_breast_cancer()
X_train, X_test, y_train, y_test =\
    train_test_split(ds.data, ds.target, stratify=ds.target, random_state=42)

clf = DecisionTreeClassifier(max_depth=4, random_state=0)
clf.fit(X_train, y_train)

print(clf.feature_importances_)

df = DataFrame(
    {'feature':ds.feature_names, 'importance':clf.feature_importances_})
print(df)

fig, ax = plt.subplots()
fig.subplots_adjust(left=0.3)
ax.barh(df.feature, df.importance)
plt.show()

import matplotlib.pyplot as plt

from pandas import DataFrame

from sklearn.model_selection import train_test_split

from sklearn.datasets import load_breast_cancer

from sklearn.tree import DecisionTreeClassifier

ds = load_breast_cancer()

X_train, X_test, y_train, y_test =\

train_test_split(ds.data, ds.target, stratify=ds.target, random_state=42)

clf = DecisionTreeClassifier(max_depth=4, random_state=0)

clf.fit(X_train, y_train)

print(clf.feature_importances_)

df = DataFrame(

{'feature':ds.feature_names, 'importance':clf.feature_importances_})

print(df)

fig, ax = plt.subplots()

fig.subplots_adjust(left=0.3)

ax.barh(df.feature, df.importance)

plt.show()

まず、この決定木のfeature_importances_パラメーターそのものの内容は以下の通り。すべての値を合計すると1.0となる。

[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.        ]

[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. ]

また、pandas.DataFrameで特徴量名と併せて表示すると以下の通り。

                    feature  importance
0               mean radius    0.000000
1              mean texture    0.000000
2            mean perimeter    0.000000
3                 mean area    0.000000
4           mean smoothness    0.000000
5          mean compactness    0.000000
6            mean concavity    0.000000
7       mean concave points    0.000000
8             mean symmetry    0.000000
9    mean fractal dimension    0.000000
10             radius error    0.010197
11            texture error    0.048398
12          perimeter error    0.000000
13               area error    0.000000
14         smoothness error    0.002416
15        compactness error    0.000000
16          concavity error    0.000000
17     concave points error    0.000000
18           symmetry error    0.000000
19  fractal dimension error    0.000000
20             worst radius    0.726829
21            worst texture    0.045816
22          worst perimeter    0.000000
23               worst area    0.000000
24         worst smoothness    0.014158
25        worst compactness    0.000000
26          worst concavity    0.018188
27     worst concave points    0.122113
28           worst symmetry    0.011885
29  worst fractal dimension    0.000000

feature importance

0 mean radius 0.000000

1 mean texture 0.000000

2 mean perimeter 0.000000

3 mean area 0.000000

4 mean smoothness 0.000000

5 mean compactness 0.000000

6 mean concavity 0.000000

7 mean concave points 0.000000

8 mean symmetry 0.000000

9 mean fractal dimension 0.000000

10 radius error 0.010197

11 texture error 0.048398

12 perimeter error 0.000000

13 area error 0.000000

14 smoothness error 0.002416

15 compactness error 0.000000

16 concavity error 0.000000

17 concave points error 0.000000

18 symmetry error 0.000000

19 fractal dimension error 0.000000

20 worst radius 0.726829

21 worst texture 0.045816

22 worst perimeter 0.000000

23 worst area 0.000000

24 worst smoothness 0.014158

25 worst compactness 0.000000

26 worst concavity 0.018188

27 worst concave points 0.122113

28 worst symmetry 0.011885

29 worst fractal dimension 0.000000

特徴量重要度の計算過程

特徴量重要度の計算過程を視覚的に追ってみる。まず深さ4までの決定木をgraphvizで視覚化すると以下の通り。

この決定木の葉以外のノードについて、gini不純度とサンプル数を掛け合わせた値をwIとして決定木を描きなおすと以下のとおり。また、gは着目するノードとその左右の子ノードのwIの差で、総サンプル数で無次元化しない重みによる情報利得と等価。

(3) $\begin{equation*} g = n_p I(t_P) - n_L I(t_L) - n_R I(t_R) \end{equation*}$

まずmax_depth=1のとき、worst radiusについてg = 138.6130となり、特徴量重要度はこの1つに対して1となる。

次にmax_depth=2とすると、3つのノードについて特徴量とgの値、重要度は以下の通り。

特徴量	g	重要度
worst radius	138.6130	0.809982
wors concave points	23.28798	0.136083
texture error	9.229907	0.053935
	171.13089	1

実際の計算結果は以下の通りで符合している。

                    feature  importance
0               mean radius    0.000000
1              mean texture    0.000000
2            mean perimeter    0.000000
3                 mean area    0.000000
4           mean smoothness    0.000000
5          mean compactness    0.000000
6            mean concavity    0.000000
7       mean concave points    0.000000
8             mean symmetry    0.000000
9    mean fractal dimension    0.000000
10             radius error    0.000000
11            texture error    0.053935
12          perimeter error    0.000000
13               area error    0.000000
14         smoothness error    0.000000
15        compactness error    0.000000
16          concavity error    0.000000
17     concave points error    0.000000
18           symmetry error    0.000000
19  fractal dimension error    0.000000
20             worst radius    0.809981
21            worst texture    0.000000
22          worst perimeter    0.000000
23               worst area    0.000000
24         worst smoothness    0.000000
25        worst compactness    0.000000
26          worst concavity    0.000000
27     worst concave points    0.136084
28           worst symmetry    0.000000
29  worst fractal dimension    0.000000

feature importance

0 mean radius 0.000000

1 mean texture 0.000000

2 mean perimeter 0.000000

3 mean area 0.000000

4 mean smoothness 0.000000

5 mean compactness 0.000000

6 mean concavity 0.000000

7 mean concave points 0.000000

8 mean symmetry 0.000000

9 mean fractal dimension 0.000000

10 radius error 0.000000

11 texture error 0.053935

12 perimeter error 0.000000

13 area error 0.000000

14 smoothness error 0.000000

15 compactness error 0.000000

16 concavity error 0.000000

17 concave points error 0.000000

18 symmetry error 0.000000

19 fractal dimension error 0.000000

20 worst radius 0.809981

21 worst texture 0.000000

22 worst perimeter 0.000000

23 worst area 0.000000

24 worst smoothness 0.000000

25 worst compactness 0.000000

26 worst concavity 0.000000

27 worst concave points 0.136084

28 worst symmetry 0.000000

29 worst fractal dimension 0.000000

同じように計算していき、max_depth=4の時は以下の通り。ただしここで、worst textureが2回登場していることに注意。1つ目は深さ2の左から2番目、もう1つは深さ3の右から2番目で、それぞれのノードが分割されたときのクラスは異なっている。worst textureの重要度を計算する際には、この2つのgを加えている。

特徴量	g	重要度
worst radius	138.6130	0.726829
worst concave points	23.28798	0.122113
texture error	9.229907	0.048398
radius error	1.944615	0.010197
worst texture	8.737504	0.045816
worst concavity	3.168669	0.018188
smoothness error	0.460808	0.002416
worst smoothness	2.7	0.001416
worst symmetry	2.266668	0.011885
	190.7091	1

実際の計算結果は以下の通りで符合している。

                    feature  importance
0               mean radius    0.000000
1              mean texture    0.000000
2            mean perimeter    0.000000
3                 mean area    0.000000
4           mean smoothness    0.000000
5          mean compactness    0.000000
6            mean concavity    0.000000
7       mean concave points    0.000000
8             mean symmetry    0.000000
9    mean fractal dimension    0.000000
10             radius error    0.010197
11            texture error    0.048398
12          perimeter error    0.000000
13               area error    0.000000
14         smoothness error    0.002416
15        compactness error    0.000000
16          concavity error    0.000000
17     concave points error    0.000000
18           symmetry error    0.000000
19  fractal dimension error    0.000000
20             worst radius    0.726829
21            worst texture    0.045816
22          worst perimeter    0.000000
23               worst area    0.000000
24         worst smoothness    0.014158
25        worst compactness    0.000000
26          worst concavity    0.018188
27     worst concave points    0.122113
28           worst symmetry    0.011885
29  worst fractal dimension    0.000000

feature importance

0 mean radius 0.000000

1 mean texture 0.000000

2 mean perimeter 0.000000

3 mean area 0.000000

4 mean smoothness 0.000000

5 mean compactness 0.000000

6 mean concavity 0.000000

7 mean concave points 0.000000

8 mean symmetry 0.000000

9 mean fractal dimension 0.000000

10 radius error 0.010197

11 texture error 0.048398

12 perimeter error 0.000000

13 area error 0.000000

14 smoothness error 0.002416

15 compactness error 0.000000

16 concavity error 0.000000

17 concave points error 0.000000

18 symmetry error 0.000000

19 fractal dimension error 0.000000

20 worst radius 0.726829

21 worst texture 0.045816

22 worst perimeter 0.000000

23 worst area 0.000000

24 worst smoothness 0.014158

25 worst compactness 0.000000

26 worst concavity 0.018188

27 worst concave points 0.122113

28 worst symmetry 0.011885

29 worst fractal dimension 0.000000

DecisionTreeClassifierの可視化環境

2020-05-27 / tau / コメントする

概要

Pythonのscikit-learnで提供される決定木のクラス分類モデルDecisionTreeClassifierの実行結果を可視化する環境について。

Graphvizとgraphvizパッケージ

この方法は、決定木の画像がPDFとして生成され、デフォルトのPDFリーダーが自動的に起動して確認できる。画像ファイルを利用する場合、PDFから切り出すか、以下のpydotplusパッケージを利用する。

import graphviz
from sklearn.datasets import make_moons
from sklearn.tree import DecisionTreeClassifier, export_graphviz

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

treeclf = DecisionTreeClassifier(random_state=0)
treeclf.fit(X, y)

dot_data = export_graphviz(treeclf, out_file=None,
    feature_names=["f0", "f1"], class_names=["c0", "c1"])
graph = graphviz.Source(dot_data)
graph.render("image", view=True)

import graphviz

from sklearn.datasets import make_moons

from sklearn.tree import DecisionTreeClassifier, export_graphviz

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

treeclf = DecisionTreeClassifier(random_state=0)

treeclf.fit(X, y)

dot_data = export_graphviz(treeclf, out_file=None,

feature_names=["f0", "f1"], class_names=["c0", "c1"])

graph = graphviz.Source(dot_data)

graph.render("image", view=True)

Graphvizとpydotplosパッケージ

この方法は、決定木の画像がファイルとして生成・保存される。画像を確認するためにファイルが保存されたディレクトリでファイルを開く手順が必要になるが、得られたファイルをそのまま活用することができる。

pydotplusのインストール

pydotplusをインストールする。

> pip install pydotplus

1	> pip install pydotplus

これだけでは次のようなエラーが出る。

pydotplus.graphviz.InvocationException: GraphViz's executables not found

1	pydotplus.graphviz.InvocationException: GraphViz's executables not found

Graphvizのインストール

Graphvizのサイトから実行ファイル(msiファイル)をダウンロード、インストールする。

実行方法1：Graphvizの実行位置を指定

以下のコード例13行目のように、Graphvizの実行プログラムの位置を指定。

import pydotplus as pdp
from sklearn.datasets import make_moons
from sklearn.tree import DecisionTreeClassifier, export_graphviz

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

treeclf = DecisionTreeClassifier(random_state=0)
treeclf.fit(X, y)

dot_data = export_graphviz(treeclf, out_file=None,
    feature_names=["f0", "f1"], class_names=["c0", "c1"])
graph = pdp.graph_from_dot_data(dot_data)
graph.progs = {'dot': u"C:\\Program Files (x86)\\Graphviz2.38\\bin\\dot.exe"}
graph.write_png("graph.png")

import pydotplus as pdp

from sklearn.datasets import make_moons

from sklearn.tree import DecisionTreeClassifier, export_graphviz

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

treeclf = DecisionTreeClassifier(random_state=0)

treeclf.fit(X, y)

dot_data = export_graphviz(treeclf, out_file=None,

feature_names=["f0", "f1"], class_names=["c0", "c1"])

graph = pdp.graph_from_dot_data(dot_data)

graph.progs = {'dot': u"C:\\Program Files (x86)\\Graphviz2.38\\bin\\dot.exe"}

graph.write_png("graph.png")

実行方法2：Graphvizへのパスを環境変数に登録

環境変数に上記のGraphvizのパスを指定する。

デスクトップのPCアイコンを右クリック→プロパティ
システム・ウィンドウ→システムの詳細設定
システムのプロパティダイアログ→環境変数ボタン
環境変数ダイアログのシステム環境変数→Pathを指定して編集ボタン
環境変数名の編集ダイアログ→新規ボタン
Graphvizへのパス(例えばC:\Program Files (x86)\Graphviz2.38\bin\)を入力してOK
以下、各ダイアログでOK

環境変数を設定しておくと、毎回パスを指定しなくてよい。

import pydotplus as pdp
from sklearn.datasets import make_moons
from sklearn.tree import DecisionTreeClassifier, export_graphviz

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

treeclf = DecisionTreeClassifier(random_state=0)
treeclf.fit(X, y)

dot_data = export_graphviz(treeclf, out_file=None,
    feature_names=["f0", "f1"], class_names=["c0", "c1"])
graph = pdp.graph_from_dot_data(dot_data)
graph.write_png("graph.png")

import pydotplus as pdp

from sklearn.datasets import make_moons

from sklearn.tree import DecisionTreeClassifier, export_graphviz

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

treeclf = DecisionTreeClassifier(random_state=0)

treeclf.fit(X, y)

dot_data = export_graphviz(treeclf, out_file=None,

feature_names=["f0", "f1"], class_names=["c0", "c1"])

graph = pdp.graph_from_dot_data(dot_data)

graph.write_png("graph.png")

dtreeviz

dtreevizのインストール

dtreevizをインストールする。

> pip install dtreeviz

1	> pip install dtreeviz

実行方法

Graphvizの実行方法2で環境変数を追加。

from sklearn.datasets import make_moons
from sklearn.tree import DecisionTreeClassifier
from dtreeviz.trees import dtreeviz

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

treeclf = DecisionTreeClassifier(random_state=0)
treeclf.fit(X, y)

viz = dtreeviz(treeclf, X, y, target_name="Classes",
    feature_names=["f0", "f1"], class_names=["c0", "c1"])

viz.view()

from sklearn.datasets import make_moons

from sklearn.tree import DecisionTreeClassifier

from dtreeviz.trees import dtreeviz

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

treeclf = DecisionTreeClassifier(random_state=0)

treeclf.fit(X, y)

viz = dtreeviz(treeclf, X, y, target_name="Classes",

feature_names=["f0", "f1"], class_names=["c0", "c1"])

viz.view()