import numpy as np
import matplotlib.pyplot as plt
import scipy
from scipy import stats
from numpy.random import randn # 正規乱数
import cvxopt
from cvxopt import matrix
%matplotlib inline
# import japanize_matplotlib # Google Colabの場合
# Anacondaの場合は下記
import matplotlib
from matplotlib import font_manager
matplotlib.rc("font", family="BIZ UDGothic")
a = randn(1)
b = randn(1)
n = 100
X = randn(n, 2)
y = np.sign(a*X[:, 0] + b*X[:, 1] + 0.1*randn(n))
y = y.reshape(-1, 1) # 形を明示してわたす必要がある
def svm_1(X, y, C):
eps = 0.0001
n = X.shape[0]
P = np.zeros((n, n))
for i in range(n):
for j in range(n):
P[i, j] = np.dot(X[i, :], X[j, :]) * y[i] * y[j]
# パッケージにあるmatrix関数を使って指定する必要がある
P = matrix(P + np.eye(n)*eps)
A = matrix(- y.T.astype(np.float))
b = matrix(np.array([0]).astype(np.float))
h = matrix(np.array([C]*n + [0]*n).reshape(-1, 1).astype(np.float))
G = matrix(np.concatenate([np.diag(np.ones(n)), np.diag(-np.ones(n))]))
q = matrix(np.array([-1]*n).astype(np.float))
res = cvxopt.solvers.qp(P, q, A=A, b=b, G=G, h=h) # ソルバーの実行
alpha = np.array(res["x"]) # xが本文中のalphaに対応
beta = ((alpha*y).T@X).reshape(2, 1)
index = (eps < alpha[:, 0]) & (alpha[:, 0] < C - eps)
beta_0 = np.mean(y[index] - X[index, :]@beta)
return {"beta":beta, "beta_0":beta_0}
def f(x):
return - res["beta_0"]/res["beta"][1] - x*res["beta"][0]/res["beta"][1]
a = randn(1)
b = randn(1)
n = 100
X = randn(n, 2)
y = np.sign(a*X[:, 0] + b*X[:, 1] + 0.1*randn(n))
y = y.reshape(-1, 1) # 形を明示してわたす必要がある
for i in range(n):
if y[i] == 1:
plt.scatter(X[i, 0], X[i, 1], c="red")
else :
plt.scatter(X[i, 0], X[i, 1], c="blue")
res = svm_1(X, y, C=10)
x_seq = np.arange(-3, 3, 0.5)
plt.plot(x_seq, f(x_seq))
def K_linear(x, y):
return x.T@y
# ちなみに多項式カーネルは以下のように構成する
def K_poly(x, y):
return (1 + x.T@y) ** 2
# 実行
a = 3
b = -1
n = 200
X = randn(n, 2)
y = np.sign(a*X[:, 0] + b*X[:, 1]**2 + 0.3*randn(n))
y = y.reshape(-1, 1)
def plot_kernel(K, line): # 引数lineで線の種類を指定する
res = svm_2(X, y, 1, K)
alpha = res["alpha"][:, 0]
beta_0 = res["beta_0"]
def f(u, v):
S = beta_0
for i in range(X.shape[0]):
S = S + alpha[i]*y[i]*K(X[i, :], [u, v])
return S[0]
uu = np.arange(-2, 2, 0.1)
vv = np.arange(-2, 2, 0.1)
ww = []
for v in vv:
w = []
for u in uu:
w.append(f(u, v))
ww.append(w)
plt.contour(uu, vv, ww, levels=0, linestyles=line)
def svm_2(X, y, C, K):
eps = 0.0001
n = X.shape[0]
P = np.zeros((n, n))
for i in range(n):
for j in range(n):
P[i, j] = K(X[i, :], X[j, :])*y[i]*y[j]
# パッケージにあるmatrix関数を使って指定する必要がある
P = matrix(P + np.eye(n)*eps)
A = matrix(-y.T.astype(np.float))
b = matrix(np.array([0]).astype(np.float))
h = matrix(np.array([C]*n + [0]*n).reshape(-1, 1).astype(np.float))
G = matrix(np.concatenate([np.diag(np.ones(n)), np.diag(-np.ones(n))]))
q = matrix(np.array([-1]*n).astype(np.float))
res = cvxopt.solvers.qp(P, q, A=A, b=b, G=G, h=h)
alpha = np.array(res["x"]) # xが本文中のalphaに対応
beta = ((alpha*y).T@X).reshape(2, 1)
index = (eps < alpha[:, 0]) & (alpha[:, 0] < C - eps)
beta_0 = np.mean(y[index] - X[index, :]@beta)
return {"alpha":alpha, "beta":beta, "beta_0":beta_0}
uu = np.arange(-2, 2, 0.1)
vv = np.arange(-2, 2, 0.1)
ww = []
for v in vv:
w = []
for u in uu:
w.append(f(u, v))
ww.append(w)
plt.contour(uu, vv, ww, levels=0, linestyles=line)
for i in range(n):
if y[i] == 1:
plt.scatter(X[i, 0], X[i, 1], c="red")
else:
plt.scatter(X[i, 0], X[i, 1], c="blue")
plot_kernel(K_poly, line="dashed")
plot_kernel(K_linear, line="solid")
import sklearn
from sklearn import svm
x = randn(200, 2)
x[0:100, ] = x[0:100, ] + 2
x[100:150, ] = x[100:150, ] - 2
y = np.concatenate(([1 for i in range(150)], [2 for i in range(50)]))
train = np.random.choice(200, 100, replace=False)
test = list(set(range(200)) - set(train))
res_svm = svm.SVC(kernel="rbf", gamma=1, C=100) # チューニングなしのSVM
res_svm.fit(x[train, ], y[train]) # 実行
res_svm.predict(x[test, ]) # テストデータの予測結果
import mlxtend
from mlxtend.plotting import plot_decision_regions
plot_decision_regions(x, y, clf=res_svm)
from sklearn.model_selection import GridSearchCV
grid = {"C": [0.1, 1, 10, 100, 1000], "gamma" : [0.5, 1, 2, 3, 4]}
tune = GridSearchCV(svm.SVC(), grid, cv=10)
tune.fit(x[train, ], y[train])
from sklearn.datasets import load_iris
iris = load_iris()
iris.target_names
x = iris.data
y = iris.target
train = np.random.choice(150, 120, replace=False)
test = np.ones(150, dtype=bool)
test[train] = False
iris_svm = svm.SVC(kernel="rbf", gamma=1, C=10)
iris_svm.fit(x[train, ], y[train])
# 2章で定義した関数table_count (再掲)
def table_count(m, u, v):
n = u.shape[0]
count = np.zeros([m, m])
for i in range(n):
count[int(u[i]), int(v[i])] += 1
return (count)
y_pre = iris_svm.predict(x[test, ])
table_count(3, y[test], y_pre)