機械学習【ラビットチャレンジ】

概要

本記事は「JDLA E資格」の取得に必要なJDLA認定プログラムのひとつであるラビットチャレンジの受講レポートとして投稿したものです。

E資格について
https://www.jdla.org/certificate/engineer/

目次

機械学習
線形回帰
線形回帰モデル
線形回帰モデル・ハンズオン
非線形回帰モデル
非線形回帰モデル・ハンズオン
ロジスティック回帰モデル
ロジスティック回帰モデル・ハンズオン
主成分分析
主成分分析・ハンズオン
k近傍法(kNN)
k平均法(k-means)
k近傍法・ハンズオン
SVM（サポートベクターマシン）
SVM・ハンズオン

機械学習

・人がプログラムするのは認識の仕方ではなく学習の仕方（数学で記述）

機械学習モデリングプロセス

1. 問題設定
2. データ選定
3. データの前処理
4. 機械学習モデルの選定
5. モデルの学習(パラメータ推定)
6. モデルの評価

・実務上重要なのは1.の問題設定であることも多いが、ここでは4.がメインテーマとなる。
・GIGO(Garbage in, garbage out)： 「ゴミを入れるとゴミが出てくる」

機械学習の分類

機械学習は大きく分けて教師あり学習と教師なし学習があるが、ここでは主に教師あり学習を扱う。

線形回帰

線形： ざっくり言うと比例

回帰問題

・ある入力(離散あるいは連続値)から出力(連続値)を予測する問題
直線で予測→線形回帰
線で予測→非線形回帰

回帰で扱うデータ

・入力（説明変数・特徴量）： m次元のベクトル(m=1の場合はスカラー)
・出力（目的変数）： スカラー

線形回帰モデル

・回帰問題を解くための機械学習モデルの1つ
・教師あり学習(教師データから学習)
・入力とm次元パラメータの線形結合を出力するモデル
・慣例として予測値にはハットを付ける(正解データとは異なる)

・説明変数が1次元の場合(m=1)： 単回帰モデル
・データへの仮定： データは回帰直線に誤差が加わり観測されていると仮定

・説明変数が多次元の場合(m>1)： 線形重回帰モデル

・学習用と検証用にデータを分割する

\begin{align}
\\
\hat{y} &= \boldsymbol{w}^T \boldsymbol{x} + w_0 \\
&= \sum_{j=1}^{m} w_j x_j + w_0
\end{align}

最小二乗法

線形回帰モデルのパラメータは最小二乗法で推定する

・線形回帰モデルのパラメータは最小二乗法で推定
・平均二乗誤差(残差平方和)： データとモデル出力の二乗誤差の和
・最小二乗法： 学習データの平均二乗誤差を最小とするパラメータを探索
・勾配が0になる点を求める
・二乗する理由： 最小値が求めやすい（微分して0になる点）（ただし外れ値に弱い）

線形回帰モデル・ハンズオン

・ボストンハウジングデータを利用

単回帰分析

コード

# from モジュール名 import クラス名（もしくは関数名や変数名）
from sklearn.datasets import load_boston
from pandas import DataFrame
import numpy as np
## sklearnモジュールからLinearRegressionをインポート
from sklearn.linear_model import LinearRegression

# ボストンデータを"boston"というインスタンスにインポート
boston = load_boston()

# 説明変数らをDataFrameへ変換
df = DataFrame(data=boston.data, columns = boston.feature_names)

# 目的変数をDataFrameへ追加
df['PRICE'] = np.array(boston.target)

# 説明変数
data = df.loc[:, ['RM']].values

# 目的変数
target = df.loc[:, 'PRICE'].values

# オブジェクト生成
model = LinearRegression()

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

# 予測
model.predict([[1]])

結果

array([-25.5685118])

・部屋数1のときマイナスの家賃予測値が返ってきている。
・元のデータに部屋数1の場合が含まれないことが考えられる→外挿問題

重回帰分析

コード

# 説明変数
data2 = df.loc[:, ['CRIM', 'RM']].values
# 目的変数
target2 = df.loc[:, 'PRICE'].values
# オブジェクト生成
model2 = LinearRegression()
# fit関数でパラメータ推定
model2.fit(data2, target2)
# 予測
model2.predict([[0.2, 7]])

結果

array([29.43977562])

回帰係数と切片の値を確認

# 単回帰の回帰係数と切片を出力
print('推定された回帰係数: %.3f, 推定された切片 : %.3f' % (model.coef_, model.intercept_))
# 重回帰の回帰係数と切片を出力
print(model.coef_)
print(model.intercept_)

結果

