支持向量機核函數的實現

一：回顧SVM中的SMO算法

https://www.cnblogs.com/ssyfj/p/13363526.html

二：核函數的了解

（一）西瓜書（粗略了解）

（二）統計學習方法（詳細）

（三）推文：支持向量機原理(三)線性不可分支持向量機與核函數

（四）推文：核函數和核矩陣

（五）機器學習實戰（代碼加強理解）

三：重點知識

對於在低維線性不可分的數據，在映射到了高維以后，就變成線性可分。這個映射過程分為隱式轉換和顯示轉換兩種情況。

（一）顯式映射（x->顯示構造高維特征φ(x)->內積計算φ(x).T@φ(x))

線性不可分的低維特征數據，我們可以將其映射到高維，就能線性可分。現在我們將它運用到我們的SVM的算法上。回顧線性可分SVM的優化目標函數：

注意到上式低維特征僅僅以內積xi∙xj的形式出現，如果我們定義一個低維特征空間到高維特征空間的映射ϕ(x)（比如從2維到5維的映射），將所有特征映射到一個更高的維度，讓數據線性可分，我們就可以繼續按前面的SMO方法來優化目標函數，求出分離超平面和分類決策函數了。

也就是說現在的SVM的優化目標函數變成：

可以看到，和線性可分SVM的優化目標函數的區別僅僅是將內積x_i∙x_j替換為ϕ(x_i)∙ϕ(x_j)。

這就是顯式映射，我們需要真正意義上將數據映射到高維特征空間中去，然后對這個轉換后的高維向量進行內積運算。

看起來似乎這樣我們就已經完美解決了線性不可分SVM的問題了，但是事實是不是這樣呢？

我們看看，假如是一個2維特征的數據，我們可以將其映射到5維來做特征的內積，如果原始空間是三維，可以映射到到19維空間，似乎還可以處理。

但是如果我們的低維特征是100個維度，1000個維度呢？

那么我們要將其映射到超級高的維度來計算特征的內積。

這時候映射成的高維維度是爆炸性增長的，這個計算量實在是太大了，而且如果遇到無窮維的情況，就根本無從計算了。

顯示轉換存在以下兩個問題：

1.通常我們不知道映射ϕ的具體形式； 2.映射后的空間維數可能非常高，甚至無限維，直接計算開銷太大，十分困難.

（二）隱式映射（核函數）---見統計學習方法p136

利用核函數，不需要將輸入空間中的樣本映射到新的空間中去，在輸入空間中就可以直接計算出內積了。

對輸入空間向高維空間的一種隱式映射（注意，這也是低維空間到高維空間的一種映射<理解意義上的>）。不需要給出那個映射，而在輸入空間中就可以計算內積，即核函數方法。【核函數可以簡化映射空間中的內積運算】

（三）顯式映射與隱式映射案例對比（加深理解）---見西瓜書筆記

顯示映射：需要構建高維空間，需要大量計算

隱式映射：使用核函數，直接將X作為輸入，避免中間構造高維特征空間的計算。所有的計算都是在原始特征空間中進行的，並沒有構建高維空間。

（四）核技巧---統計學習方法

在學習核預測中只定義核函數K(x,z),而不是顯式的定義映射函數。因為直接計算K(x,z)比較容易，而通過φ(x)核φ(z)計算K(x,z)並不容易（也就是說明我們即便是想要通過顯式映射到高維特征空間進行划分數據集並不容易）：

因為對於給定的核K(x,z)或者我們想要的高維中划分超平面選取，高維特征空間和映射函數的取法並不唯一，可以選取不同的特征空間，即便是在同一特征空間里，也可以取不同的映射

四：核函數代碼實現

（一）實現核函數轉換函數

