Svm算法原理及實現

      Svm（support Vector Mac）又稱爲支持向量機，是一種二分類的模型。固然若是進行修改以後也是能夠用於多類別問題的分類。支持向量機能夠分爲線性核非線性兩大類。其主要思想爲找到空間中的一個更夠將全部數據樣本劃開的超平面，而且使得本本集中全部數據到這個超平面的距離最短。

1、基於最大間隔分隔數據

1.1支持向量與超平面

    在瞭解svm算法以前，咱們首先須要瞭解一下線性分類器這個概念。好比給定一系列的數據樣本，每一個樣本都有對應的一個標籤。爲了使得描述更加直觀，咱們採用二維平面進行解釋，高維空間原理也是同樣。舉個簡單子：以下圖所示是一個二維平面，平面上有兩類不一樣的數據，分別用圓圈和方塊表示。咱們能夠很簡單地找到一條直線使得兩類數據正好可以徹底分開。可是能將據點徹底劃開直線不止一條，那麼在如此衆多的直線中咱們應該選擇哪一條呢？從直觀感受上看圖中的幾條直線，是否是要更好一些呢？是的咱們就是但願尋找到這樣的直線，使得距離這條直線最近的點到這條直線的距離最短。這讀起來有些拗口，咱們從圖三中直觀來解釋這一句話就是要求的兩條外面的線之間的間隔最大。這是能夠理解的，由於假如數據樣本是隨機出現的，那麼這樣分割以後數據點落入到其類別一側的機率越高那麼最終預測的準確率也會越高。在高維空間中這樣的直線稱之爲超平面，由於當維數大於三的時候咱們已經沒法想象出這個平面的具體樣子。那些距離這個超平面最近的點就是所謂支持向量，實際上若是肯定了支持向量也就肯定了這個超平面，找到這些支持向量以後其餘樣本就不會起做用了。

                                  圖 1                                                 圖2

 

1.2尋找最大間隔

1.2.1點到超平面的距離公式

      既然這樣的直線是存在的，那麼咱們怎樣尋找出這樣的直線呢？與二維空間相似，超平面的方程也能夠寫成一下形式：

                                                                                             （1.1）

有了超平面的表達式以後以後，咱們就能夠計算樣本點到平面的距離了。假設爲樣本的中的一個點，其中表示爲第個特徵變量。那麼該點到超平面的距離就能夠用以下公式進行計算：

                                                                         （1.2）

 

其中||W||爲超平面的範數，常數b相似於直線方程中的截距。

上面的公式能夠利用解析幾何或高中平面幾何知識進行推導，這裏不作進一步解釋。

1.2.2最大間隔的優化模型

    如今咱們已經知道了如何去求數據點到超平面的距離，在超平面肯定的狀況下，咱們就可以找出全部支持向量，而後計算出間隔margin。每個超平面都對應着一個margin，咱們的目標就是找出全部margin中最大的那個值對應的超平面。所以用數學語言描述就是肯定w、b使得margin最大。這是一個優化問題其目標函數能夠寫成：

                                         （1.3）

其中表示數據點的標籤，且其爲-1或1。距離用計算，這是就能體會出-1和1的好處了。若是數據點在平面的正方向(即+1類)那麼是一個正數，而當數據點在平面的負方向時(即-1類)，依然是一個正數，這樣就可以保證始終大於零了。注意到當w和b等比例放大時，d的結果是不會改變的。所以咱們能夠令全部支持向量的u爲1，而其餘點的u大1這是能夠辦經過調節w和b求到的。所以上面的問題能夠簡化爲：                     （1.4）

爲了後面計算的方便，咱們將目標函數等價替換爲：

                                                               （1.5）

這是一個有約束條件的優化問題，一般咱們能夠用拉格朗日乘子法來求解。拉格朗日乘子法的介紹能夠參考這篇博客。應用拉格朗日乘子法以下：

令                                （1.6）

求L關於求偏導數得：                          （1.7）

將(1.7)代入到(1.6)中化簡得：

                                      （1.8）

原問題的對偶問題爲：

                                              （1.9）