推定された回帰係数: 9.102, 推定された切片 : -34.671
[9.10210898]
-34.67062077643857

非線形回帰モデル

基底展開法

・回帰関数として、基底関数と呼ばれる既知の非線形関数とパラメータベクトルの線型結合を使用
・未知パラメータは線形回帰モデルと同様に最小2乗法や最尤法により推定
・パラメータについては線形のまま（xをΦで飛ばすだけ）

y_i = w_0 + \sum_{j=1}^{m} w_j \phi_j(x_i) + \varepsilon_i

未学習(underfitting)と過学習(overfitting)

・未学習： 学習データに対して、十分小さな誤差が得られないモデル
→モデルの表現力が低いため、表現力の高いモデルを利用する

・過学習： 小さな誤差は得られたけど、テスト集合誤差との差が大きいモデル
→学習データの数を増やす
→不要な基底関数(変数)を削除して表現力を抑止
→正則化法を利用して表現力を抑止

・不要な基底関数を削除
→組み合わせ爆発が起こり難しいため、正則化法を用いる

正則化法（罰則化法）

「モデルの複雑さに伴って、その値が大きくなる正則化項(罰則項)を課した関数」を最小化

■Ridge推定量
・L2ノルムを利用
・縮小推定（パラメータを0に近づけるよう推定）

■Lasso推定量
・L1ノルムを利用
・スパース推定（いくつかのパラメータを正確に0に推定）

ホールドアウト法

・データを学習用とテスト用の2つに分割
・データが少ないと問題が起こりやすい（検証用に外れ値が入ってしまうなど）

クロスバリデーション(交差検証)

・データ分割を繰り返し、CV値（精度の平均）を取ってより精度の高いモデルを選択する

非線形回帰モデル・ハンズオン

コード

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

# seaborn設定
sns.set()
# 背景変更
sns.set_style("darkgrid", {'grid.linestyle': '--'})
# 大きさ(スケール変更)
sns.set_context("paper")

n=100

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

def linear_func(x):
    z = x
    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

# ノイズ付きデータを描画

plt.scatter(data, target)

plt.title('NonLinear Regression')
plt.legend(loc=2)

コード

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, label='data')
plt.plot(data, p_lin, color='darkorange', marker='', linestyle='-', linewidth=1, markersize=6, label='linear regression')
plt.legend()
print(clf.score(data, target))

コード

from sklearn.kernel_ridge import KernelRidge

clf = KernelRidge(alpha=0.0002, kernel='rbf')
clf.fit(data, target)

p_kridge = clf.predict(data)

plt.scatter(data, target, color='blue', label='data')

plt.plot(data, p_kridge, color='orange', linestyle='-', linewidth=3, markersize=6, label='kernel ridge')
plt.legend()

コード

# Ridge

from sklearn.metrics.pairwise import rbf_kernel
from sklearn.linear_model import Ridge

kx = rbf_kernel(X=data, Y=data, gamma=50)
# KX = rbf_kernel(X, x)

# clf = LinearRegression()
clf = Ridge(alpha=30)
clf.fit(kx, target)

p_ridge = clf.predict(kx)

plt.scatter(data, target,label='data')
for i in range(len(kx)):
    plt.plot(data, kx[i], color='black', linestyle='-', linewidth=1, markersize=3, label='rbf', alpha=0.2)

# plt.plot(data, p, color='green', marker='o', linestyle='-', linewidth=0.1, markersize=3)
plt.plot(data, p_ridge, color='green', linestyle='-', linewidth=1, markersize=3,label='ridge regression')
# plt.legend()

print(clf.score(kx, target))

コード

from sklearn.preprocessing import PolynomialFeatures
from sklearn.pipeline import Pipeline

# PolynomialFeatures(degree=1)

deg = [1,2,3,4,5,6,7,8,9,10]
for d in deg:
    regr = Pipeline([
        ('poly', PolynomialFeatures(degree=d)),
        ('linear', LinearRegression())
    ])
    regr.fit(data, target)
    # make predictions
    p_poly = regr.predict(data)
    # plot regression result
    plt.scatter(data, target, label='data')
    plt.plot(data, p_poly, label='polynomial of degree %d' % (d))

コード

# Lasso

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

kx = rbf_kernel(X=data, Y=data, gamma=5)
# KX = rbf_kernel(X, x)

# lasso_clf = LinearRegression()
lasso_clf = Lasso(alpha=10000, max_iter=1000)
lasso_clf.fit(kx, target)

p_lasso = lasso_clf.predict(kx)

plt.scatter(data, target)

# plt.plot(data, p, color='green', marker='o', linestyle='-', linewidth=0.1, markersize=3)
plt.plot(data, p_lasso, color='green', linestyle='-', linewidth=3, markersize=3)