#實現核函數（核轉換函數）
def kernelTrans(data_X,data_Xi,kTup):
    """
 因為我們的核函數是對數據集進行的隱式映射，后面求解α等，並沒有對數據集本身修改過， 所以我們只需要調用一次核轉換函數即可求得全局的核矩陣解 :param data_X: 表示傳入的全部數據集 :param data_Xi: 表示傳入的第i個樣本數據 :param kTup: 是一個包含核函數信息的元組 :return: 返回核矩陣的第i列 """
    m,n = data_X.shape
    K = np.zeros((m,1)) #這就是我們要返回的核矩陣的第i列初始狀態

    #根據我們輸入的核函數信息判斷選取的是哪一個核函數，並獲取傳入參數
    if kTup[0] == 'lin':    #線性核函數，即沒有使用核函數
        K = data_X@data_Xi  #(m,1)
    elif kTup[0] == 'rbf':     #高斯核函數
        for j in range(m):  #對數據集每一個樣本進行高斯處理
            deltaRow = data_X[j,:] - data_Xi    #獲取差值
 K[j] = deltaRow@deltaRow.T #獲取高斯函數分子部分，無論上面求解的deltaRow正負，這里內積以后，全是正值
        K = np.exp(K/(-1*kTup[1]**2))   #對整體K進行除方差操作，獲取高斯處理后的值
    else:   #如果還有其他的核函數，接着使用elif即可
        raise NameError("the kernel is not recognized")#對於沒有的核函數，我們拋出異常

    return K

（二）修改數據結構，獲取全局核矩陣

#首先我們定義一個數據結構（類），來保存所有的重要值
class optStruct:
    def __init__(self,data_X,data_Y,C,toler,kTup):   #輸入參數分別是數據集、類別標簽、常數C用於軟間隔、和容錯率toler
        self.X = data_X
        self.label = data_Y
        self.C = C
        self.toler = toler  #就是軟間隔中的ε，調節最大間隔大小
        self.m = data_X.shape[0]
        self.alphas = np.zeros(self.m)    #存放每個樣本點的α值
        self.b = 0  #存放閾值
        self.eCache = np.zeros((self.m,2))  #用於緩存誤差，每個樣本點對應一個Ei值，第一列為標識符，標志是否為有效值，第二列存放有效值
        self.K = np.zeros((self.m,self.m))
        #這里獲取全局核矩陣即可
        for i in range(self.m): self.K[:,i] = kernelTrans(data_X,data_X[i],kTup).flatten()

注意：核矩陣是對稱陣，由（一）可以知道，我們求解的最終K核矩陣是對稱陣，即半正定矩陣（實對稱矩陣A稱為半正定矩陣），所以我們是符合核矩陣要求的。

（三）修改內循環函數---因為這里面使用了核函數