該對偶問題的KKT條件爲

                                （1.10）

    到此，彷佛問題就可以完美地解決了。可是這裏有個假設：數據必須是百分之百可分的。可是實際中的數據幾乎都不那麼「乾淨」，或多或少都會存在一些噪點。爲此下面咱們將引入了鬆弛變量來解決這種問題。

1.2.3鬆弛變量

    由上一節的分析咱們知道實際中不少樣本數據都不可以用一個超平面把數據徹底分開。若是數據集中存在噪點的話，那麼在求超平的時候就會出現很大問題。從圖三中課看出其中一個藍點誤差太大，若是把它做爲支持向量的話所求出來的margin就會比不算入它時要小得多。更糟糕的狀況是若是這個藍點落在了紅點之間那麼就找不出超平面了。

                                   

                                                         圖 3

所以引入一個鬆弛變量ξ來容許一些數據能夠處於分隔面錯誤的一側。這時新的約束條件變爲:

  （1.11）

式中ξi的含義爲容許第i個數據點容許偏離的間隔。若是讓ξ任意大的話，那麼任意的超平面都是符合條件的了。因此在原有目標的基礎之上，咱們也儘量的讓ξ的總量也儘量地小。因此新的目標函數變爲：

（1.12）

（1.13）

其中的C是用於控制「最大化間隔」和「保證大部分的點的函數間隔都小於1」這兩個目標的權重。將上述模型完整的寫下來就是：

（1.14）

新的拉格朗日函數變爲：

（1.15）

接下來將拉格朗日函數轉化爲其對偶函數，首先對分別求ξ的偏導，並令其爲0,結果以下：

（1.16）

代入原式化簡以後獲得和原來同樣的目標函數：

（1.17）

可是因爲咱們獲得而，所以有因此對偶問題寫成：

（1.18）

通過添加鬆弛變量的方法，咱們如今可以解決數據更加混亂的問題。經過修改參數C，咱們能夠獲得不一樣的結果而C的大小到底取多少比較合適，須要根據實際問題進行調節。

1.2.4核函數

    以上討論的都是在線性可分狀況進行討論的，可是實際問題中給出的數據並非都是線性可分的，好比有些數據多是如圖4樣子。      

                                                          圖4

那麼這種非線性可分的數據是否就不能用svm算法來求解呢？答案是否認的。事實上，對於低維平面內不可分的數據，放在一個高維空間中去就有可能變得可分。以二維平面的數據爲例，咱們能夠經過找到一個映射將二維平面的點放到三維平面之中。理論上任意的數據樣本都可以找到一個合適的映射使得這些在低維空間不能劃分的樣本到高維空間中以後可以線性可分。咱們再來看一下以前的目標函數：

（1.19）

定義一個映射使得將全部映射到更高維空間以後等價於求解上述問題的對偶問題：

（1.20）

這樣對於線性不可分的問題就解決了，如今只須要找出一個合適的映射便可。當特徵變量很是多的時候在，高維空間中計算內積的運算量是很是龐大的。考慮到咱們的目的並非爲找到這樣一個映射而是爲了計算其在高維空間的內積，所以若是咱們可以找到計算高維空間下內積的公式，那麼就可以避免這樣龐大的計算量，咱們的問題也就解決了。實際上這就是咱們要找的核函數，即兩個向量在隱式映射後的空間中的內積。下面的一個簡單例子能夠幫助咱們更好地理解核函數。

經過以上例子，咱們能夠很明顯地看到核函數是怎樣運做的。上述問題的對偶問題能夠寫成以下形式：

（1.21）

那麼怎樣的函數才能夠做爲核函數呢？下面的一個定理能夠幫助咱們判斷。

Mercer定理：任何半正定的函數均可以做爲核函數。其中所謂半正定函數是指擁有訓練集數據集合，咱們定義一個矩陣的元素，這個矩陣是的矩陣，若是這個矩陣是半正定的，那麼就稱爲半正定函數。

值得注意的是，上述定理中所給出的條件是充分條件而非充要條件。由於有些非正定函數也能夠做爲核函數。

