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

筆者：川嶋宏彰

本連載では、機械学習の基礎となるさまざまな手法の仕組みや、それらの手法のPythonでの利用方法を解説していきます。今回は「確率モデル」による機械学習である、ガウス分布を用いた教師あり学習と教師なし学習の手法を紹介します。

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

図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'] # ★
# 対象とするペンギン種
target_species = ['Adelie', 'Gentoo'] # ★
# 今回用いる特徴量をクラスラベルと共に取り出す
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'] # ★

# 対象とするペンギン種

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

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

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()

図4　各クラス（ペンギン種）の体重分布をプロットするPythonコード

import numpy as np
import scipy

X = df2[features].values  # 各個体の2次元特徴量（行数=個体数）
y = df2['species'].values  # 各個体のクラスラベル
prob_y = {c: np.sum(y==c)/y.size for c in target_species}  # 事前確率
print('Prior:', prob_y)
# 平均と分散は2次元にまとめてベクトルや行列として計算しておく
mu = {c: X[y==c, :].mean(axis=0) for c in target_species}  # 平均
sigma2 = {c: X[y==c, :].var(axis=0, ddof=1) for c in target_species}  # 不偏分散
f_idx = 1  # 用いる特徴量のindex（これは1次元の場合）
fig = plt.figure(figsize=(10, 5))
ax1 = fig.add_subplot(111)
ax2 = ax1.twinx()  # 右の縦軸
# ヒストグラムをプロット (体重のみ)
sns.histplot(ax=ax1, x=X[:, f_idx], hue=y, palette=palette, alpha=0.2)
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点用意
# 確率密度関数を準備 (体重のみ)
norm_dist1d = {c: scipy.stats.multivariate_normal(mu[c][f_idx], sigma2[c][f_idx])
               for c in target_species}
for c in target_species:
    # 各クラスの確率密度関数をプロット
    sns.lineplot(ax=ax2, x=xs, y=norm_dist1d[c].pdf(xs)*prob_y[c],
                 hue=[c]*len(xs), palette=palette)
# 各データを小さな縦線でプロット
sns.rugplot(x=X[:, f_idx], hue=y, palette=palette, height=0.05)
ax2.set_ylabel('Probability density')
ax1.set_xlabel(features[f_idx])
plt.show()

import numpy as np

import scipy

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

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

prob_y = {c: np.sum(y==c)/y.size for c in target_species} # 事前確率

print('Prior:', prob_y)

# 平均と分散は2次元にまとめてベクトルや行列として計算しておく

mu = {c: X[y==c, :].mean(axis=0) for c in target_species} # 平均

sigma2 = {c: X[y==c, :].var(axis=0, ddof=1) for c in target_species} # 不偏分散

f_idx = 1 # 用いる特徴量のindex（これは1次元の場合）

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

ax1 = fig.add_subplot(111)

ax2 = ax1.twinx() # 右の縦軸

# ヒストグラムをプロット (体重のみ)

sns.histplot(ax=ax1, x=X[:, f_idx], hue=y, palette=palette, alpha=0.2)

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点用意

# 確率密度関数を準備 (体重のみ)

norm_dist1d = {c: scipy.stats.multivariate_normal(mu[c][f_idx], sigma2[c][f_idx])

for c in target_species}

for c in target_species:

# 各クラスの確率密度関数をプロット

sns.lineplot(ax=ax2, x=xs, y=norm_dist1d[c].pdf(xs)*prob_y[c],

hue=[c]*len(xs), palette=palette)

# 各データを小さな縦線でプロット

sns.rugplot(x=X[:, f_idx], hue=y, palette=palette, height=0.05)

ax2.set_ylabel('Probability density')

ax1.set_xlabel(features[f_idx])

plt.show()

図6　1次元での決定境界を求めるPythonコード

# 曲線の上下が変化するおおよその点を求める
diff = norm_dist1d['Adelie'].pdf(xs)*prob_y['Adelie'] -
       norm_dist1d['Gentoo'].pdf(xs)*prob_y['Gentoo']
# 符号の変化点を見つけるために隣り合う要素の積を計算
ddiff = diff[1:] * diff[:-1]
print('boundary:', xs[np.where(ddiff < 0)[0][0]])

# 曲線の上下が変化するおおよその点を求める

diff = norm_dist1d['Adelie'].pdf(xs)*prob_y['Adelie'] -

norm_dist1d['Gentoo'].pdf(xs)*prob_y['Gentoo']

# 符号の変化点を見つけるために隣り合う要素の積を計算

ddiff = diff[1:] * diff[:-1]

print('boundary:', xs[np.where(ddiff < 0)[0][0]])

図8　2次元ガウス分布とその等高線を表示するPythonコード