#三：實現內循環函數，相比於外循環，這里包含了主要的更新操作
def innerL(i,oS):   #由外循環提供i值（具體選取要違背kkT<這里實現>,使用交替遍歷<外循環中實現>）---提供α_1的索引
    Ei = calcEk(oS,i)   #計算E1值，主要是為了下面KKT條件需要使用到

    #如果下面違背了KKT條件，則正常進行α、Ek、b的更新,重點：后面單獨說明下面是否滿足違反KKT條件
    if ((oS.label[i]*Ei < -oS.toler) and (oS.alphas[i] < oS.C)) or\
       ((oS.label[i]*Ei > oS.toler) and (oS.alphas[i] > 0)):  #注意：對於硬間隔，我們直接和1對比，對於軟間隔，我們要和1 +或- ε對比
        #開始在內循環中，選取差值最大的α_2下標索引
        j,Ej = selectJ(i,oS,Ei)
        #因為后面要修改α_1與α_2的值，但是后面修改閾值b的時候需要用到新舊兩個值，所以我們需要在更新α值之前進行保存舊值
        alphaIold = oS.alphas[i].copy()
        alphaJold = oS.alphas[j].copy()

        #分析約束條件（是對所有α都適用），一會對我們新的α_2進行截取糾正，注意：α_1是由α_2推出的，所以不需要進行驗證了。
        #如果y_1!=y_2異號時：
        if oS.label[i] != oS.label[j]:
            L = max(0,alphaJold-alphaIold)
            H = min(oS.C,oS.C+alphaJold-alphaIold)
        else:   #如果y_1==y_2同號時
            L = max(0,alphaJold+alphaIold-oS.C)
            H = min(oS.C,alphaJold+alphaIold)
        #上面就是將α_j調整到L,H之間
        if L==H:    #如果L==H，之間返回0，跳出這次循環，不進行改變(單值選擇，沒必要）
            return 0

        #計算η值=k_11+k_22-2k_12
        eta = oS.K[i,i] + oS.K[j,j] - 2.0*oS.K[i,j] #eta性質可以知道是>=0的，所以我們只需要判斷是否為0即可 if eta <= 0:
            print("eta <= 0")
            return 0

        #當上面所有條件都滿足以后，我們開始正式修改α_2值，並更新對應的Ek值
        oS.alphas[j] += oS.label[j]*(Ei-Ej)/eta
        oS.alphas[j] = clipAlpha(oS.alphas[j],H,L)
        updateEk(oS,j)

        #查看α_2是否有足夠的變化量，如果沒有足夠變化量，我們直接返回，不進行下面更新α_1,注意：因為α_2變化量較小，所以我們沒有必要非得把值變回原來的舊值
        if abs(oS.alphas[j] - alphaJold) < 0.00001:
            print("J not move enough")
            return 0

        #開始更新α_1值,和Ek值
        oS.alphas[i] += oS.label[i]*oS.label[j]*(alphaJold-oS.alphas[j])
        updateEk(oS,i)

 #開始更新閾值b,正好使用到了上面更新的Ek值 b1 = oS.b - Ei - oS.label[i] * (oS.alphas[i] - alphaIold) * oS.K[i,i] - oS.label[j] * ( oS.alphas[j] - alphaJold) * oS.K[i,j] b2 = oS.b - Ej - oS.label[i] * (oS.alphas[i] - alphaIold) * oS.K[i,j] - oS.label[j] * ( oS.alphas[j] - alphaJold) * oS.K[j,j]

        #根據統計學習方法中閾值b在每一步中都會進行更新，
        #1.當新值alpha_1不在界上時(0<alpha_1<C)，b_new的計算規則為：b_new=b1
        #2.當新值alpha_2不在界上時(0 < alpha_2 < C)，b_new的計算規則為：b_new = b2
        #3.否則當alpha_1和alpha_2都不在界上時，b_new = 1/2(b1+b2)
        if oS.alphas[i] > 0 and oS.alphas[i] < oS.C:
            oS.b = b1
        elif oS.alphas[j] > 0 and oS.alphas[j] < oS.C:
            oS.b = b2
        else:
            oS.b = 1/2*(b1+b2)

        #注意：這里我們應該根據b_new更新一次Ei，但是我們這里沒有寫，因為我們將這一步提前到了最開始，即selectJ中

        #以上全部更新完畢，開始返回標識
        return 1
    return 0    #沒有違背KKT條件

（四）修改計算誤差函數---含核函數

#計算每個樣本點k的Ek值，就是計算誤差值=預測值-標簽值
def calcEk(oS,k):
    # 根據西瓜書6.24，我們可以知道預測值如何使用α值進行求解
    fxk = np.multiply(oS.alphas,oS.label).T@(oS.K[:,k])+oS.b #np.multiply之后還是(m,1),(oS.X@oS.X[k,:])之后是(m,1),通過轉置(1,m)@(m,1)-->實數后+b即可得到預測值fx  #獲取誤差值Ek
    Ek = fxk - oS.label[k]
    return Ek

（五）全部代碼

import numpy as np
import matplotlib.pyplot as plt
from numpy import *

#一：SMO算法中的輔助函數
#加載數據集
def loadDataSet(filename):
    dataSet = np.loadtxt(filename)
    m,n = dataSet.shape
    data_X = dataSet[:,0:n-1]
    data_Y = dataSet[:,n-1]

    return data_X,data_Y