下面是一些經常使用的核函數：

                                             表1 經常使用核函數表

核函數名稱	核函數表達式	核函數名稱	核函數表達式
線性核		指數核
多項式核		拉普拉斯核
高斯核		Sigmoid核

    如今咱們已經瞭解了一些支持向量機的理論基礎，咱們經過對偶問題的的轉化將最開始求的問題轉化爲求的對偶問題。只要找到全部的(即找出全部支持向量)，咱們就可以肯定。而後就能夠經過計算數據點到這個超平面的距離從而判斷出該數據點的類別。

2、Smo算法原理

2.1 約束條件

     根據以上問題的分析，咱們已經將原始問題轉化爲了求的值，即求下面優化模型的解：

（2.1）

求解的值的方法有不少，Smo算法就是其中一種比較經常使用的方法。該算法是由John Platt在1996年發佈，他的思路是將大的優化問題轉化爲小的優化問題。而這些小的優化問題每每更容易求解，而且對它們進行順序求解的結果和將它們做爲總體求解的結果徹底一致可是Smo算法的時間要小得多。

     Smo算法的原理爲：每次任意抽取兩個乘子和，而後固定和之外的其它乘子，使得目標函數只是關於和的函數。而後增大其中一個乘子同時減小另一個。這樣，不斷的從一堆乘子中任意抽取兩個求解，不斷的迭代求解子問題，最終達到求解原問題的目的。

     而原對偶問題的子問題的目標函數能夠表達成：

（2.2）

其中：

（2.3）

這裏之因此算兩個是由於的限制，若是隻改變其中的一個量，那麼這個約束條件可能就不成立了。要解決這個問題，咱們必須得選取這樣的兩個乘子。那麼怎樣肯定這樣的和呢？這是咱們首先要考慮的問題，在《機器學習實戰》這本書中，做者首先給出了一種簡化版的方法，遍歷每個而後在剩餘的中隨機選取一個進行優化。雖然樣也可以解決問題，可是運算量太大，所以考慮找一種更好的方法來尋找對。

爲了表述方便，定義一個特徵到輸出結果的輸出函數：

（2.4）

該對偶問題中KKT條件爲：

（2.5）

根據上述問題的KKT條件能夠得出目標函數中的的含義以下：

一、 ，代表是正常分類，在邊界外；

二、，代表是支持向量，在邊界上；

三、，代表在兩邊界之間。

最優解須要知足KKT條件，所以須要知足以上的三個條件都知足。而不知足這三個條件的狀況也有三種：

一、<=1可是<C則是不知足的，而本來=C;

二、>=1可是>0則是不知足的，而本來=0;

三、=1可是=0或者=C則代表不知足的，而本來應該是0<<C.

也就是說若是存在不知足這些KKT條件的，咱們就要更新它，這就是約束條件之一。其次，還受到約束條件的限制，所以假設選擇的兩個因子爲和，他們在更新前分別爲在更新後爲，爲了保證上述約束條件成立必需要保證下列等式的成立：

（2.6）

其中爲常數。

2.2參數優化

    由於兩個因子很差同時求解，因此可先求第二個乘子的解，而後再用的解表示的解。爲了求解，得先肯定的取值範圍。假設它的上下邊界分別爲H和L，那麼有：（2.6）

接下來，綜合和這兩個約束條件，求取的取值範圍。

當時，根據可得，因此有。

當時，一樣根據可得：，因此有。

回顧第二個約束條件 ：（2.7）

令其兩邊同時乘y1，可得：

                          . （2.8）

其中：.

所以能夠用表示，即：               （2.9）

令                    （2.10）

通過轉化以後可得：

                   ,.        (2.11)

那麼如何來選擇乘子呢？對於第一個乘子，咱們能夠按照3種不知足KTT的條件來尋找。對於第二個乘子，咱們能夠尋找知足條件的乘子。

而b在知足如下條件時須要更新：

                       (2.12)

