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核的時候需要設的小很多，不然得到的結果很差。

參考：