#隨機選取一個數J，為后面內循環選取α_2做輔助（如果α選取不滿足條件，就選擇這個方法隨機選取）
def selectJrand(i,m):   #主要就是根據α_1的索引i,從所有數據集索引中隨機選取一個作為α_2的索引
    j = i
    while j==i:
        j = np.int(np.random.uniform(0,m))  #從0~m中隨機選取一個數，是進行整數化的
    print("random choose index for α_2:%d"%(j))
    return j    #由於這里返回隨機數，所以后面結果 可能導致不同

def clipAlpha(aj,H,L):  #根據我們的SVM算法中的約束條件的分析，我們對獲取的aj,進行了截取操作
    if aj > H:
        aj = H
    if aj < L:
        aj = L
    return aj

#二：SMO的支持函數
#實現核函數（核轉換函數）
def kernelTrans(data_X,data_Xi,kTup):
    """
    因為我們的核函數是對數據集進行的隱式映射，后面求解α等，並沒有對數據集本身修改過，
    所以我們只需要調用一次核轉換函數即可求得全局的核矩陣解
    :param data_X:  表示傳入的全部數據集
    :param data_Xi: 表示傳入的第i個樣本數據
    :param kTup: 是一個包含核函數信息的元組
    :return: 返回核矩陣的第i列
    """
    m,n = data_X.shape
    K = np.zeros((m,1)) #這就是我們要返回的核矩陣的第i列初始狀態

    #根據我們輸入的核函數信息判斷選取的是哪一個核函數，並獲取傳入參數
    if kTup[0] == 'lin':    #線性核函數，即沒有使用核函數
        K = data_X@data_Xi  #(m,1)
    elif kTup[0] == 'rbf':     #高斯核函數
        for j in range(m):  #對數據集每一個樣本進行高斯處理
            deltaRow = data_X[j,:] - data_Xi    #獲取差值
            K[j] = deltaRow@deltaRow.T  #獲取高斯函數分子部分
        K = np.exp(K/(-1*kTup[1]**2))   #對整體K進行除方差操作，獲取高斯處理后的值
    else:   #如果還有其他的核函數，接着使用elif即可
        raise NameError("the kernel is not recognized")#對於沒有的核函數，我們拋出異常

    return K

#首先我們定義一個數據結構（類），來保存所有的重要值
class optStruct:
    def __init__(self,data_X,data_Y,C,toler,kTup):   #輸入參數分別是數據集、類別標簽、常數C用於軟間隔、和容錯率toler
        self.X = data_X
        self.label = data_Y
        self.C = C
        self.toler = toler  #就是軟間隔中的ε，調節最大間隔大小
        self.m = data_X.shape[0]
        self.alphas = np.zeros(self.m)    #存放每個樣本點的α值
        self.b = 0  #存放閾值
        self.eCache = np.zeros((self.m,2))  #用於緩存誤差，每個樣本點對應一個Ei值，第一列為標識符，標志是否為有效值，第二列存放有效值
        self.K = np.zeros((self.m,self.m))
        #這里獲取全局核矩陣即可
        for i in range(self.m):
            self.K[:,i] = kernelTrans(data_X,data_X[i],kTup).flatten()

#計算每個樣本點k的Ek值，就是計算誤差值=預測值-標簽值
def calcEk(oS,k):
    # 根據西瓜書6.24，我們可以知道預測值如何使用α值進行求解
    fxk = np.multiply(oS.alphas,oS.label).T@(oS.K[:,k])+oS.b #np.multiply之后還是(m,1),(oS.X@oS.X[k,:])之后是(m,1),通過轉置(1,m)@(m,1)-->實數后+b即可得到預測值fx
    #獲取誤差值Ek
    Ek = fxk - oS.label[k]
    return Ek