且更新後的和以下：

                      (2.13）

每次更新完兩個乘子以後，都須要從新計算b以及對應的E。最後更新完全部的,y和b,這樣模型也就出來了，從而能夠計算出開始所說的分類函數：

（2.14）

3、Svm的python實現與應用

     第二節已經對Smo算法進行了充分的說明並進行了推導，如今一切都準備好了。接下來須要作的就是實現這些算法了，這裏我參考了《機器學習實戰》這本書中的代碼，並利用該程序對一個問題進行了求解。因爲代碼數量過大，所以這裏就再也不列出，而是放在附錄中。有興趣的朋友能夠直接下載，也能夠去官網下載源代碼。筆者在讀這些代碼的過程當中，也遇到了許多困難，大部分都根據本身的狀況進行了註釋。

3.1問題描述

    這裏我選取的一個數據集爲聲吶數據集，該問題爲須要根據聲吶從不一樣角度返回的聲音強度來預測被測物體是岩石仍是礦井。數據集中共有208個數據，60個輸入變量和1個輸出變量。每一個數據的前60個元素爲不一樣方向返回的聲音強度，最後一個元素爲本次用於聲吶測試的物體類型。其中標籤M表示礦井，標籤R爲岩石。

3.2數據預處理

    所給數據中沒有缺失值和異常值，所以不須要對數據集進行清洗。注意到所給數據集的標籤爲字母類型，而svm算法的標準標籤爲「-1」和「1」的形式，因此須要對標籤進行轉化，用「-1」、「1」分別代替M和R。因爲該數據集中所給標籤相同的數據都放在了一塊兒，所以先對數據順序進行了打亂。代碼以下：

def loadDataSet(fileName):    dataMat=[];labelMat=[];data=[]    fr=open(fileName)    for line in fr.readlines():        line=line.replace(',','\t')#去除逗號        line=line.replace('M','-1')#對標籤進行替換        line=line.replace('R','1')        lineArr=line.strip('\n').split('\t')#分割數據        data.append([float(lineArr[i]) for i in range(len(lineArr))])    random.shuffle(data)  #隨機打亂列表    for i in range(len(data)):        dataMat.append(data[i][0:len(data)-1])        labelMat.append(data[i][-1])    return dataMat,labelMat

3.3模型訓練及測試

    首先測試了一下將全部數據即做爲訓練集又做爲測試集，而後用Smo模型進行訓練找到全部的支持向量。最後根據訓練好的模型進行求解，最終測試的準確率平均爲54%。若是把數據集分紅測試集和訓練集，發現測試的準確率也在這附近。而根據網上數據統計該數據集測試的準確率最高爲84%，通常狀況下爲百分之六十幾。數據集自己是形成測試準確率低的一個緣由，可是另一個更加劇要的緣由大概是參數的選擇問題。如何選取合適的參數是一個值得思考的問題，在接下來的學習過程當中我也會注意一下參數選取這個問題。到這裏，關於svm的算法就告一段落了。

#svm算法的實現
from numpy import*
import random
from time import*
def loadDataSet(fileName):#輸出dataArr(m*n),labelArr(1*m)其中m爲數據集的個數
    dataMat=[];labelMat=[]
    fr=open(fileName)
    for line in fr.readlines():
        lineArr=line.strip().split('\t')#去除製表符，將數據分開
        dataMat.append([float(lineArr[0]),float(lineArr[1])])#數組矩陣
        labelMat.append(float(lineArr[2]))#標籤
    return dataMat,labelMat

def selectJrand(i,m):#隨機找一個和i不一樣的j
    j=i
    while(j==i):
        j=int(random.uniform(0,m))
    return j

def clipAlpha(aj,H,L):#調整大於H或小於L的alpha的值
    if aj>H:
        aj=H
    if aj<L:
        aj=L
    return aj

def smoSimple(dataMatIn,classLabels,C,toler,maxIter):
    dataMatrix=mat(dataMatIn);labelMat=mat(classLabels).transpose()#轉置
    b=0;m,n=shape(dataMatrix)#m爲輸入數據的個數，n爲輸入向量的維數
    alpha=mat(zeros((m,1)))#初始化參數，肯定m個alpha
    iter=0#用於計算迭代次數
    while (iter<maxIter):#當迭代次數小於最大迭代次數時（外循環）
        alphaPairsChanged=0#初始化alpha的改變量爲0
        for i in range(m):#內循環
            fXi=float(multiply(alpha,labelMat).T*\
                      (dataMatrix*dataMatrix[i,:].T))+b#計算f(xi)
            Ei=fXi-float(labelMat[i])#計算f(xi)與標籤之間的偏差
            if ((labelMat[i]*Ei<-toler)and(alpha[i]<C))or\
                    ((labelMat[i]*Ei>toler)and(alpha[i]>0)):#若是能夠進行優化
                j=selectJrand(i,m)#隨機選擇一個j與i配對
                fXj=float(multiply(alpha,labelMat).T*\
                          (dataMatrix*dataMatrix[j,:].T))+b#計算f(xj)
                Ej=fXj-float(labelMat[j])#計算j的偏差
                alphaIold=alpha[i].copy()#保存原來的alpha(i)
                alphaJold=alpha[j].copy()
                if(labelMat[i]!=labelMat[j]):#保證alpha在0到c之間
                    L=max(0,alpha[j]-alpha[i])
                    H=min(C,C+alpha[j]-alpha[i])
                else:
                    L=max(0,alpha[j]+alpha[i]-C)
                    H=min(C,alpha[j]+alpha[i])
                if L==H:print('L=H');continue
                eta=2*dataMatrix[i,:]*dataMatrix[j,:].T-\
                    dataMatrix[i,:]*dataMatrix[i,:].T-\
                    dataMatrix[j,:]*dataMatrix[j,:].T
                if eta>=0:print('eta=0');continue
                alpha[j]-=labelMat[j]*(Ei-Ej)/eta
                alpha[j]=clipAlpha(alpha[j],H,L)#調整大於H或小於L的alpha
                if (abs(alpha[j]-alphaJold)<0.0001):
                    print('j not move enough');continue
                alpha[i]+=labelMat[j]*labelMat[i]*(alphaJold-alpha[j])
                b1=b-Ei-labelMat[i]*(alpha[i]-alphaIold)*\
                    dataMatrix[i,:]*dataMatrix[i,:].T-\
                    labelMat[j]*(alpha[j]-alphaJold)*\
                    dataMatrix[i,:]*dataMatrix[j,:].T#設置b
                b2=b-Ej-labelMat[i]*(alpha[i]-alphaIold)*\
                    dataMatrix[i,:]*dataMatrix[i,:].T-\
                    labelMat[j]*(alpha[j]-alphaJold)*\
                    dataMatrix[j,:]*dataMatrix[j,:].T
                if (0<alpha[i])and(C>alpha[j]):b=b1
                elif(0<alpha[j])and(C>alpha[j]):b=b2
                else:b=(b1+b2)/2
                alphaPairsChanged+=1
                print('iter:%d i:%d,pairs changed%d'%(iter,i,alphaPairsChanged))
        if (alphaPairsChanged==0):iter+=1
        else:iter=0
        print('iteraction number:%d'%iter)
    return b,alpha
#定義徑向基函數
def kernelTrans(X, A, kTup):#定義核轉換函數（徑向基函數）
    m,n = shape(X)
    K = mat(zeros((m,1)))
    if kTup[0]=='lin': K = X * A.T   #線性核K爲m*1的矩陣
    elif kTup[0]=='rbf':
        for j in range(m):
            deltaRow = X[j,:] - A
            K[j] = deltaRow*deltaRow.T
        K = exp(K/(-1*kTup[1]**2)) #divide in NumPy is element-wise not matrix like Matlab
    else: raise NameError('Houston We Have a Problem -- \
    That Kernel is not recognized')
    return K

class optStruct:
    def __init__(self,dataMatIn, classLabels, C, toler, kTup):  # Initialize the structure with the parameters
        self.X = dataMatIn
        self.labelMat = classLabels
        self.C = C
        self.tol = toler
        self.m = shape(dataMatIn)[0]
        self.alphas = mat(zeros((self.m,1)))
        self.b = 0
        self.eCache = mat(zeros((self.m,2))) #first column is valid flag
        self.K = mat(zeros((self.m,self.m)))
        for i in range(self.m):
            self.K[:,i] = kernelTrans(self.X, self.X[i,:], kTup)

def calcEk(oS, k):
    fXk = float(multiply(oS.alphas,oS.labelMat).T*oS.K[:,k] + oS.b)
    Ek = fXk - float(oS.labelMat[k])
    return Ek

def selectJ(i, oS, Ei):
    maxK = -1; maxDeltaE = 0; Ej = 0
    oS.eCache[i] = [1,Ei]
    validEcacheList = nonzero(oS.eCache[:,0].A)[0]
    if (len(validEcacheList)) > 1:
        for k in validEcacheList:   #loop through valid Ecache values and find the one that maximizes delta E
            if k == i: continue #don't calc for i, waste of time
            Ek = calcEk(oS, k)
            deltaE = abs(Ei - Ek)
            if (deltaE > maxDeltaE):
                maxK = k; maxDeltaE = deltaE; Ej = Ek
        return maxK, Ej
    else:   #in this case (first time around) we don't have any valid eCache values
        j = selectJrand(i, oS.m)
        Ej = calcEk(oS, j)
    return j, Ej

def updateEk(oS, k):#after any alpha has changed update the new value in the cache
    Ek = calcEk(oS, k)
    oS.eCache[k] = [1,Ek]

def innerL(i, oS):
    Ei = calcEk(oS, i)
    if ((oS.labelMat[i]*Ei < -oS.tol) and (oS.alphas[i] < oS.C)) or ((oS.labelMat[i]*Ei > oS.tol) and (oS.alphas[i] > 0)):
        j,Ej = selectJ(i, oS, Ei) #this has been changed from selectJrand
        alphaIold = oS.alphas[i].copy(); alphaJold = oS.alphas[j].copy()
        if (oS.labelMat[i] != oS.labelMat[j]):
            L = max(0, oS.alphas[j] - oS.alphas[i])
            H = min(oS.C, oS.C + oS.alphas[j] - oS.alphas[i])
        else:
            L = max(0, oS.alphas[j] + oS.alphas[i] - oS.C)
            H = min(oS.C, oS.alphas[j] + oS.alphas[i])
        if L==H: print("L==H"); return 0
        eta = 2.0 * oS.K[i,j] - oS.K[i,i] - oS.K[j,j] #changed for kernel
        if eta >= 0: print("eta>=0"); return 0
        oS.alphas[j] -= oS.labelMat[j]*(Ei - Ej)/eta
        oS.alphas[j] = clipAlpha(oS.alphas[j],H,L)
        updateEk(oS, j) #added this for the Ecache
        if (abs(oS.alphas[j] - alphaJold) < 0.00001): print("j not moving enough"); return 0
        oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJold - oS.alphas[j])#update i by the same amount as j
        updateEk(oS, i) #added this for the Ecache                    #the update is in the oppostie direction
        b1 = oS.b - Ei- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.K[i,i] - oS.labelMat[j]*(oS.alphas[j]-alphaJold)*oS.K[i,j]
        b2 = oS.b - Ej- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.K[i,j]- oS.labelMat[j]*(oS.alphas[j]-alphaJold)*oS.K[j,j]
        if (0 < oS.alphas[i]) and (oS.C > oS.alphas[i]): oS.b = b1
        elif (0 < oS.alphas[j]) and (oS.C > oS.alphas[j]): oS.b = b2
        else: oS.b = (b1 + b2)/2.0
        return 1
    else: return 0