print(lasso_clf.score(kx, target))

コード

from sklearn import model_selection, preprocessing, linear_model, svm

# SVR-rbf
clf_svr = svm.SVR(kernel='rbf', C=1e3, gamma=0.1, epsilon=0.1)
clf_svr.fit(data, target)
y_rbf = clf_svr.fit(data, target).predict(data)
 
# plot

plt.scatter(data, target, color='darkorange', label='data')
plt.plot(data, y_rbf, color='red', label='Support Vector Regression (RBF)')
plt.legend()
plt.show()

コード

from sklearn.model_selection import train_test_split
x_train, x_test, y_train, y_test = train_test_split(data, target, test_size=0.1, random_state=0)

from keras.callbacks import EarlyStopping, TensorBoard, ModelCheckpoint

cb_cp = ModelCheckpoint('/ディレクトリ/checkpoints/weights.{epoch:02d}-{val_loss:.2f}.hdf5', verbose=1, save_weights_only=True)
cb_tf  = TensorBoard(log_dir='/ディレクトリ/tensorBoard', histogram_freq=0)

def relu_reg_model():
    model = Sequential()
    model.add(Dense(10, input_dim=1, activation='relu'))
    model.add(Dense(1000, activation='relu'))
    model.add(Dense(1000, activation='relu'))
    model.add(Dense(1000, activation='relu'))
    model.add(Dense(1000, activation='relu'))
    model.add(Dense(1000, activation='relu'))
    model.add(Dense(1000, activation='relu'))
    model.add(Dense(1000, activation='relu'))
    model.add(Dense(1000, activation='linear'))
    model.add(Dense(1))

    model.compile(loss='mean_squared_error', optimizer='adam')
    return model

from keras.models import Sequential
from keras.layers import Input, Dense, Dropout, BatchNormalization
from keras.wrappers.scikit_learn import KerasRegressor

# use data split and fit to run the model
estimator = KerasRegressor(build_fn=relu_reg_model, epochs=100, batch_size=5, verbose=1)

history = estimator.fit(x_train, y_train, callbacks=[cb_cp, cb_tf], validation_data=(x_test, y_test))

y_pred = estimator.predict(x_train)

plt.title('NonLiner Regressions via DL by ReLU')
plt.plot(data, target, 'o')
plt.plot(data, true_func(data), '.')
plt.plot(x_train, y_pred, "o", label='predicted: deep learning')
# plt.legend(loc=2)

print(lasso_clf.coef_)

結果

[-0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0.
 -0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0.
  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]

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

■分類問題(クラス分類)
・ある入力(数値)からクラスに分類する問題

■分類で扱うデータ
・入力： m次元のベクトル(m=1の場合はスカラ)
・出力： 0 or 1の値

シグモイド関数

・シグモイド関数に入力することで0～1の値にしてくれる
・微分するとシグモイド関数自身で表現できるので便利

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

■シグモイド関数の出力をY=1になる確率に対応させる
・シグモイド関数に入力→出力が確率に対応

最尤推定

・ロジスティック回帰モデルではベルヌーイ分布を利用する

■尤度関数
・データは固定し、パラメータを変化させる
・尤度関数を最大化するようなパラメータを選ぶ推定方法を最尤推定という
・対数をとると微分の計算が簡単

■勾配降下法(Gradient descent)
・シグモイド関数が原因で解析解を求めることが困難→少しずつ最適解に近づけていく

ロジスティック回帰モデル・ハンズオン

・タイタニックの乗客データから生死を予測

コード

# from モジュール名 import クラス名（もしくは関数名や変数名）
import pandas as pd
from pandas import DataFrame
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.linear_model import LogisticRegression

# titanic data csvファイルの読み込み（パスは適宜変更）
titanic_df = pd.read_csv('/content/drive/My Drive/study_ai_ml_google/data/titanic_train.csv')

# 予測に不要と考えるカラムをドロップ
titanic_df.drop(['PassengerId', 'Name', 'Ticket', 'Cabin'], axis=1, inplace=True)

# Ageカラムのnullを中央値で補完
titanic_df['AgeFill'] = titanic_df['Age'].fillna(titanic_df['Age'].mean())

# 運賃だけのリストを作成
data1 = titanic_df.loc[:, ["Fare"]].values

#生死フラグのみのリストを作成
label1 = titanic_df.loc[:,["Survived"]].values

model = LogisticRegression()

model.fit(data1, label1)

# 予測
model.predict([[61]])

# 確率
model.predict_proba([[61]])

結果

array([0])
array([[0.50358033, 0.49641967]])