#內循環的啟發式方法，獲取最大差值|Ei-Ej|對應的Ej的索引J
def selectJ(i,oS,Ei):   #注意我們要傳入第一個α對應的索引i和誤差值Ei,后面會用到
    maxK = -1   #用於保存臨時最大索引
    maxDeltaE = 0   #用於保存臨時最大差值--->|Ei-Ej|
    Ej = 0  #保存我們需要的Ej誤差值

    #重點：這里我們是把SMO最后一步（根據最新閾值b，來更新Ei)提到第一步來進行了，所以這一步是非常重要的
    oS.eCache[i] = [1,Ei]

    #開始獲取各個Ek值，比較|Ei-Ej|獲取Ej的所有
    #獲取所有有效的Ek值對應的索引
    validECacheList = np.where(oS.eCache[:,0]!=0)[0]    #根據誤差緩存中第一列非0，獲取對應的有效誤差值
    if len(validECacheList) > 1:   #如果有效誤差緩存長度大於1(因為包括Ei),則正常進行獲取j值，否則使用selectJradn方法選取一個隨機J值
        for k in validECacheList:
            if k == i:  #相同則不處理
                continue
            #開始計算Ek值，進行對比，獲取最大差值
            Ek = calcEk(oS,k)
            deltaE = abs(Ei - Ek)
            if deltaE > maxDeltaE:  #更新Ej及其索引位置
                maxK = k
                maxDeltaE = deltaE
                Ej = Ek
        return maxK,Ej  #返回我們找到的第二個變量α_2的位置
    else:   #沒有有效誤差緩存，則隨機選取一個索引，進行返回
        j = selectJrand(i,oS.m)
        Ej = calcEk(oS,j)
        return j,Ej

#實現更新Ek操作,因為除了最后我們需要更新Ei之外，我們在內循環中計算α_1與α_2時還是需要用到E1與E2，
#因為每次的E1與E2由於上一次循環中更新了α值，所以這一次也是需要更新E1與E2值，所以單獨實現一個更新Ek值的方法還是有必要的
def updateEk(oS,k):
    Ek = calcEk(oS,k)
    oS.eCache[k] = [1,Ek]   #第一列1，表示為有效標識

