支援向量機(SVM)和python實現(三)
阿新 • • 發佈:2018-12-08
6. python實現
根據前面的一步步推導獲得的結果,我們就可以使用python來實現SVM了 這裡我們使用iris資料集進行驗證,由於該資料集有4維,不容易在二維平面上表示,我們先使用LDA對其進行降維,又因為該資料集有3類樣本,我們編寫的SVM是二分類的,所以我們將獲取的第二個樣本的label設為1,其他兩類樣本的label設為-1
# -*- coding: gbk -*-
import numpy as np
import matplotlib.pyplot as plt
import math
from sklearn.datasets import load_iris
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
class SVM():
def __init__(self, C, kernel, kernel_arg, e=0.001):
'''
kernel_arg
kernel的型別: 'linear': 1
'poly': d(d>1且為整數)
'gaussian': σ(σ>0)
'lapras': σ(σ>0)
'sigmoid': beta,theta(beta>0,theta<0)
kernel_arg若不符合要求將按照預設引數進行計算
C為目標函式非線性部分的權重
e為誤差
'''
self.kernel = kernel
self.kernel_arg = kernel_arg
self.C = C
self.e = e
self.bias = 0
def kernel_function(self, x1, x2):
if self.kernel == 'linear':
return np.dot(x1, x2)
elif self.kernel == 'poly':
if isinstance(self.kernel_arg, int) == False :
self.kernel_arg = 2
return np.dot(x1, x2)**self.kernel_arg
elif self.kernel == 'gaussian':
if isinstance(self.kernel_arg, float) == False:
self.kernel_arg = 0.5
return math.exp(-np.linalg.norm(x1 - x2)**2 / (2 * self.kernel_arg**2))
elif self.kernel == 'lapras':
if isinstance(self.kernel_arg, float) == False:
self.kernel_arg = 0.5
return math.exp(-np.linalg.norm(x1 - x2) / self.kernel_arg)
elif self.kernel == 'sigmoid':
if len(self.kernel_arg) != 2:
self.kernel_arg = [0.5, -0.5]
if self.kernel_arg[0] <= 0:
self.kernel_arg[0] = 0.5
if self.kernel_arg[1] >= 0:
self.kernel_arg[1] = 0.5
return math.tanh(self.kernel_arg[0] * np.dot(x1, x2) + self.kernel_arg[1])
def fit(self, train_x, train_y, max_iter=1000):
self.train_x = np.array(train_x)
self.train_y = np.array(train_y)
self.alpha = np.zeros(train_x.shape[0])
iter = 0
while(iter < max_iter):
print('iter = {}'.format(iter))
index1, index2 = self.SMO_get_alpha()
if index1 == -1:
print('結束迭代, iter = {}'.format(iter))
break
train_result = self.SMO_train(index1, index2)
if train_result == True:
print('結束迭代, iter = {}'.format(iter))
break
iter += 1
def SMO_get_alpha(self):
for i in range(self.alpha.shape[0]):
if 0 < self.alpha[i] < self.C:
if self.train_y[i] * self.f(self.train_x[i]) != 1:
index2 = self.choose_another_alpha(i)
return i, index2
for i in range(self.alpha.shape[0]):
if self.alpha[i] == 0:
if self.train_y[i] * self.f(self.train_x[i]) < 1:
index2 = self.choose_another_alpha(i)
return i, index2
elif self.alpha[i] == self.C:
if self.train_y[i] * self.f(self.train_x[i]) > 1:
index2 = self.choose_another_alpha(i)
return i, index2
return -1, -1
def f(self, x):
result = 0
for i in range(self.alpha.shape[0]):
result += self.alpha[i] * self.train_y[i] * \
self.kernel_function(self.train_x[i], x)
return result + self.bias
def error(self, index):
return self.f(self.train_x[index]) - self.train_y[index]
def choose_another_alpha(self, index):
result_index = 0
temp_diff_error = 0
for i in range(self.alpha.shape[0]):
diff_error = np.abs(self.error(index) - self.error(i))
if diff_error > temp_diff_error:
temp_diff_error = diff_error
result_index = i
return result_index
def SMO_train(self, index1, index2):
old_alpha = self.alpha.copy()
x1 = self.train_x[index1]
y1 = self.train_y[index1]
x2 = self.train_x[index2]
y2 = self.train_y[index2]
eta = self.kernel_function(
x1, x1) + self.kernel_function(x2, x2) - 2 * self.kernel_function(x1, x2)
alpha2 = old_alpha[index2] + y2 * \
(self.error(index1) - self.error(index2)) / eta
if y1 != y2:
L = max(0, old_alpha[index2] - old_alpha[index1])
H = min(self.C, self.C + old_alpha[index2] - old_alpha[index1])
else:
L = max(0, old_alpha[index1] + old_alpha[index2] - self.C)
H = min(self.C, old_alpha[index1] + old_alpha[index2])
if alpha2 > H:
alpha2 = H
elif alpha2 < L:
alpha2 = L
alpha1 = old_alpha[index1] + y1 * y2 * (old_alpha[index2] - alpha2)
self.alpha[index1] = alpha1
self.alpha[index2] = alpha2
b1 = -self.error(index1) \
- y1 * self.kernel_function(x1, x1) * (alpha1 - old_alpha[index1]) \
- y2 * self.kernel_function(x1, x2) * (alpha2 - old_alpha[index2]) \
+ self.bias
b2 = -self.error(index2) \
- y1 * self.kernel_function(x1, x2) * (alpha1 - old_alpha[index1]) \
- y2 * self.kernel_function(x2, x2) * (alpha2 - old_alpha[index2]) \
+ self.bias
if 0 < alpha1 < self.C:
self.bias = b1
elif 0 < alpha2 < self.C:
self.bias = b2
else:
self.bias = (b1 + b2) / 2
print('E = {}'.format(np.linalg.norm(old_alpha - self.alpha)))
if np.linalg.norm(old_alpha - self.alpha) < self.e:
return True
else:
return False
def predict_one(self, x):
if self.f(x) > 0:
return 1
else:
return -1
def predict(self, x_group):
return np.array([self.predict_one(x) for x in x_group])
if __name__ == '__main__':
iris = load_iris()
X = iris.data
y = iris.target
lda = LinearDiscriminantAnalysis(n_components=2)
lda.fit(X, y)
X = lda.transform(X)
train_data = X[0::2, :]
train_label = y[0::2]
for i in range(train_label.shape[0]):
if train_label[i] != 1:
train_label[i] = -1
else:
train_label[i] = 1
test_data = X[1::2, :]
test_label = y[1::2]
for i in range(test_label.shape[0]):
if test_label[i] != 1:
test_label[i] = -1
else:
test_label[i] = 1
svm = SVM(5, 'gaussian', 0.5, 0.001)
svm.fit(train_data, train_label)
predict_label = svm.predict(test_data)
a = predict_label - test_label
print(a)
# print(svm.alpha)
count = 0
for i in range(a.shape[0]):
if a[i] == 0:
count += 1
print(count / test_label.shape[0])
points = []
for i in np.linspace(-10.0, 10.0, num=400):#獲取超平面上的點為了作圖
for j in np.linspace(-5.0, 5.0, num=200):
x_ij = np.array([i, j])
if -0.05 < svm.f(x_ij) < 0.05:
tmp = [i, j]
points.append(tmp)
points = np.array(points)
print(points)
plt.scatter(X[:, 0], X[:, 1], marker='o', c=y)
plt.scatter(points[:, 0], points[:, 1], marker='o')
plt.show()
我們用了高斯核獲取樣本和劃分曲線,由於不知道怎麼畫出超平面,我只好取巧不斷將點代入方程,判斷結果是否與0很接近,然後把和0接近的點畫出來,另外我發現在用高斯核的時候e設0.01就好了,但是在用poly核的時候需要設的小很多,不然得到的結果很差。