【E資格】機械学習

 ここでは、機械学習(マシーンラーニング;ML)について整理する。目的は以下の2点である。

①機械学習の基本的な手法を理解し実装する
②機械学習モデリングの流れを理解

機械学習ではアプリケーションにさせたいことを性能指標に関して経験から学習させる。その際、人間がプログラムするのは認識の仕方ではなく学習の仕方である。

1. 線形回帰モデル

線形回帰

 回帰問題とは、ある入力から出力を予測する問題である。直線で予測する場合を線形回帰、曲線で予測する場合を非線形回帰という。この際の入力の各要素を説明変数/特微量(m次元のベクトル)といい、出力を目的変数(スカラー値)という。

$$ \begin{eqnarray}説明変数&:&{\bf x} = (x_1, x_2, \cdot\cdot\cdot, x_m)^T \in \mathbb{R}^m \\
目的変数&:&y \in \mathbb{R}^1 \end{eqnarray}$$

 回帰問題を解くための機械学習モデルの一つとして線形回帰モデルがある。線形回帰モデルでは入力とm次元パラメータを線型結合し、出力する。線型結合は「入力とパラメータの内積+バイアス」を行う処理を指す。

$$ \begin{eqnarray}パラメータ&:&{\bf w} = (w_1, w_2, \cdot\cdot\cdot, w_m)^T \in \mathbb{R}^m \\
線型結合&:&\hat{y} = {\bf w}^T{\bf x} + w_0 = \sum^m_{j=1}w_jx_j + w_0 \end{eqnarray}

ここでパラメータとは、モデルに含まれる未知の変数を指す。学習における重みを指し、機械学習や深層学習において重要な指標である。

データ分割

 学習を行う上で必要なことは、「限りのあるデータを有効に活用して、現実の未知のデータに対する学習制度を高める」ことである。このためにはデータの分割が行われる。データの分割では全てのデータを学習用データと検証用データに分割する。

  • 学習用データ:機械学習モデルの学習に利用するデータ
  • 検証用データ:学習済みモデルの制度を検証するためのデータ

まずは学習データを用いてモデル内の推定パラメータ(重み)を計算する。その後、その推定パラメータを用いて、検証データに対する誤差を求める。

 線形回帰モデルのパラメータは最小二乗法で推定を行う。最小二乗法では平均二乗誤差を最小とするパラメータを探索する。平均二乗誤差(Mean Squared Error; MSE)はデータとモデル出力の二乗誤差の和である。

$$MSE_{{\rm train}}=\frac{1}{n_{{\rm train}}}\sum_{i=1}^{n_{{\rm train}}}(\hat{y}_i^{({\rm train})}-y_i^{({\rm train})})^2$$

学習データの平均二乗誤差の最小化は、その勾配が0になる点を求める。MSEを最小にするwは

$$\hat{{\bf w}} = \arg\min_{{\bf w}\in\mathbb{R}^{m+1}}{\rm MSE}_{{\rm train}} $$

と表現できる。これを計算すると、

$$ \hat{{\bf w}}=(X^{({\rm train})T}X^{({\rm train})})^{-1}X^{({\rm train})T}{\bf y}^{({\rm train})}$$

となる。以上により、学習データが平均二乗誤差が最小となるパラメータ(重み)が決定される。このパラメータを用いることで学習後の線形回帰モデルによる予測を行うことが可能となる。

$$\begin{eqnarray}回帰係数&:&\hat{{\bf w}}=(X^{({\rm train})T}X^{({\rm train})})^{-1}X^{({\rm train})T}{\bf y}^{({\rm train})} \\
予測値&:&\hat{{\bf y}}=X_{n_{{\rm new}}\times m+1}(X^{({\rm train})T}X^{({\rm train})})^{-1}X^{({\rm train})T}{\bf y}^{({\rm train})}\end{eqnarray}$$

実装演習

ボストンの住宅データモデルセットを線形回帰モデルで分析する。
部屋数が4で犯罪率が0.3の物件はいくらになるか?

ソースコード例

from sklearn.datasets import load_boston
from sklearn.linear_model import LinearRegression
from pandas import DataFrame
import numpy as np

#ボストンデータのインポート
boston = load_boston()

#DataFrameに変換
df = DataFrame(data=boston.data, columns = boston.feature_names)
#目的変数をDataFrameに追加
df['PRICE'] = np.array(boston.target)

#説明変数(部屋数と犯罪率)
data = df.loc[:, ['CRIM', 'RM']].values
#目的変数(物件の価格)
target = df.loc[:, 'PRICE'].values

