機械学習ことはじめ（Vol.75掲載）

著者：川嶋宏彰

本連載では、機械学習の基礎となるさまざまな手法の仕組みや、それらの手法のPythonでの利用方法を解説していきます。今回は、入力データがどのようなクラス（カテゴリ）であるかを予測する分類の問題を扱います。

シェルスクリプトマガジン Vol.75は以下のリンク先でご購入できます。

図2　データセットの読み込みから、散布図のプロットまでを実施するPythonコード

import seaborn as sns
import matplotlib.pyplot as plt

plt.rcParams['font.size'] = 14
penguins = sns.load_dataset('penguins')
# 取り出す特徴量
features = ['bill_depth_mm', 'body_mass_g']  # （★）
# features = ['bill_depth_mm', 'bill_length_mm']  # 取り出す特徴量を変える
# 対象とするペンギン種
target_species = ['Adelie', 'Gentoo']  # （★）
# target_species = ['Adelie', 'Chinstrap', 'Gentoo']  # 3種のペンギンにする
# 今回用いる特徴量をクラスラベルと共に取り出す
df = penguins[['species'] + features].copy()
df.dropna(inplace=True)  # NaN が含まれる行は削除
# 今回用いるペンギン種のみとする
df2 = df[df['species'].isin(target_species)].copy()
print(df2.shape)  # (274, 3) と表示される
# 現在のパレットから用いる種 (target_species) に合わせたパレットを作成
palette = {c: sns.color_palette()[k] for k, c
           in enumerate(df['species'].unique()) if c in target_species}
fig = plt.figure(figsize=(5, 5))
sns.scatterplot(data=df2, x=features[0], y=features[1],
                hue='species', palette=palette)
plt.show()

import seaborn as sns

import matplotlib.pyplot as plt

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

penguins = sns.load_dataset('penguins')

# 取り出す特徴量

features = ['bill_depth_mm', 'body_mass_g'] # （★）

# features = ['bill_depth_mm', 'bill_length_mm'] # 取り出す特徴量を変える

# 対象とするペンギン種

target_species = ['Adelie', 'Gentoo'] # （★）

# target_species = ['Adelie', 'Chinstrap', 'Gentoo'] # 3種のペンギンにする

# 今回用いる特徴量をクラスラベルと共に取り出す

df = penguins[['species'] + features].copy()

df.dropna(inplace=True) # NaN が含まれる行は削除

# 今回用いるペンギン種のみとする

df2 = df[df['species'].isin(target_species)].copy()

print(df2.shape) # (274, 3) と表示される

# 現在のパレットから用いる種 (target_species) に合わせたパレットを作成

palette = {c: sns.color_palette()[k] for k, c

in enumerate(df['species'].unique()) if c in target_species}

fig = plt.figure(figsize=(5, 5))

sns.scatterplot(data=df2, x=features[0], y=features[1],

hue='species', palette=palette)

plt.show()

図3　ロジスティック回帰の準備をするPythonコード

from sklearn.linear_model import LogisticRegression
import numpy as np

X = df2[features].values   # 各個体の2次元特徴量（行数 = 個体数）
y = df2['species'].values  # 各個体のクラスラベル
clf_lr = LogisticRegression(penalty='none')  # 分類モデルの準備

from sklearn.linear_model import LogisticRegression

import numpy as np

X = df2[features].values # 各個体の2次元特徴量（行数 = 個体数）

y = df2['species'].values # 各個体のクラスラベル

clf_lr = LogisticRegression(penalty='none') # 分類モデルの準備

図4　特徴量一つ（体重）のみでロジスティック回帰をするPythonコード

f_idx = 1  # 用いる特徴量のインデックス
# 指定した列（特徴量）だけ取り出して学習
clf_lr.fit(X[:, [f_idx]], y)
# 学習されたパラメータを表示
w0 = clf_lr.intercept_[0]
w1 = clf_lr.coef_[0][0]
print(f'w_0 = {w0:.2f}, w_1 = {w1:.4f}, b = -w0/w1 = {-w0/w1:.2f}')
# 予測値および分類結果の可視化用グリッド
xmin = np.min(X[:, f_idx])
xmax = np.max(X[:, f_idx])
xs = np.linspace(xmin - (xmax - xmin) / 10,
                 xmax + (xmax - xmin) / 10, 100)  # 等間隔に100点用意