#smoP函數用於計算超平的alpha,b
def smoP(dataMatIn, classLabels, C, toler, maxIter,kTup=('lin', 0)):    #完整的Platter SMO
    oS = optStruct(mat(dataMatIn),mat(classLabels).transpose(),C,toler, kTup)
    iter = 0#計算循環的次數
    entireSet = True; alphaPairsChanged = 0
    while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):
        alphaPairsChanged = 0
        if entireSet:   #go over all
            for i in range(oS.m):
                alphaPairsChanged += innerL(i,oS)
                print("fullSet, iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged))
            iter += 1
        else:#go over non-bound (railed) alphas
            nonBoundIs = nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0]
            for i in nonBoundIs:
                alphaPairsChanged += innerL(i,oS)
                print("non-bound, iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged))
            iter += 1
        if entireSet: entireSet = False #toggle entire set loop
        elif (alphaPairsChanged == 0): entireSet = True
        print("iteration number: %d" % iter)
    return oS.b,oS.alphas
#calcWs用於計算權重值w
def calcWs(alphas,dataArr,classLabels):#計算權重W
    X = mat(dataArr); labelMat = mat(classLabels).transpose()
    m,n = shape(X)
    w = zeros((n,1))
    for i in range(m):
        w += multiply(alphas[i]*labelMat[i],X[i,:].T)
    return w

#值得注意的是測試準確與k1和C的取值有關。
def testRbf(k1=1.3):#給定輸入參數K1
    #測試訓練集上的準確率
    dataArr,labelArr = loadDataSet('testSetRBF.txt')#導入數據做爲訓練集
    b,alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000, ('rbf', k1)) #C=200 important
    datMat=mat(dataArr); labelMat = mat(labelArr).transpose()
    svInd=nonzero(alphas.A>0)[0]#找出alphas中大於0的元素的位置
    #此處須要說明一下alphas.A的含義
    sVs=datMat[svInd] #獲取支持向量的矩陣，由於只要alpha中不等於0的元素都是支持向量
    labelSV = labelMat[svInd]#支持向量的標籤
    print("there are %d Support Vectors" % shape(sVs)[0])#輸出有多少個支持向量
    m,n = shape(datMat)#數據組的矩陣形狀表示爲有m個數據，數據維數爲n
    errorCount = 0#計算錯誤的個數
    for i in range(m):#開始分類，是函數的核心
        kernelEval = kernelTrans(sVs,datMat[i,:],('rbf', k1))#計算原數據集中各元素的核值
        predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + b#計算預測結果y的值
        if sign(predict)!=sign(labelArr[i]): errorCount += 1#利用符號判斷類別
        ### sign（a）爲符號函數：若a>0則輸出1，若a<0則輸出-1.###
    print("the training error rate is: %f" % (float(errorCount)/m))
    #二、測試測試集上的準確率
    dataArr,labelArr = loadDataSet('testSetRBF2.txt')
    errorCount = 0
    datMat=mat(dataArr)#labelMat = mat(labelArr).transpose()此處能夠不用
    m,n = shape(datMat)
    for i in range(m):
        kernelEval = kernelTrans(sVs,datMat[i,:],('rbf', k1))
        predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + b
        if sign(predict)!=sign(labelArr[i]): errorCount += 1
    print("the test error rate is: %f" % (float(errorCount)/m))
def main():
    t1=time()
    dataArr,labelArr=loadDataSet('testSet.txt')
    b,alphas=smoP(dataArr,labelArr,0.6,0.01,40)
    ws=calcWs(alphas,dataArr,labelArr)
    testRbf()
    t2=time()
    print("程序所用時間爲%ss"%(t2-t1))

if __name__=='__main__':
    main()

 

                                                   後記

    這是第一次寫博客，其中不免會出錯，所以但願你們可以批評指正。首先很是感謝網上的一些朋友，在理解svm這算法他們給了我不少啓發，在公式推導中給了我不少參考的地方。本文主要參考的資料是《機器學習實戰》和《驚呼！理解svm的三種境界》這篇博客。對於svm雖然學的時間不長，可是對它一直有種特別的感受。第一次據說svm是在作一個驗證碼識別問題的時候，但那時候我使用的是KNN算法，儘管效果還不錯，可是我一直但願可以用svm算法來完成這個題目。原本此次是打算把這個問題一塊兒解決的，可是因爲時間關係，沒有來得及作。只能等下次有空閒的時候再來作這個問題了。