model = LinearRegression()

#fit関数でパラメータ推定を行う
model.fit(data,target)

#predict関数で条件にあったパラメータ推定
print(model.predict([[0.3,4]]))

出力

[4.24007956]

動かしてみるとわかるが、この場合の学習では犯罪数と部屋の価格には相関はなく、部屋数に大きく依存している。そもそも指標として学習させるべきなのかも事前にある程度は考慮してもいいのかもしれない(現実には、人間ではバイアスがかかって判断できないことを機械に計算させているのであるから全てトライ&エラーにするべきだとは思う)。

2. 非線形回帰モデル

現実の現象は非線形である場合が多い。ここでは線形モデルによる非線形回帰について整理する。

基底展開法

基底展開法は、基底関数(φ(x))とハイパーパラメータベクトルの線型結合を使用する。ここで基底関数は既知の非線形関数であり、多項式関数やガウス型基底関数、スプライン関数などがある。未知パラメータは線形回帰モデルと同様に最小二乗法や最尤法により決定する。

$$ y_i = w_0 + \sum_{j=1}^m w_j\phi_j({\bf x}_i)+\epsilon_i $$

以下に示すとおり、基本的な考え方は線形回帰と同じである。

$$\begin{eqnarray}
説明変数&:&{\bf x}_i =(x_{i1}, x_{i2}, \cdot\cdot\cdot, x_{im}) \in \mathbb{R}^m \\
非線形関数ベクトル&:& {\bf\phi}({\bf x}_i) = (\phi_1({\bf x}_i), \phi_2({\bf x}_i), \cdot\cdot\cdot, \phi_k({\bf x}_i))^T \in \mathbb{R}^k\\
非線形関数の計画行列&:& \Phi^{({\rm train})} = ({\bf \phi}({\bf x}_1), {\bf \phi}({\bf x}_2), \cdot\cdot\cdot, {\bf \phi}({\bf x}_n))^T \in \mathbb{R}^{n\times k}\\
最尤法による予測値&:& \hat{{\bf y}}=\Phi(\Phi^{({\rm train})T}\Phi^{({\rm train})})^{-1}\Phi^{({\rm train})T}){\bf y}^{({\rm train})}
\end{eqnarray}$$

単回帰でのXがΦに変更されたイメージである。

未学習と過学習

  • 未学習:モデルの表現力が低い
  • 過学習:学習データに対する誤差は小さいが、テストデータに対する誤差が大きいモデル

 未学習と過学習はモデルの表現力に応じて生じる現象であり、未学習に関しては表現力の高いモデルを利用することが対策となる。現実的には過学習がよくある現象であり、さまざまな対策方がある。

 過学習の対策として正則化法(罰則化法)についてまとめる。正則化法とはモデルの複雑さに伴って、値が大きくなる正則化項(罰則化項)を課した関数を最小化する方法である。

$$ S_{\gamma} = ({\bf y} - \Phi^{n\times k}{\bf w})^T({\bf y} - \Phi {\bf w}) + \gamma R({\bf w}) $$

ここで、第1項は基底関数の数(k)が増えると残渣は減るが、モデルは複雑になることを表し、第2項はモデルの複雑さに伴う罰則を表している。具体的な手法としてはL2ノルム(Ridge推定)と、L1ノルム(Lasso推定)がある。

データ分割法

 線形回帰モデルでは、単純なデータ分割法のみ説明したが、モデルが複雑になるほど、より効率よく学習を進めて汎化させる必要がある。有限のデータを学習前に学習用とテスト用の2つに分割し、以降買えない方法をホールドアウト法という。これではデータ量が少ない場合は、学習/テストのどちらかが不十分になる。
 一方で、イテレーションを行う検証法を交差検証法(クロスバリデーション)という。例えば、
前データを5分割して学習用とテスト用を4:1の割合でモデルを作成する(計5つ)
→それぞれのモデルで精度を計算し、平均を「CV値」とする
→同じことを複数の「回帰モデル」で行い、最もCV値が小さいモデルを採用する
という流れで検証を行う。
 最後に全てのチューニングパラメータの組み合わせて評価値を算出する方式をグリッドサーチという。例えば(λ, m) がハイパーパラメータだとすると、各要素(λi, mi)に対して探査を行う方式である。

実装演習

  まずは、真の関数からノイズと伴うデータを生成する。

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

n=200

def true_func(x):
    z = 1-48*x+218*x**2-315*x**3+145*x**4
    return z

#真の関数からノイズを伴うデータを生成
data = np.random.rand(n).astype(np.float32)
data = np.sort(data)
target = true_func(data)