fig = plt.figure(figsize=(10, 5))
# 元データのプロット
sns.scatterplot(x=X[:, f_idx], y=(y=='Gentoo'), hue=y, palette=palette)
# 予測値（事後確率）の曲線をプロット
y_pred = clf_lr.predict_proba(xs.reshape(-1, 1))[:, 1]
plt.plot(xs, y_pred, c=palette['Gentoo'])
plt.ylim(-0.1, 1.1)
plt.xlabel(f'x ({features[f_idx]})')
plt.ylabel('y')
plt.grid()
plt.show()

f_idx = 1 # 用いる特徴量のインデックス

# 指定した列（特徴量）だけ取り出して学習

clf_lr.fit(X[:, [f_idx]], y)

# 学習されたパラメータを表示

w0 = clf_lr.intercept_[0]

w1 = clf_lr.coef_[0][0]

print(f'w_0 = {w0:.2f}, w_1 = {w1:.4f}, b = -w0/w1 = {-w0/w1:.2f}')

# 予測値および分類結果の可視化用グリッド

xmin = np.min(X[:, f_idx])

xmax = np.max(X[:, f_idx])

xs = np.linspace(xmin - (xmax - xmin) / 10,

xmax + (xmax - xmin) / 10, 100) # 等間隔に100点用意

fig = plt.figure(figsize=(10, 5))

# 元データのプロット

sns.scatterplot(x=X[:, f_idx], y=(y=='Gentoo'), hue=y, palette=palette)

# 予測値（事後確率）の曲線をプロット

y_pred = clf_lr.predict_proba(xs.reshape(-1, 1))[:, 1]

plt.plot(xs, y_pred, c=palette['Gentoo'])

plt.ylim(-0.1, 1.1)

plt.xlabel(f'x ({features[f_idx]})')

plt.ylabel('y')

plt.grid()

plt.show()

図6　特徴量二つ（くちばしの高さと体重）でロジスティック回帰をするPythonコード

import matplotlib as mpl
from sklearn.preprocessing import StandardScaler

def plot_decision_boundary(X, y, clf, xylabels=None, palette=None, fig=None, ngrid=50):
  """ 分類器 clf の決定境界を描画 """
  if fig is None: fig = plt.figure(figsize=(5, 5))
  else: plt.figure(fig.number)
  # 2次元空間にグリッド点を準備
  xmin = X.min(axis=0)  # 各列の最小値
  xmax = X.max(axis=0)  # 各列の最大値
  xstep = [(xmax[j]-xmin[j]) / ngrid for j in range(2)]  # グリッドのステップ幅
  xmin = [xmin[j] - 8*xstep[j] for j in range(2)]  # 少し広めに
  xmax = [xmax[j] + 8*xstep[j] for j in range(2)]
  aranges = [np.arange(xmin[j], xmax[j] + xstep[j], xstep[j]) for j in range(2)]
  x0grid, x1grid = np.meshgrid(*aranges)
  # 各グリッド点でクラスを判定
  y_pred = clf.predict(np.c_[x0grid.ravel(), x1grid.ravel()])
  y_pred = y_pred.reshape(x0grid.shape)  # 2次元に
  y_pred = np.searchsorted(np.unique(y_pred), y_pred)  # 値をインデックスに
  clist = palette.values() if type(palette) is dict else palette
  cmap = mpl.colors.ListedColormap(clist) if palette is not None else None
  plt.contourf(x0grid, x1grid, y_pred, alpha=0.3, cmap=cmap)
  sns.scatterplot(x=X[:, 0], y=X[:, 1], hue=y, palette=palette)
  plt.legend()
  plt.xlim([xmin[0], xmax[0]])
  plt.ylim([xmin[1], xmax[1]])
  if xylabels is not None:
    plt.xlabel(xylabels[0])
    plt.ylabel(xylabels[1])
  return fig