from mpl_toolkits.mplot3d import axes3d
from matplotlib import cm
# 共分散行列を求めておく
Sigma = {c: np.cov(X[y==c, :].T, ddof=1) for c in target_species}
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) が返る
def gaussian_2dgrid(m, S, x0grid, x1grid):
    """ 2次元ガウス分布の密度関数の値を2次元メッシュで計算 """
    gaussian = scipy.stats.multivariate_normal(m, S)
    return gaussian.pdf(np.c_[x0grid.ravel(), x1grid.ravel()]).reshape(x0grid.shape)
c = 'Adelie'  # プロットするクラスを設定
xgrid_2d = gen_2dgrid(X)  # xgrid_2d: (x0grid, x1grid) のような二つ組
px_adelie = gaussian_2dgrid(mu[c], Sigma[c], *xgrid_2d)  # *xgrid_2d で2引数に展開
fig = plt.figure(figsize=(8, 5))  # 3次元プロット
ax = fig.add_subplot(111, projection='3d')
# 2次元ガウシアンの3次元的プロット
cmap = cm.coolwarm  # カラーマップを設定
ax.plot_surface(*xgrid_2d, px_adelie, cmap=cmap)
# 等高線のプロット
z_offset = -np.max(px_adelie)
ax.contour(*xgrid_2d, px_adelie, zdir='z', offset=z_offset, cmap=cmap)
ax.set_zlim(z_offset, ax.get_zlim()[1])
ax.view_init(elev=40, azim=-60)  # 3次元プロットの表示角度の設定
plt.show()

from mpl_toolkits.mplot3d import axes3d

from matplotlib import cm

# 共分散行列を求めておく

Sigma = {c: np.cov(X[y==c, :].T, ddof=1) for c in target_species}

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) が返る

def gaussian_2dgrid(m, S, x0grid, x1grid):

""" 2次元ガウス分布の密度関数の値を2次元メッシュで計算 """

gaussian = scipy.stats.multivariate_normal(m, S)

return gaussian.pdf(np.c_[x0grid.ravel(), x1grid.ravel()]).reshape(x0grid.shape)

c = 'Adelie' # プロットするクラスを設定

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

px_adelie = gaussian_2dgrid(mu[c], Sigma[c], *xgrid_2d) # *xgrid_2d で2引数に展開

fig = plt.figure(figsize=(8, 5)) # 3次元プロット

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

# 2次元ガウシアンの3次元的プロット

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

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

# 等高線のプロット

z_offset = -np.max(px_adelie)

ax.contour(*xgrid_2d, px_adelie, zdir='z', offset=z_offset, cmap=cmap)

ax.set_zlim(z_offset, ax.get_zlim()[1])

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

plt.show()

図10　2次元ガウス分布の等高線と元データの散布図を表示するPythonコード

def plot_ellipses(ms, Ss, pys, xgrids, xylabels=None, fig=None):
    """ 各クラスの確率だ円をプロット """
    if fig is None: fig = plt.figure(figsize=(5, 5))  # 新たに作成
    else: plt.figure(fig.number)  # 既存のfigureを利用
    levels = None
    # クラス名を辞書キーから取得
    for c in ms.keys() if type(ms) is dict else range(len(ms)):
        cs = plt.contour(*xgrids, gaussian_2dgrid(ms[c],
                         Ss[c], *xgrids)*pys[c], cmap=cmap,
                         levels=levels)
        levels = cs.levels  # 二つ目以降は一つ目のlevelsを利用
        # 密度（山の標高）をだ円に表示
        # plt.clabel(cs, inline=True, fontsize=8)
        if xylabels is not None:
        plt.xlabel(xylabels[0])
        plt.ylabel(xylabels[1])
    return fig
plot_ellipses(mu, Sigma, prob_y, xgrid_2d, xylabels=features)
sns.scatterplot(x=X[:, 0], y=X[:, 1], hue=y, palette=palette)
plt.show()

def plot_ellipses(ms, Ss, pys, xgrids, xylabels=None, fig=None):

""" 各クラスの確率だ円をプロット """

if fig is None: fig = plt.figure(figsize=(5, 5)) # 新たに作成

else: plt.figure(fig.number) # 既存のfigureを利用

levels = None

# クラス名を辞書キーから取得

for c in ms.keys() if type(ms) is dict else range(len(ms)):

cs = plt.contour(*xgrids, gaussian_2dgrid(ms[c],

Ss[c], *xgrids)*pys[c], cmap=cmap,

levels=levels)

levels = cs.levels # 二つ目以降は一つ目のlevelsを利用

# 密度（山の標高）をだ円に表示

# plt.clabel(cs, inline=True, fontsize=8)

if xylabels is not None:

plt.xlabel(xylabels[0])

plt.ylabel(xylabels[1])

return fig

plot_ellipses(mu, Sigma, prob_y, xgrid_2d, xylabels=features)

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

plt.show()

図12　2次の識別による決定境界を表示するPythonコード