#ノイズを加える
noise = 0.5 * np.random.randn(n)
target = target + noise

#plot
plt.scatter(data,target, c='blue')
plt.title('NonLinear Regression')
plt.grid()
plt.savefig("data_with_noise.png")

 非線形回帰を行うべきなのは明らかであるが、線形回帰を行った結果を示す。

#線形回帰
from sklearn.linear_model import LinearRegression

clf = LinearRegression()
data = data.reshape(-1,1)
target = target.reshape(-1,1)
clf.fit(data, target)

p_lin = clf.predict(data)

plt.scatter(data, target, c='blue', label='data')
plt.plot(data, p_lin, c='darkorange', label = 'linear regression')
plt.legend()
plt.savefig("linear_regression.png")
print(clf.score(data, target))

次に基底展開法で多項式関数を用いた場合の結果を示す。

#基底展開法(多項式関数)
from sklearn.preprocessing import PolynomialFeatures
from sklearn.pipeline import Pipeline

deg = [1,3,5,7,9]
plt.scatter(data, target, label='data',c='blue')
for d in deg:
    regr = Pipeline([
        ('poly', PolynomialFeatures(degree=d)),
        ('linear', LinearRegression())
    ])
    regr.fit(data,target)
    p_poly = regr.predict(data)
    plt.plot(data, p_poly, label='polynomial of degree %d' %(d))
plt.grid()
plt.legend()
plt.savefig("polynomial.png")

最後にLasso推定を行った結果を示す。

Lasso推定
from sklearn.metrics.pairwise import rbf_kernel
from sklearn.linear_model import Lasso

kx = rbf_kernel(X=data, Y=data, gamma=5)

lasso_clf = Lasso(alpha=0.01,max_iter=100000)
lasso_clf.fit(kx, target)

p_lasso = lasso_clf.predict(kx)

plt.scatter(data, target, c='blue')
plt.plot(data, p_lasso)
plt.grid()
plt.savefig("Lasso.png")

print(lasso_clf.score(kx, target))

3. ロジスティック回帰モデル

ロジスティック回帰モデル

 ロジスティック回帰モデルは、分類問題を解くための教師あり機械学習モデルである。分類問題とは入力からクラスに分類する問題である。このため入力(説明変数、特微量)はm次元のベクトル。出力(目的変数)は0または1となる。

$$\begin{eqnarray}
説明変数&:&{\bf x} = (x_1, x_2, \cdot\cdot\cdot, x_m)^T \in \mathbb{R}^m\\
目的変数&:& y\in \{0,1\}
\end{eqnarray}$$

出力を0または1とするために、入力とm次元パラメータの線型結合をシグモイド関数に入力する。シグモイド関数の入力は実数、出力は0〜1となる単調増加関数である。

$$\begin{eqnarray}
\sigma(x) &=& \frac{1}{1+\exp(-ax)}\\
\frac{\partial \sigma(x)}{\partial x}&=&a\sigma(x)(1-\sigma(x))
\end{eqnarray}$$

シグモイド関数の出力が確率に対応している。

$$ P(Y=1|{\bf x}) = \sigma(w_0+w_1 x_1+\cdot\cdot\cdot+w_m x_m) $$

データYは確率が0.5以上ならば1・未満なら0と予測する。

最尤推定

 最尤推定とは、データからそのデータを推定したであろう尤もらしい分布(パラメータ)を推定することを指す。つまり、尤度関数(L(w))を最大化するようなパラメータを選ぶ推定方法である。尤度関数はデータを固定し、パラメータのみを変数とし、尤もらしさを表す数値を出力する関数である。

$$\begin{eqnarray}
P(y_1, y_2, \cdot\cdot\cdot, y_n | w_0, w_1, \cdot\cdot\cdot, w_m) &=& \prod^n_{i=1}p_i^{y_i}(1-p_i)^{1-y_i}\\
&=&\prod_{i=1}^n\sigma({\bf w}^T{\bf x}_i)^{y_i}(1-\sigma({\bf w}^T{\bf x}_i))^{1-y_i}\\
&=&L({\bf w})
\end{eqnarray}$$

尤度関数P(y|w)は確率変数が独立を仮定することで確率の積に分解可能となる。これにより、尤度関数はパラメータのみに依存する関数L(w)となる。尤度関数を最大化する上で対数を取ることで微分の計算を簡単にしたものを対数尤度関数という。対数関数は単調増加であるため、対数尤度関数が最大になる点と尤度関数が最大になる点は同じである。