#三：實現內循環函數，相比於外循環，這里包含了主要的更新操作
def innerL(i,oS):   #由外循環提供i值（具體選取要違背kkT<這里實現>,使用交替遍歷<外循環中實現>）---提供α_1的索引
    Ei = calcEk(oS,i)   #計算E1值，主要是為了下面KKT條件需要使用到

    #如果下面違背了KKT條件，則正常進行α、Ek、b的更新,重點：后面單獨說明下面是否滿足違反KKT條件
    if ((oS.label[i]*Ei < -oS.toler) and (oS.alphas[i] < oS.C)) or\
       ((oS.label[i]*Ei > oS.toler) and (oS.alphas[i] > 0)):  #注意：對於硬間隔，我們直接和1對比，對於軟間隔，我們要和1 +或- ε對比
        #開始在內循環中，選取差值最大的α_2下標索引
        j,Ej = selectJ(i,oS,Ei)
        #因為后面要修改α_1與α_2的值，但是后面修改閾值b的時候需要用到新舊兩個值，所以我們需要在更新α值之前進行保存舊值
        alphaIold = oS.alphas[i].copy()
        alphaJold = oS.alphas[j].copy()

        #分析約束條件（是對所有α都適用），一會對我們新的α_2進行截取糾正，注意：α_1是由α_2推出的，所以不需要進行驗證了。
        #如果y_1!=y_2異號時：
        if oS.label[i] != oS.label[j]:
            L = max(0,alphaJold-alphaIold)
            H = min(oS.C,oS.C+alphaJold-alphaIold)
        else:   #如果y_1==y_2同號時
            L = max(0,alphaJold+alphaIold-oS.C)
            H = min(oS.C,alphaJold+alphaIold)
        #上面就是將α_j調整到L,H之間
        if L==H:    #如果L==H，之間返回0，跳出這次循環，不進行改變(單值選擇，沒必要）
            return 0

        #計算η值=k_11+k_22-2k_12
        eta = oS.K[i,i] + oS.K[j,j] - 2.0*oS.K[i,j]  #eta性質可以知道是>=0的，所以我們只需要判斷是否為0即可
        if eta <= 0:
            print("eta <= 0")
            return 0

        #當上面所有條件都滿足以后，我們開始正式修改α_2值，並更新對應的Ek值
        oS.alphas[j] += oS.label[j]*(Ei-Ej)/eta
        oS.alphas[j] = clipAlpha(oS.alphas[j],H,L)
        updateEk(oS,j)

        #查看α_2是否有足夠的變化量，如果沒有足夠變化量，我們直接返回，不進行下面更新α_1,注意：因為α_2變化量較小，所以我們沒有必要非得把值變回原來的舊值
        if abs(oS.alphas[j] - alphaJold) < 0.00001:
            print("J not move enough")
            return 0

        #開始更新α_1值,和Ek值
        oS.alphas[i] += oS.label[i]*oS.label[j]*(alphaJold-oS.alphas[j])
        updateEk(oS,i)

        #開始更新閾值b,正好使用到了上面更新的Ek值
        b1 = oS.b - Ei - oS.label[i] * (oS.alphas[i] - alphaIold) * oS.K[i,i] - oS.label[j] * (
                    oS.alphas[j] - alphaJold) * oS.K[i,j]

        b2 = oS.b - Ej - oS.label[i] * (oS.alphas[i] - alphaIold) * oS.K[i,j] - oS.label[j] * (
                    oS.alphas[j] - alphaJold) * oS.K[j,j]

        #根據統計學習方法中閾值b在每一步中都會進行更新，
        #1.當新值alpha_1不在界上時(0<alpha_1<C)，b_new的計算規則為：b_new=b1
        #2.當新值alpha_2不在界上時(0 < alpha_2 < C)，b_new的計算規則為：b_new = b2
        #3.否則當alpha_1和alpha_2都不在界上時，b_new = 1/2(b1+b2)
        if oS.alphas[i] > 0 and oS.alphas[i] < oS.C:
            oS.b = b1
        elif oS.alphas[j] > 0 and oS.alphas[j] < oS.C:
            oS.b = b2
        else:
            oS.b = 1/2*(b1+b2)

        #注意：這里我們應該根據b_new更新一次Ei，但是我們這里沒有寫，因為我們將這一步提前到了最開始，即selectJ中

        #以上全部更新完畢，開始返回標識
        return 1
    return 0    #沒有違背KKT條件