from sklearn.discriminant_analysis import QuadraticDiscriminantAnalysis
# plot_decision_boundary()は
# 注釈のリンク先を参照して別途定義
clf_qda = QuadraticDiscriminantAnalysis()
clf_qda.fit(X, y)  # 学習
# 決定境界の表示
fig = plot_decision_boundary(X, y, clf_qda,
                             features, palette)
# 確率だ円の表示
plot_ellipses(mu, Sigma, prob_y, xgrid_2d, fig=fig)
plt.show()

from sklearn.discriminant_analysis import QuadraticDiscriminantAnalysis

# plot_decision_boundary()は

# 注釈のリンク先を参照して別途定義

clf_qda = QuadraticDiscriminantAnalysis()

clf_qda.fit(X, y) # 学習

# 決定境界の表示

fig = plot_decision_boundary(X, y, clf_qda,

features, palette)

# 確率だ円の表示

plot_ellipses(mu, Sigma, prob_y, xgrid_2d, fig=fig)

plt.show()

図15　ガウス分布によるデータの生成をするPythonコード

fig = plt.figure(figsize=(5, 5))
Ngen_each = 120  # 各クラスで120個体生成する場合
# 各クラスの割合を変化させてもよい
X_gen = np.vstack(
            [np.random.multivariate_normal(mu[c], Sigma[c], Ngen_each)
            for c in target_species])
y_gen = np.hstack([[c] * Ngen_each for c in target_species])
sns.scatterplot(x=X_gen[:, 0], y=X_gen[:, 1], hue=y_gen, palette=palette)
plt.xlabel(features[0])
plt.ylabel(features[1])
plt.show()

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

Ngen_each = 120 # 各クラスで120個体生成する場合

# 各クラスの割合を変化させてもよい

X_gen = np.vstack(

[np.random.multivariate_normal(mu[c], Sigma[c], Ngen_each)

for c in target_species])

y_gen = np.hstack([[c] * Ngen_each for c in target_species])

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

plt.xlabel(features[0])

plt.ylabel(features[1])

plt.show()

図18　混合ガウス分布によるソフトクラスタリングをするPythonコード

from sklearn.mixture import GaussianMixture
from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler
k = 2  # クラスタ数（非階層的クラスタリングなので決める必要がある）
use_gmm = True  # Falseにするとk-meansになる
if use_gmm:
    gmm = GaussianMixture(n_components=k, covariance_type='full')
    # full、diagはスケーリング無しも可（sphericalでは必要）
    Xu = X
    # GMM推定による各データのクラスタID
    y_pred = gmm.fit_predict(Xu)
else:
    kmeans = KMeans(n_clusters=k)
    # 標準化によるスケーリング
    Xu = StandardScaler().fit_transform(X)
    # k-meansによる各データのクラスタID
    y_pred = kmeans.fit_predict(Xu)
fig = plt.figure(figsize=(11, 5))
fig.add_subplot(121)
sns.scatterplot(x=X[:, 0], y=X[:, 1], color='k')
plt.xlabel(features[0])
plt.ylabel(features[1])
plt.title('Unlabeled data')
fig.add_subplot(122)
sns.scatterplot(x=Xu[:, 0], y=Xu[:, 1],
                hue=y_pred, palette='bright')
method_str = 'GMM' if use_gmm else 'k-means'
plt.title(f'{method_str} clustering')
if use_gmm:
    plot_ellipses(gmm.means_, gmm.covariances_,
                  gmm.weights_, xgrid_2d, fig=fig)
plt.show()

from sklearn.mixture import GaussianMixture

from sklearn.cluster import KMeans

from sklearn.preprocessing import StandardScaler

k = 2 # クラスタ数（非階層的クラスタリングなので決める必要がある）

use_gmm = True # Falseにするとk-meansになる

if use_gmm:

gmm = GaussianMixture(n_components=k, covariance_type='full')

# full、diagはスケーリング無しも可（sphericalでは必要）

Xu = X

# GMM推定による各データのクラスタID

y_pred = gmm.fit_predict(Xu)

else:

kmeans = KMeans(n_clusters=k)

# 標準化によるスケーリング

Xu = StandardScaler().fit_transform(X)

# k-meansによる各データのクラスタID

y_pred = kmeans.fit_predict(Xu)

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

fig.add_subplot(121)

sns.scatterplot(x=X[:, 0], y=X[:, 1], color='k')

plt.xlabel(features[0])

plt.ylabel(features[1])

plt.title('Unlabeled data')

fig.add_subplot(122)

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

hue=y_pred, palette='bright')

method_str = 'GMM' if use_gmm else 'k-means'

plt.title(f'{method_str} clustering')

if use_gmm:

plot_ellipses(gmm.means_, gmm.covariances_,

gmm.weights_, xgrid_2d, fig=fig)

plt.show()