Xs = StandardScaler().fit_transform(X)  # 標準化によるスケーリング
clf_lr.fit(Xs, y)  # 二つの特徴量を両方用いて学習
print('w0:', clf_lr.intercept_)
print('w1, w2:', clf_lr.coef_)
xylabels = [s + ' (standardized)' for s in features]  # 軸ラベル変更
plot_decision_boundary(Xs, y, clf_lr, xylabels, palette, ngrid=100)  # 決定境界
plt.show()

import matplotlib as mpl

from sklearn.preprocessing import StandardScaler

def plot_decision_boundary(X, y, clf, xylabels=None, palette=None, fig=None, ngrid=50):

""" 分類器 clf の決定境界を描画 """

if fig is None: fig = plt.figure(figsize=(5, 5))

else: plt.figure(fig.number)

# 2次元空間にグリッド点を準備

xmin = X.min(axis=0) # 各列の最小値

xmax = X.max(axis=0) # 各列の最大値

xstep = [(xmax[j]-xmin[j]) / ngrid for j in range(2)] # グリッドのステップ幅

xmin = [xmin[j] - 8*xstep[j] for j in range(2)] # 少し広めに

xmax = [xmax[j] + 8*xstep[j] for j in range(2)]

aranges = [np.arange(xmin[j], xmax[j] + xstep[j], xstep[j]) for j in range(2)]

x0grid, x1grid = np.meshgrid(*aranges)

# 各グリッド点でクラスを判定

y_pred = clf.predict(np.c_[x0grid.ravel(), x1grid.ravel()])

y_pred = y_pred.reshape(x0grid.shape) # 2次元に

y_pred = np.searchsorted(np.unique(y_pred), y_pred) # 値をインデックスに

clist = palette.values() if type(palette) is dict else palette

cmap = mpl.colors.ListedColormap(clist) if palette is not None else None

plt.contourf(x0grid, x1grid, y_pred, alpha=0.3, cmap=cmap)

sns.scatterplot(x=X[:, 0], y=X[:, 1], hue=y, palette=palette)

plt.legend()

plt.xlim([xmin[0], xmax[0]])

plt.ylim([xmin[1], xmax[1]])

if xylabels is not None:

plt.xlabel(xylabels[0])

plt.ylabel(xylabels[1])

return fig

Xs = StandardScaler().fit_transform(X) # 標準化によるスケーリング

clf_lr.fit(Xs, y) # 二つの特徴量を両方用いて学習

print('w0:', clf_lr.intercept_)

print('w1, w2:', clf_lr.coef_)

xylabels = [s + ' (standardized)' for s in features] # 軸ラベル変更

plot_decision_boundary(Xs, y, clf_lr, xylabels, palette, ngrid=100) # 決定境界

plt.show()

図9　ロジスティック回帰による各点での予測値をプロットするPythonコード

from matplotlib import cm
from mpl_toolkits.mplot3d import axes3d

# メッシュで
def gen_2dgrid(X):
  """ 2次元メッシュグリッドの生成 """
  d = X.shape[1]
  xmin = X.min(axis=0)  # 各列の最小値
  xmax = X.max(axis=0)  # 各列の最大値
  xstep = [(xmax[j]-xmin[j]) / 100 for j in range(d)]  # グリッドのステップ幅
  xmin = [xmin[j] - 10*xstep[j] for j in range(d)]  # 少し広めに
  xmax = [xmax[j] + 10*xstep[j] for j in range(d)]
  aranges = [np.arange(xmin[j], xmax[j] + xstep[j], xstep[j]) for j in range(2)]
  return np.meshgrid(*aranges)  # d=2のときは (x0grid, x1grid) が返る
# 2次元の場合の予測値 (proba)
def proba_2dgrid(clf, x0grid, x1grid):
  """ 2次元メッシュの各点での予測値を計算 """    
  return clf.predict_proba(np.c_[x0grid.ravel(), x1grid.ravel()])[:, 1]
                          .reshape(x0grid.shape)