#四：開始外循環，由於我們在內循環中實現了KKT條件的判斷，所以這里我們只需要進行交替遍歷即可
#交替遍歷一種方式是在所有的數據集上進行單遍掃描，另一種是在非邊界上（不在邊界0或C上的值）進行單遍掃描
# 交替遍歷：
# 交替是通過一個外循環來選擇第一個alpha值的，並且其選擇過程會在兩種方式之間交替：
# 一種方式是在所有數據集上進行單遍掃描，
# 另一種方式則是在非邊界alpha中實現單遍掃描，所謂非邊界alpha指的是那些不等於邊界0或C的alpha值。
# 對整個數據集的掃描相當容易，
# 而實現非邊界alpha值的掃描時，首先需要建立這些alpha值的列表，然后對這個表進行遍歷。
# 同時，該步驟會跳過那些已知不變的alpha值。
def smoP(data_X,data_Y,C,toler,maxIter,kTup):
    oS = optStruct(data_X,data_Y,C,toler,kTup)
    iter = 0
    entireSet = True    #標志是否應該遍歷整個數據集
    alphaPairsChanged = 0   #標志一次循環中α更新的次數
    #開始進行迭代
    #當iter >= maxIter或者((alphaPairsChanged == 0) and not entireSet)退出循環
    #前半個判斷條件很好理解，后面的判斷條件中，表示上一次循環中，是在整個數據集中遍歷，並且沒有α值更新過，則退出
    while iter < maxIter and ((alphaPairsChanged > 0) or entireSet):
        alphaPairsChanged = 0
        if entireSet:   #entireSet是true，則在整個數據集上進行遍歷
            for i in range(oS.m):
                alphaPairsChanged += innerL(i,oS)   #調用內循環
            print("full dataset, iter: %d i:%d,pairs changed:%d"%(iter,i,alphaPairsChanged))
            iter += 1   #無論是否更新過，我們都計算迭代一次
        else:   #遍歷非邊界值
            nonBounds = np.where((oS.alphas>0) & (oS.alphas<C))[0]  #獲取非邊界值中的索引
            for i in nonBounds: #開始遍歷
                alphaPairsChanged += innerL(i,oS)
            print("non bound, iter: %d i:%d,pairs changed:%d"%(iter,i,alphaPairsChanged))
            iter += 1   #無論是否更新過，我們都計算迭代一次

        #下面實現交替遍歷
        if entireSet:
            entireSet = False
        elif alphaPairsChanged == 0:    #如果是在非邊界上，並且α更新過。則entireSet還是False,下一次還是在非邊界上進行遍歷。可以認為這里是傾向於非邊界遍歷，因為非邊界遍歷的樣本更符合內循環中的違反KKT條件
            entireSet = True

        print("iteration number: %d"%iter)

    return oS.b,oS.alphas

def calcWs(alphas,data_X,data_Y):
    #根據西瓜書6.37求W
    m,n = data_X.shape
    w = np.zeros(n)
    for i in range(m):
        w += alphas[i]*data_Y[i]*data_X[i].T

    return w

#繪制圖像
def plotFigure(data_X,data_Y,alphas):
    m,n = data_X.shape
    # 進行數據集分類操作
    cls_1x = data_X[np.where(data_Y==1)]
    cls_2x = data_X[np.where(data_Y!=1)]

    plt.scatter(cls_1x[:,0].flatten(), cls_1x[:,1].flatten(), s=30, c='b')
    plt.scatter(cls_2x[:,0].flatten(), cls_2x[:,1].flatten(), s=30, c='g', marker='s')

    # 畫出支持向量點
    for i in range(m):
        if alphas[i] > 0.0:
            plt.scatter(data_X[i, 0], data_X[i, 1], s=150, c='none', alpha=0.7, linewidth=1.5, edgecolor='r')

    plt.xlim((-1.5, 1.5))
    plt.ylim((-1, 1.5))
    plt.show()

def testRbf(data_X,data_Y,C,toler,maxIter,kTup):
    #獲取閾值b和α值
    b,alphas = smoP(data_X,data_Y,C,toler,maxIter,kTup)
    plotFigure(data_X,data_Y,alphas)
    # print(b)
    # print(alphas)
    #獲取支持向量，即α>0點
    svIdx = np.where(alphas>0)[0]
    sVs = data_X[svIdx]
    labelSV = data_Y[svIdx]
    print("there are %d support vectors"%sVs.shape[0])
    m,n = data_X.shape
    errorCount = 0

    #獲取數據異常錯誤率
    for i in range(m):
        #先進行預測，預測公式見西瓜書6.24
        kernelEval = kernelTrans(sVs,data_X[i],kTup)    #使用支持向量點，和對應的預測數據點，進行預測
        pred = [e-mail](labelSV,alphas[svIdx]) + b #進行預測
        if np.sign(pred) != np.sign(data_Y[i]):
            errorCount += 1
    print("the trainning set error rate is: %f"%(errorCount/m))

    #測試集預測
    data_X,data_Y = loadDataSet("testSetRBF2.txt")
    m,n = data_X.shape
    errorCount = 0

    #獲取數據異常錯誤率
    for i in range(m):
        #先進行預測，預測公式見西瓜書6.24
        kernelEval = kernelTrans(sVs,data_X[i],kTup)    #使用支持向量點，和對應的預測數據點，進行預測
        pred = [e-mail](labelSV,alphas[svIdx]) + b #進行預測
        if np.sign(pred) != np.sign(data_Y[i]):
            errorCount += 1
    print("the testing set error rate is: %f"%(errorCount/m))