・運賃が61のとき死亡すると予測（確率： 0.50358033）

主成分分析

・次元圧縮する。ただし情報損失は少なくしたい
・そのために分散の大きな射影軸を見つける
・ノルムが1となる制約を入れる
→ラグランジュ関数を最大にする係数ベクトルを探索(微分して0になる点)

E(a_j) = a_j^T Var(\bar X) a_j - \lambda (a_j^T a_j - 1) \\
\frac{\partial E(a_j)}{\partial a_j} = 2 Var(\bar X) a_j - 2 \lambda a_j = 0 \\
Var(\bar X) a_j = \lambda a_j \\

・射影先の分散は固有値と一致

Var(s_1) = a_1^T Var(\bar X) a_1 = \lambda_1 a_1^T a_1 = \lambda_1

■寄与率
・第k主成分の分散の全分散に対する割合(第k主成分が持つ情報量の割合)
■累積寄与率
・第1-k主成分まで圧縮した際の全分散に対する割合

主成分分析・ハンズオン

コード

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
%matplotlib inline

cancer_df = pd.read_csv('CSVのパスを記入')

# 目的変数の抽出
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)

# 標準化
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)

# ロジスティック回帰で学習
logistic = LogisticRegressionCV(cv=10, random_state=0)
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('Confustion matrix:\n{}'.format(confusion_matrix(y_true=y_test, y_pred=logistic.predict(X_test_scaled))))

結果

Train score: 0.988
Test score: 0.972
Confustion matrix:
[[89  1]
 [ 3 50]]

・元の空間上でロジスティック回帰をすると約97%のテストデータを正確に分類できたことがわかる

コード

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

# PCA
# 次元数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))
# X_train_pca shape: (426, 2)

# 寄与率
evr = pca.explained_variance_ratio_
print('explained variance ratio: {}'.format(evr))
# explained variance ratio: [ 0.43315126  0.19586506]

# 累積寄与率
print('total explained variance ratio:', sum(evr))

# 散布図にプロット
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') # 第1主成分をx軸
plt.ylabel('PC 2') # 第2主成分をy軸

結果

explained variance ratio: [0.43315126 0.19586506]
total explained variance ratio: 0.6290163245532953

・2次元に圧縮すると約63%にまで累積寄与率が落ちてしまうが、人間の目で確認できる次元にまで落とし込むことができた

k近傍法(kNN)

・最近傍のデータを個取ってきて、それらがもっとも多く所属するクラスに識別する
・kを大きくすると決定境界は滑らかになる

k平均法(k-means)

・教師なし学習
・クラスタリング手法
・与えられたデータをk個のクラスタに分類する
・各クラスタ中心の初期値を設定し、各々の最近傍中心との距離から中心を再計算してそれを収束するまで繰り返す。

k-means++

・k-meansの初期値の設定に改良を加えたもの
・重みつき確率分布を用いて新しいクラスタ中心をランダムにk個選ぶことを繰り返す

k近傍法・ハンズオン

コード

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn import cluster, preprocessing, datasets
from sklearn.cluster import KMeans

wine = datasets.load_wine()
X = wine.data
X.shape
# (178, 13)
y = wine.target
y.shape
# (178,)
wine.target_names
# array(['class_0', 'class_1', 'class_2'], dtype='<U7')

model = KMeans(n_clusters=3)
labels = model.fit_predict(X)
df = pd.DataFrame({'labels': labels})
def species_label(theta):
    if theta == 0:
        return wine.target_names[0]
    if theta == 1:
        return wine.target_names[1]
    if theta == 2:
        return wine.target_names[2]
df['species'] = [species_label(theta) for theta in wine.target]
pd.crosstab(df['labels'], df['species'])

■結果

・初期値が近いとうまくクラスタリングできない場合があるため、改良版のk-means++などが考案されている

SVM（サポートベクターマシン）

・n次元空間上の集合をn-1次元の超平面で2クラスに分割することを考える
・この超平面と最も近いデータ（サポートベクトル）とのマージンを最大化する手法

■ソフトマージン
・誤差を許容し、その誤差にペナルティを与える

■非線形分離
・特徴空間に写像し、その空間で線形に分離する

SVM・ハンズオン

線形分離可能

訓練データ生成

%matplotlib inline
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)

学習・予測

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.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.quiver(0, 0, 0.1, 0.35, width=0.01, scale=1, color='pink')

線形分離不可

訓練データ生成

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)

学習・予測

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=['--', '-', '--'])

ソフトマージン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)

学習・予測

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=['--', '-', '--'])

・ソフトマージンを導入して誤判別を許容することで、重なり合うデータを分離することができた