xgrid_2d = gen_2dgrid(Xs)  # xgrid_2d: (x0grid, x1grid) のような二つ組
proba_2d = proba_2dgrid(clf_lr, *xgrid_2d)  # 予測値
# 3次元プロット
fig = plt.figure(figsize=(14, 8))
ax = fig.add_subplot(111, projection='3d')
cmap = cm.winter  # カラーマップを設定
ax.plot_surface(*xgrid_2d, proba_2d, cmap=cmap)
# ax.set_box_aspect((3, 3, 1))  # matplotlib のバージョンが3.3以上で利用可能
ax.set_zlim(-0.1, 1.1)
ax.set_xlabel(xylabels[0])
ax.set_ylabel(xylabels[1])
ax.set_zlabel('y')
ax.set_zticks([0, 0.5, 1.0])
ax.view_init(elev=60, azim=-120)  # 3次元プロットの表示角度の設定
plt.show()

from matplotlib import cm

from mpl_toolkits.mplot3d import axes3d

# メッシュで

def gen_2dgrid(X):

""" 2次元メッシュグリッドの生成 """

d = X.shape[1]

xmin = X.min(axis=0) # 各列の最小値

xmax = X.max(axis=0) # 各列の最大値

xstep = [(xmax[j]-xmin[j]) / 100 for j in range(d)] # グリッドのステップ幅

xmin = [xmin[j] - 10*xstep[j] for j in range(d)] # 少し広めに

xmax = [xmax[j] + 10*xstep[j] for j in range(d)]

aranges = [np.arange(xmin[j], xmax[j] + xstep[j], xstep[j]) for j in range(2)]

return np.meshgrid(*aranges) # d=2のときは (x0grid, x1grid) が返る

# 2次元の場合の予測値 (proba)

def proba_2dgrid(clf, x0grid, x1grid):

""" 2次元メッシュの各点での予測値を計算 """

return clf.predict_proba(np.c_[x0grid.ravel(), x1grid.ravel()])[:, 1]

.reshape(x0grid.shape)

xgrid_2d = gen_2dgrid(Xs) # xgrid_2d: (x0grid, x1grid) のような二つ組

proba_2d = proba_2dgrid(clf_lr, *xgrid_2d) # 予測値

# 3次元プロット

fig = plt.figure(figsize=(14, 8))

ax = fig.add_subplot(111, projection='3d')

cmap = cm.winter # カラーマップを設定

ax.plot_surface(*xgrid_2d, proba_2d, cmap=cmap)

# ax.set_box_aspect((3, 3, 1)) # matplotlib のバージョンが3.3以上で利用可能

ax.set_zlim(-0.1, 1.1)

ax.set_xlabel(xylabels[0])

ax.set_ylabel(xylabels[1])

ax.set_zlabel('y')

ax.set_zticks([0, 0.5, 1.0])

ax.view_init(elev=60, azim=-120) # 3次元プロットの表示角度の設定

plt.show()

図12　線形SVM による学習と決定境界の表示をするPythonコード

from sklearn.svm import SVC

C = 10  # 大きいほどハードマージンに近い
clf_lin = SVC(kernel='linear', C=C)
clf_lin.fit(Xs, y)
# 決定境界とサポートベクターを可視化
plot_decision_boundary(Xs, y, clf_lin, xylabels, palette, ngrid=100)
plt.scatter(clf_lin.support_vectors_[:, 0], clf_lin.support_vectors_[:, 1],
            s=50, facecolors='none', edgecolors='r')
plt.title(f'SVM(linear) C = {C}')
plt.show()

from sklearn.svm import SVC

C = 10 # 大きいほどハードマージンに近い

clf_lin = SVC(kernel='linear', C=C)

clf_lin.fit(Xs, y)

# 決定境界とサポートベクターを可視化

plot_decision_boundary(Xs, y, clf_lin, xylabels, palette, ngrid=100)

plt.scatter(clf_lin.support_vectors_[:, 0], clf_lin.support_vectors_[:, 1],

s=50, facecolors='none', edgecolors='r')