data_X,data_Y = loadDataSet("testSetRBF.txt")
C = 200
toler = 0.00001
maxIter = 10000
kTup = ('rbf',1.3)

testRbf(data_X,data_Y,C,toler,maxIter,kTup)

（六）結果顯示

kTup=('rbf',1.3)時：

 

kTup=('rbf',0.1)時：

 

五：手寫數字識別問題

（一）與KNN對比

與KNN算法對比，KNN算法效果不錯，但是想要保留所有的訓練樣本，在預測的時候也是需要所有的樣本，但是我們使用SVM時，我們只需要保留支持向量即可，保留的樣本少了很多，效果不會太差。

注意:我們主要實現了SVM對數字識別的二分類問題,主要對數字1和9進行了分類

（二）代碼實現獲取數據集

#將每一個數字文件轉換為矩陣向量
def img2vector(filename):
    data = []
    with codecs.open(filename,'r') as fp:
        for i in range(32):
            linestr = fp.readline() #讀取一行數據
            for j in range(32):
                data.append(int(linestr[j])) #添加數據
        fp.close()
    return np.array(data)

def loadImages(path):
    # 讀取數據
    hwLabels = []
    filelist = listdir(path)  # 獲取所有文件目錄
    m = len(filelist)
    data = np.zeros((m, 1024))

    # 先獲取標簽值
    for i in range(m):
        filename = filelist[i]
        classNumStr = int(filename.split('_')[0])
        if classNumStr == 9: hwLabels.append(-1)
        else: hwLabels.append(1)
        data[i, :] = img2vector("%s/%s" % (path, filename))

    return data, hwLabels

（三）測試數據集錯誤率

def testDights(data_X,data_Y,C,toler,maxIter,kTup):
    #獲取閾值b和α值
    b,alphas = smoP(data_X,data_Y,C,toler,maxIter,kTup)
    #獲取支持向量，即α>0點
    svIdx = np.where(alphas>0)[0]
    sVs = data_X[svIdx]
    labelSV = data_Y[svIdx]
    print("there are %d support vectors"%sVs.shape[0])
    m,n = data_X.shape
    errorCount = 0

    #獲取數據異常錯誤率
    for i in range(m):
        #先進行預測，預測公式見西瓜書6.24
        kernelEval = kernelTrans(sVs,data_X[i],kTup)    #使用支持向量點，和對應的預測數據點，進行預測
        pred = [e-mail](labelSV,alphas[svIdx]) + b #進行預測
        if np.sign(pred) != np.sign(data_Y[i]):
            errorCount += 1
    print("the trainning set error rate is: %f"%(errorCount/m))

    #測試集預測
    data_X, data_Y = loadImages("digits/testDigits")
    m,n = data_X.shape
    errorCount = 0

    #獲取數據異常錯誤率
    for i in range(m):
        #先進行預測，預測公式見西瓜書6.24
        kernelEval = kernelTrans(sVs,data_X[i],kTup)    #使用支持向量點，和對應的預測數據點，進行預測
        pred = [e-mail](labelSV,alphas[svIdx]) + b #進行預測
        if np.sign(pred) != np.sign(data_Y[i]):
            errorCount += 1
    print("the testing set error rate is: %f"%(errorCount/m))

（四）進行測試

data_X,data_Y = loadImages("digits/trainingDigits")
data_Y = np.array(data_Y)   #注意我們loadImages返回的data_Y是列表，我們要轉化為numpy數組類型
C = 200
toler = 0.0001
maxIter = 10000
kTup = ('rbf',10)

testDights(data_X,data_Y,C,toler,maxIter,kTup)