$$\begin{eqnarray}
E(w_0, w_1, \cdot\cdot\cdot, w_m) &=& -\log L(w_0, w_1, \cdot\cdot\cdot, w_m)\\
&=& -\sum_{i=1}^n\{y_i \log p_i +(1-y_i)\log (1-p_i)\}
\end{eqnarray}$$

尤度関数の微分=0の解析解はないため、イテレーションで近傍解を求める必要がある。ここで用いられるのが勾配降下法(Gradient descent)である。

$$\begin{eqnarray}
{\bf w}(k+1) &=& {\bf w}^k -\eta\frac{\partial E({\bf w})}{\partial {\bf w}}\\
\frac{\partial E({\bf w})}{\partial {\bf w}}&=&-\sum^n_{i=1}(y_i-p_i){\bf x}_i
\end{eqnarray}$$

ここでηは学習率を表す。パラメータが更新されなくなった場合、そこが(探査範囲では)最適解ということになる。
 ただし、勾配降下法では、パラメータを更新するのにN個全てのデータに対する和を求める必要がある。そこで計算時間の削減のために、データを確率的に選んでパラメータを更新する方法を確率的勾配降下法という。

$${\bf w}(k+1) = {\bf w}^k + \eta(y_i-p_i){\bf x}_i$$

これにより、勾配降下法でパラメータを1回更新することと同じ計算量でパラメータをn回更新することが可能となり、効率よく最適な解を探査可能となる。

実装演習