plt.title(f'SVM(linear) C = {C}')

plt.show()

図17　非線形SVMとパラメータCの関係を示すPythonコード

from sklearn.datasets import make_moons

X_moon, y_moon = make_moons(n_samples=100, noise=0.2, random_state=6)
palette_o = {k: sns.color_palette()[k] for k in range(2)}
# moonデータの散布図をプロット
plt.figure(figsize=(5, 5))
sns.scatterplot(x=X_moon[:, 0], y=X_moon[:, 1],
                hue=y_moon, palette=palette_o)
plt.show()
# 異なるCでSVM（RBFカーネル）の学習と決定境界の表示
for C in [0.1, 10, 10000]:
  clf_rbf = SVC(C=C, kernel='rbf')  # kernel='rbf'は省略可能
  clf_rbf.fit(X_moon, y_moon)
  # 決定境界をプロット
  plot_decision_boundary(X_moon, y_moon,
                         clf_rbf, None, palette_o, ngrid=100)
  plt.title(f'SVM(rbf) C = {C}')
  plt.show()

from sklearn.datasets import make_moons

X_moon, y_moon = make_moons(n_samples=100, noise=0.2, random_state=6)

palette_o = {k: sns.color_palette()[k] for k in range(2)}

# moonデータの散布図をプロット

plt.figure(figsize=(5, 5))

sns.scatterplot(x=X_moon[:, 0], y=X_moon[:, 1],

hue=y_moon, palette=palette_o)

plt.show()

# 異なるCでSVM（RBFカーネル）の学習と決定境界の表示

for C in [0.1, 10, 10000]:

clf_rbf = SVC(C=C, kernel='rbf') # kernel='rbf'は省略可能

clf_rbf.fit(X_moon, y_moon)

# 決定境界をプロット

plot_decision_boundary(X_moon, y_moon,

clf_rbf, None, palette_o, ngrid=100)

plt.title(f'SVM(rbf) C = {C}')

plt.show()

図20　グリッドサーチをするPythonコード

from sklearn.model_selection import cross_val_score

best_score = -1
for gamma in [0.1, 0.5, 1, 2, 10]:  # 5通り
  for C in [0.1, 1, 10, 100, 1000, 10000]:  # 6通り
    print(f'gamma: {gamma},\tC: {C}', end='')
    svm = SVC(gamma=gamma, C=C)
    cv_scores = cross_val_score(svm, X_moon, y_moon, cv=5)
    score = np.mean(cv_scores)
    print(f'\t | average: {score:.2} <- {cv_scores}')
    if score > best_score:
      best_score = score
      best_params = {'C': C, 'gamma': gamma}
print('best_score:', best_score)
print('best_params:', best_params)
clf_rbf = SVC(**best_params)
clf_rbf.fit(X_moon, y_moon)
plot_decision_boundary(X_moon, y_moon,
                       clf_rbf, None, palette_o, ngrid=100)
title = f"C = {best_params['C']}, gamma = {best_params['gamma']}"
plt.title(title)
plt.show()

from sklearn.model_selection import cross_val_score

best_score = -1

for gamma in [0.1, 0.5, 1, 2, 10]: # 5通り

for C in [0.1, 1, 10, 100, 1000, 10000]: # 6通り

print(f'gamma: {gamma},\tC: {C}', end='')

svm = SVC(gamma=gamma, C=C)

cv_scores = cross_val_score(svm, X_moon, y_moon, cv=5)

score = np.mean(cv_scores)

print(f'\t | average: {score:.2} <- {cv_scores}')

if score > best_score:

best_score = score

best_params = {'C': C, 'gamma': gamma}

print('best_score:', best_score)

print('best_params:', best_params)

clf_rbf = SVC(**best_params)

clf_rbf.fit(X_moon, y_moon)

plot_decision_boundary(X_moon, y_moon,

clf_rbf, None, palette_o, ngrid=100)

title = f"C = {best_params['C']}, gamma = {best_params['gamma']}"

plt.title(title)

plt.show()