タイタニックの乗客データを利用し、ロジスティック回帰モデルを作成する。
年齢が30歳の男性の乗客は生き残ることが可能であるか?
(源泉データ:https://www.kaggle.com/c/titanic/data?select=train.csv

import pandas as pd
from pandas import DataFrame

#csvの読み込み
titanic_df = pd.read_csv('../../Stage2_ML/study_ai_ml_google/data/titanic_train.csv')

###データ整形
#予測に不要なカラムのドロップ
titanic_df.drop(['PassengerId', 'Name', 'Ticket', 'Cabin'], axis=1, inplace=True)
#nullを中央値で補完
titanic_df['AgeFill'] = titanic_df['Age'].fillna(titanic_df['Age'].mean())

###実装
from sklearn.linear_model import LogisticRegression
#Genderを数値に変換(female:0, male:1)
titanic_df['Gender'] = titanic_df['Sex'].map({'female': 0, 'male': 1}).astype(int)
#ロジスティック回帰
data = titanic_df.loc[:, ["AgeFill", "Gender"]].values
label = titanic_df.loc[:,["Survived"]].values
model = LogisticRegression()

model.fit(data, label)
print(model.predict([[30,1]]))
###出力
[0]

出力としては0(生存できない)ことになる。ただし、今回は年齢と性別のみしか学習させていないため、情報としては不十分である。例えばチケットのクラスや、家族構成によっても異なるかもしれない。本来は相関性を見ながら、必要十分と思われるデータを入力に学習させるべきである。

4. 主成分分析

主成分分析

 主成分分析とは、多変量データの持つ構造をより少数個の指標に圧縮する分析手法である。具体的には線形変換後の変数の分散が最大となる射影軸を探索する手法である。分散が最大となるということはももっとも多くデータの情報が残ることを意味しており、情報を損失せずに次元数を削減することを意味している。

$$\begin{eqnarray}
データ行列&:&\bar{\bf {X}}=({\bf x}_1-\bar{{\bf x}},\cdot\cdot\cdot,{\bf x}_n-\bar{{\bf x}})^T\in \mathbb{R}^{n\times m}\\
線形変換後のベクトル&:&s_j =(s_{1j},\cdot\cdot\cdot,s_{nj})^T=\bar{X}a_j;\ a_j\in\mathbb{R}^m
\end{eqnarray}$$

ただし、jは射影軸のインデックスを表す。
 第k主成分の分散の全分散に対する割合を寄与率という。平たく言うと、第k主成分がもつ情報量の割合である。また、1〜k主成分まで圧縮した際の情報損失量の割合を累積寄与率という。

$$\begin{eqnarray}
寄与率&:& c_k = \frac{\lambda_k}{\sum_{i=1}^m \lambda_i}\\
累積寄与率&:& = r_k = \frac{\sum_{j=1}^k \lambda_j}{\sum_{i=1}^m\lambda_i}
\end{eqnarray}$$

実装演習

乳がん検査データを利用してロジスティック回帰モデルを作成。
32次元のデータを2次元上に次元圧縮した際に、うまく判別できるかを確認する。
(源泉データ:https://www.kaggle.com/uciml/breast-cancer-wisconsin-data

import pandas as pd
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegressionCV
from sklearn.metrics import confusion_matrix
from sklearn.decomposition import PCA
import matplotlib.pyplot as plt


cancer_df = pd.read_csv('../../Stage2_ML/study_ai_ml_google/data/cancer.csv')

#不要なデータの削除
cancer_df.drop('Unnamed: 32', axis=1, inplace=True)

#目的変数の抽出
y = cancer_df.diagnosis.apply(lambda d:1 if d == 'M' else 0)

#説明変数の抽出
X = cancer_df.loc[:, 'radius_mean':]

#学習用とテスト用でデータを分割
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=0)
#標準化
scalar = StandardScaler()
X_train_scaled = scalar.fit_transform(X_train)
X_test_scaled = scalar.transform(X_test)
#ロジスティック回帰で学習
logistic = LogisticRegressionCV(cv=10, random_state=0, max_iter=5000)
logistic.fit(X_train_scaled, y_train)
#検証
print('Train score:{:.3f}'.format(logistic.score(X_train_scaled, y_train)))
print('Test score: {:.3f}'.format(logistic.score(X_test_scaled, y_test)))
print('Confusion matrix:\n{}'.format(confusion_matrix(y_true=y_test, y_pred=logistic.predict(X_test_scaled))))

###PCA
pca = PCA(n_components=30)
pca.fit(X_train_scaled)
plt.bar([n for n in range(1, len(pca.explained_variance_ratio_)+1)], pca.explained_variance_ratio_)
plt.savefig("pca01.png")
plt.clf()

#次元数2まで圧縮
pca = PCA(n_components=2)
X_train_pca = pca.fit_transform(X_train_scaled)
print('X_train_pca shape: {}'.format(X_train_pca.shape))
#寄与率
print('explained variance ratio: {}'.format(pca.explained_variance_ratio_))
#散布図にプロット
temp = pd.DataFrame(X_train_pca)
temp['Outcome'] = y_train.values
b = temp[temp['Outcome'] == 0]
m = temp[temp['Outcome'] == 1]
plt.scatter(x=b[0], y=b[1], marker='o')
plt.scatter(x=m[0], y=m[1], marker='^')
plt.xlabel('PC 1')
plt.ylabel('PC 2')
plt.savefig("pca02.png")

#検証結果(出力)
Train score:0.988
Test score: 0.972
Confusion matrix:
[[89  1]
 [ 3 50]]

第1・2主成分がデータのほとんどの情報を持っていることがわかる。

#2次元に圧縮後の出力
X_train_pca shape: (426, 2)
explained variance ratio: [0.43315126 0.19586506]

上記より、第1主成分の寄与率は0.433、第2主成分の寄与率は0.196となる。可視化しても2次元で判別がある程度は可能であることがわかる。

5. アルゴリズム

k近傍法(kNN)

 k近傍法(kNN)は分類問題のための機械学習手法である。分類方法は、最近傍のデータをk個取ってきて、それらがもっとも多く所属するクラスに識別する。イメージとしては、近傍のデータを用いて多数決を行う手法である。
 kを大きくするとクラスの境界(決定境界)は滑らかになる。

実装演習(k近傍法)

人口データを分類する(仮想的にここでは乱数でクラス分類を行う)。

 訓練データの作成

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

###訓練データ生成
def gen_data():
    x0 = np.random.normal(size=50).reshape(-1,2) - 1.
    x1 = np.random.normal(size=50).reshape(-1,2) + 1.
    x_train = np.concatenate([x0, x1])
    y_train = np.concatenate([np.zeros(25), np.ones(25)]).astype(np.int)
    return x_train, y_train

X_train, ys_train = gen_data()
plt.scatter(X_train[:,0], X_train[:,1], c=ys_train)
plt.savefig("kNN_train.png")
plt.clf()

予測

###予測
def distance(x1, x2):
    return np.sum((x1 - x2)**2, axis=1)

def knc_predict(n_neighbors, x_train, y_train, X_test):
    y_pred = np.empty(len(X_test), dtype=y_train.dtype)
    for i, x in enumerate(X_test):
        distances = distance(x, X_train)
        nearest_index = distances.argsort()[:n_neighbors]
        mode, _ = stats.mode(y_train[nearest_index])
        y_pred[i] = mode
    return y_pred

def plt_result(x_train, y_train, y_pred,k):
    xx0, xx1 = np.meshgrid(np.linspace(-5,5,100),np.linspace(-5,5,100))
    xx = np.array([xx0, xx1]).reshape(2,-1).T
    plt.scatter(x_train[:,0], x_train[:,1], c=y_train)
    plt.contour(xx0, xx1, y_pred.reshape(100, 100).astype(dtype=np.float),alpha=0.2, levels=np.linspace(0,1,3))
    plt.savefig('kNN_results_k%d.png' %(k))
    plt.clf()

n_neighbors = 3
xx0, xx1 = np.meshgrid(np.linspace(-5,5,100), np.linspace(-5,5,100))
X_test = np.array([xx0, xx1]).reshape(2,-1).T

y_pred = knc_predict(n_neighbors, X_train, ys_train, X_test)
plt_result(X_train, ys_train, y_pred, n_neighbors)

 【実行結果】



k=3の場合の予測結果

k=10の場合の予測結果

k-means

 k-平均法(k-means)は与えられたデータをk個のクラスタに分けるクラスタリング手法である。特徴の似ているもの同士をグループ化するので、教師なし学習である。手順は以下のとおりである。

  1. 各クラスタ中心の初期値を設定する
  2. 各データ点に対して、各クラスタ中心との距離を計算し、もっとも距離が近いクラスタを割り当てる。
  3. 各クラスタの平均ベクトル(中心)を計算する。
  4. 収束するまで2,3の処理を繰り返す。

 ただし、中心の初期値によってクラスタリング結果は変わりうる。初期値が離れればうまくクラスタリングできるが、初期値が近いうまくクラスタリングできない。また分類するクラスタの数であるkの値によってもクラスタリング結果は変わるため注意が必要である。

実装演習(k-means)

データ生成

import numpy as np
import matplotlib.pyplot as plt

###データ生成
def gen_data():
    x1 = np.random.normal(size=(100,2)) + np.array([-5, 5])
    x2 = np.random.normal(size=(100,2)) + np.array([ 5,-5])
    x3 = np.random.normal(size=(100,2)) + np.array([ 0, 5])
    return np.vstack((x1,x2,x3))

X_train = gen_data()
plt.scatter(X_train[:,0], X_train[:,1])
plt.savefig("kmeans_train.png")
plt.clf()

学習・クラスタリング結果

###学習
def distance(x1, x2):
    return np.sum((x1-x2)**2, axis=1)

n_clusters = 3
iter_max = 1000

#中心の初期化
centers = X_train[ np.random.choice(len(X_train), n_clusters,replace = False) ]

for _ in range(iter_max):
    prev_centers = np.copy(centers)
    D = np.zeros((len(X_train), n_clusters))
    for i, x in enumerate(X_train):
        D[i] = distance(x,centers)
    cluster_index = np.argmin(D, axis = 1)
    for k in range(n_clusters):
        index_k = cluster_index == k
        centers[k] = np.mean(X_train[index_k],  axis=0)
    #収束判定
    if np.allclose(prev_centers, centers):
        break

###クラスタリング結果
def plt_result(X_train, centers, xx):
    plt.scatter(X_train[:,0], X_train[:,1],c=y_pred, cmap='spring')
    plt.scatter(centers[:,0], centers[:,1],s=200, marker='X', lw=2,c='black',edgecolor="white")
    pred = np.empty(len(xx), dtype=int)
    for i, x in enumerate(xx):
        d = distance(x, centers)
        pred[i] = np.argmin(d)
    plt.contourf(xx0, xx1, pred.reshape(100,100), alpha=0.2, cmap='spring')
    plt.savefig("kmeans_result.png")
    plt.clf()

y_pred = np.empty(len(X_train), dtype=int)
for i, x in enumerate(X_train):
    d = distance(x, centers)
    y_pred[i] = np.argmin(d)

xx0, xx1 = np.meshgrid(np.linspace(-10, 10, 100), np.linspace(-10, 10, 100))
xx = np.array([xx0, xx1]).reshape(2, -1).T

plt_result(X_train, centers, xx)

6. サポートベクターマシーン(SVM)

SVM

 SVMでは、クラスのデータが分類境界からなるべく離れるようにして「最適」な分類境界を決定する。分類境界を挟んでそれぞれのクラスがどのくらい離れているかをマージンと呼ぶ。SVMではマージンが大きいほど、良い分類境界であると考える。このような手法をマージン最大化という。訓練データを完璧に分類できる決定関数が存在する仮定を置いた場合をハードマージン、この仮定をなくした場合をソフトマージンという。 ソフトマージンでは誤差を許容するために、スラック変数を導入する。
 データの分布によっては、線形分布ではうまく分類できないケースが生じることがある。このようなデータに対して、特微ベクトルを非線形変換して、その空間での線形分類を行うカーネルトリックという。カーネルトリックでは特微空間に写像した関数(φ(x))からカーネル関数を構成する。

$$ K({\bf x}_i, {\bf x}_j) = \phi({\bf x}_i)^T\phi({\bf x}_j)$$

カーネル関数を用いることで、最適化や決定関数の計算にはφ(x) を用いる必要がなくなるため、計算が簡単になる。カーネル関数には多項式カーネル、ガウスカーネル、シグモイドカーネルなどがあり、これらを用いることで、入力データの非線形分離が可能になる。

実装演習

①線形分離可能な場合

import numpy as np
import matplotlib.pyplot as plt

###線形分離可能な場合
def gen_data():
    x0 = np.random.normal(size=50).reshape(-1, 2) - 2.
    x1 = np.random.normal(size=50).reshape(-1, 2) + 2.
    X_train = np.concatenate([x0, x1])
    ys_train = np.concatenate([np.zeros(25), np.ones(25)]).astype(np.int)
    return X_train, ys_train

X_train, ys_train = gen_data()
plt.scatter(X_train[:, 0], X_train[:, 1], c=ys_train)
plt.savefig("svm_linear_train.png")
plt.clf()

#学習
t = np.where(ys_train == 1.0, 1.0, -1.0)
n_samples = len(X_train)
K = X_train.dot(X_train.T)

eta1 = 0.01
eta2 = 0.001
n_iter = 500

H = np.outer(t, t) * K
a = np.ones(n_samples)
for _ in range(n_iter):
    grad = 1 - H.dot(a)
    a += eta1 * grad
    a -= eta2 * a.dot(t) * t
    a = np.where(a > 0, a, 0)

#予測
index = a > 1e-6
support_vectors = X_train[index]
support_vector_t = t[index]
support_vector_a = a[index]

term2 = K[index][:, index].dot(support_vector_a * support_vector_t)
b = (support_vector_t - term2).mean()

xx0, xx1 = np.meshgrid(np.linspace(-5, 5, 100), np.linspace(-5, 5, 100))
xx = np.array([xx0, xx1]).reshape(2, -1).T

X_test = xx
y_project = np.ones(len(X_test)) * b
for i in range(len(X_test)):
    for a, sv_t, sv in zip(support_vector_a, support_vector_t, support_vectors):
        y_project[i] += a * sv_t * sv.dot(X_test[i])
y_pred = np.sign(y_project)

#訓練データ
plt.scatter(X_train[:, 0], X_train[:, 1], c=ys_train)
#サポートベクトル
plt.scatter(support_vectors[:, 0], support_vectors[:, 1],
                    s=100, facecolors='none', edgecolors='k')
#領域
plt.contour(xx0, xx1, y_project.reshape(100, 100), colors='k',
                     levels=[-1, 0, 1], alpha=0.5, linestyles=['--', '-', '--'])
#マージンと決定領域
plt.quiver(0, 0, 0.1, 0.35, width=0.01, scale=1, color='pink')
plt.savefig("svm_linear_result.png")
plt.clf()

【実行結果】

訓練データ

予測結果

②線形分離不可

###線形分離不可
factor = .2
n_samples = 50
linspace = np.linspace(0, 2 * np.pi, n_samples // 2 + 1)[:-1]
outer_circ_x = np.cos(linspace)
outer_circ_y = np.sin(linspace)
inner_circ_x = outer_circ_x * factor
inner_circ_y = outer_circ_y * factor

X = np.vstack((np.append(outer_circ_x, inner_circ_x),
               np.append(outer_circ_y, inner_circ_y))).T
y = np.hstack([np.zeros(n_samples // 2, dtype=np.intp),
               np.ones(n_samples // 2, dtype=np.intp)])
X += np.random.normal(scale=0.15, size=X.shape)
x_train = X
y_train = y

plt.scatter(x_train[:,0], x_train[:,1], c=y_train)
plt.savefig("svm_nonlinear_train.png")
plt.clf()

#学習(ガウシアンカーネル)
def rbf(u, v):
        sigma = 0.8
        return np.exp(-0.5 * ((u - v)**2).sum() / sigma**2)
    
X_train = x_train
t = np.where(y_train == 1.0, 1.0, -1.0)

n_samples = len(X_train)
# RBFカーネル
K = np.zeros((n_samples, n_samples))
for i in range(n_samples):
    for j in range(n_samples):
        K[i, j] = rbf(X_train[i], X_train[j])

eta1 = 0.01
eta2 = 0.001
n_iter = 5000

H = np.outer(t, t) * K

a = np.ones(n_samples)
for _ in range(n_iter):
    grad = 1 - H.dot(a)
    a += eta1 * grad
    a -= eta2 * a.dot(t) * t
    a = np.where(a > 0, a, 0)

#予測
index = a > 1e-6
support_vectors = X_train[index]
support_vector_t = t[index]
support_vector_a = a[index]

term2 = K[index][:, index].dot(support_vector_a * support_vector_t)
b = (support_vector_t - term2).mean()

xx0, xx1 = np.meshgrid(np.linspace(-1.5, 1.5, 100), np.linspace(-1.5, 1.5, 100))
xx = np.array([xx0, xx1]).reshape(2, -1).T

X_test = xx
y_project = np.ones(len(X_test)) * b
for i in range(len(X_test)):
    for a, sv_t, sv in zip(support_vector_a, support_vector_t, support_vectors):
        y_project[i] += a * sv_t * rbf(X_test[i], sv)
y_pred = np.sign(y_project)
# 訓練データ
plt.scatter(x_train[:, 0], x_train[:, 1], c=y_train)
# サポートベクトル
plt.scatter(support_vectors[:, 0], support_vectors[:, 1],
                    s=100, facecolors='none', edgecolors='k')
# 領域
plt.contourf(xx0, xx1, y_pred.reshape(100, 100), alpha=0.2, levels=np.linspace(0, 1, 3))
# マージンと決定境界
plt.contour(xx0, xx1, y_project.reshape(100, 100), colors='k',
                     levels=[-1, 0, 1], alpha=0.5, linestyles=['--', '-', '--'])
plt.savefig("svm_nonlinear_result.png")
plt.clf()

【実行結果】

訓練データ

予測結果

③ソフトマージン

###ソフトマージンSVM
x0 = np.random.normal(size=50).reshape(-1, 2) - 1.
x1 = np.random.normal(size=50).reshape(-1, 2) + 1.
x_train = np.concatenate([x0, x1])
y_train = np.concatenate([np.zeros(25), np.ones(25)]).astype(np.int)
plt.scatter(x_train[:, 0], x_train[:, 1], c=y_train)
plt.savefig("svm_softmargin_train.png")
plt.clf()

#学習
X_train = x_train
t = np.where(y_train == 1.0, 1.0, -1.0)

n_samples = len(X_train)
# 線形カーネル
K = X_train.dot(X_train.T)

C = 1
eta1 = 0.01
eta2 = 0.001
n_iter = 1000

H = np.outer(t, t) * K

a = np.ones(n_samples)
for _ in range(n_iter):
    grad = 1 - H.dot(a)
    a += eta1 * grad
    a -= eta2 * a.dot(t) * t
    a = np.clip(a, 0, C)

#予測
index = a > 1e-8
support_vectors = X_train[index]
support_vector_t = t[index]
support_vector_a = a[index]

term2 = K[index][:, index].dot(support_vector_a * support_vector_t)
b = (support_vector_t - term2).mean()

xx0, xx1 = np.meshgrid(np.linspace(-4, 4, 100), np.linspace(-4, 4, 100))
xx = np.array([xx0, xx1]).reshape(2, -1).T

X_test = xx
y_project = np.ones(len(X_test)) * b
for i in range(len(X_test)):
    for a, sv_t, sv in zip(support_vector_a, support_vector_t, support_vectors):
        y_project[i] += a * sv_t * sv.dot(X_test[i])
y_pred = np.sign(y_project)

# 訓練データ
plt.scatter(x_train[:, 0], x_train[:, 1], c=y_train)
# サポートベクトル
plt.scatter(support_vectors[:, 0], support_vectors[:, 1],
                    s=100, facecolors='none', edgecolors='k')
# 領域
plt.contourf(xx0, xx1, y_pred.reshape(100, 100), alpha=0.2, levels=np.linspace(0, 1, 3))
# マージンと決定境界
plt.contour(xx0, xx1, y_project.reshape(100, 100), colors='k',
                     levels=[-1, 0, 1], alpha=0.5, linestyles=['--', '-', '--'])
plt.savefig("svm_softmargin_result.png")

【実行結果】

訓練データ

予測結果

まとめ

 今回は機械学習の全般について整理した。各モデルの内部構造を理解しておくと、実装演習はそつなくこなせるはずである。現実問題としては、機械学習で何かがしたい際に、どういった手段があるか、メリット・デメリットは何であるかを知っておくことが重要となる。実装演習を含めてここでまとめたことを徐々に深掘りしつつ、定着を図っていく。

コメント

タイトルとURLをコピーしました