机器学习（3）-支持向量机（SVM）含python代码_python支持向量机给定确定的c和g

                                          支持向量机（SVM）

定义：

支持向量机：英文名support vector machine，一般简称SVM。通俗来讲，它是一种二类分类模型，其基本模型定义为特征空间上的间隔最大的超平面（线性分类器），其学习策略便是间隔最大化，最终可转化为一个凸二次规划问题的求解。

超平面方程式： ，其中X为多维向量，W也为多维向量，b为实数。当f(x)=0时，该点位于超平面上，根据f(x)与0的比较，确定其所属标签。

求解过程：

1，找出支持向量（已知的含有标签数据集中，每个标签下的边界向量）

2，根据支持向量，求出最优超平面（最优的标准就是支持向量距离超平面间隔最大）

本文以实际应用为基础，对一维、二维、多维的线性数据集和非线性数据集进行求解 线性数据集：能用一个超平面进行数据划分的数据集，在一维中该超平面为一个点，二维中为一条直线，三维中为一平面 非线性数据集：无法用一个超平面进行数据划分的数据集，与线性数据集对应。

一维数据

有一维数据T（ 1,5,6,10,8 ,20,30,50,60,75），其中A（1,5,6,9,8）的标签为1，B（ 20,30,50,60,75 ）的标签为-1，如下图： 

      

支持向量：A中边界点为：1,10。且分类器必然处于A与B之间，所以10为A的支持向量。同理可得20为B的支持变量 将10,20带入计算公式 ：                         

可得：w为1维向量，设为[X]

[X] .T*[10] + b =-1

[X] .T*[20] + b=1

求得：X=0.2(既w=[0.2])，b=-3 同样根据f(x)=0时为分类器上的点，可知点[15]为分类器上的点 。

二维数据

有二维数据集T，其中A集合的标签为-1，B集合标签为1，如下图：

                            

支持向量：A的备选边界点有a1，a2，a3，B中备选边界点有b1，b2（边界点求法多维时介绍），对分别做成直线a1a2，a2a3，b1b2，分别求相反标签中距离该直线最近点到该直线的距离(几何间隔，既垂线距离)，既b1到直线a2a3的距离dis1，b2到a1a2的距离dis2，a2到b1b2的距离dis3，比较可得dis3>dis2>dis1.既可以的得到A的支持向量为a2，B的支持向量为b1，b2。

将a2,b1,b2带入计算公式： ，w的转置为2维向量，设为[[X],[Y]]

[[X],[Y]].T*[[9],[4]] + b = 1

[[X],[Y]] .T*[[8],[1]] + b = 1

[[X],[Y]] .T*[[6],[5]] + b = -1

计算结果：X=0.6，Y=-0.2，b=-3.6， w=[0.6,-0.2]

同样根据f(x)=0时为分类器上的点，可知0.6X-0.2Y-3.6=0为分类直线

多维数据

有多维数据集T，其中A集合的标签为-1，B集合标签为1，由于无法表达多维数据，暂用二维数据图：

                           

多维数据求解步骤如下：

1，几何间隔：样本x到超平面的距离，公式为：

2，最大间隔分离器：既为我们所求的做好的分离器，其定义为离支持变量最远的超平面，当我们取支持变量值为1时，可以得到我们所求的问题转化为：  
也可以转换公式为： 

3，对偶算法：（1）引进拉格朗日乘子,定义拉格朗日函数

2）根据拉格朗日对偶性,原始问题的对偶问题是极大极小问题，先求对w,b的极小值.将L(w,b,a)分别对w,b求偏导数并令其等于0，带入公式。将求极大转换为求极小，最后转换公式为：

3）由KKT条件成立得到以下公式，其中j为使aj*>0的下标之一.所以问题就变为求对偶问题的解a*,再求得原始问题的解w*,b*,从而得分离超平面及分类决策函数可以看出w*和b*都只依赖训练数据中ai*>0的样本点(xi,yi),这些实例点xi被称为支持向量. 

                                                                                    

                                                        非线性数据

核函数：

支持向量机通过某非线性变换 φ( x) ，将输入空间映射到高维特征空间。特征空间的维数可能非常高。如果支持向量机的求解只用到内积运算，而在低维输入空间又存在某个函数 K(x, x′) ，它恰好等于在高维空间中这个内积，即K( x, x′) =<φ( x) ⋅φ( x′) > 。那么支持向量机就不用计算复杂的非线性变换，而由这个函数 K(x, x′) 直接得到非线性变换的内积，使大大简化了计算。这样的函数 K(x, x′) 称为核函数。

核技巧:基本思想是通过一个非线性变换将输入空间对应于一个特征空间,使得在输入空间中的超曲面模型对应于特征空间中的超平面模型(支持向量机).在学习和预测中只定义核函数K(x,z),而不显式地定义映射函数.对于给定的核K(x,z),特征空间和映射函数的取法并不唯一.注意到在线性支持向量机的对偶问题中,目标函数和决策函数都只涉及输入实例与实例之间的内积,xi`xj可以用核函数K(xi,xj)=Ф(xi)`Ф(xj)来代替.当映射函数是非线性函数时,学习到的含有核函数的支持向量机是非线性分类模型.在实际应用中,往往依赖领域知识直接选择核函数.

常见核函数

多项式核函数：

对应的支持向量机是一个p次多项式分类器,分类决策函数为：

高斯核函数：

对应的支持向量机是高斯径向基函数(RBF)分类器 分类决策函数为：

线性核函数：

此函数对应的为线性可分数据，这实际上就是原始空间中的内积，这个核存在的主要目的是使得“映射后空间中的问题”和“映射前空间中的问题”两者在形式上统一起来了(意思是说，咱们有的时候，写代码，或写公式的时候，只要写个模板或通用表达式，然后再代入不同的核，便可以了，于此，便在形式上统一了起来，不用再分别写一个线性的，和一个非线性的)

松弛变量

通过核函数把原始数据映射到高维空间之后，能够线性分隔的概率大大增加，但是对于某些情况还是很难处理，比如数据有噪音，对于这种偏离正常位置很远的数据点，我们称之为 outlier，如下图： 

                                                           

这些outlier意味着不能满足函数间隔大于等于1的约束条件,可以对每个样本点引进一个松弛变量,使函数间隔加上松弛变量大于等于1,约束条件变为：      

 同时对每个松弛变量,支付一个代价,目标函数变为;      ,其中C>0称为惩罚 

 参数,C值越大对误分类的惩罚也越大.新目标函数包含了两层含义:使间隔尽量大,同时使误分类点的个数尽量小.

根据前面的对偶性算法，最终所得公式为：

KKT条件成立可以得到,j是满足0<aj*<C的下标之一.问题就变为选择惩罚参数C>0,求得对偶问题(凸二次规划问题)的最优解a*,代入计算w*和b*,求得分离超平面和分类决策函数.因为b的解并不唯一,所以实际计算b*时可以取所有样本点上的平均值。

支持向量:在线性不可分的情况下,将对应与ai*>0的样本点(xi,yi)的实例点xi称为支持向量。

SMO算法

SMO(sequence minimum optimal 序列最小最优化):快速求解凸二次规划问题的算法.基本思路是:如果所有变量都满足此优化问题的KKT条件,那么解就得到了.否则,选择两个变量,固定其他变量,针对这两个变量构建一个二次规划问题.不断地将原问题分解为子问题并对子问题求解,就可以求解原问题.注意子问题两个变量中只有一个是自由变量,另一个由等式约束确定. 

两个变量二次规划的求解方法:       

假设选择的两个变量是a1,a2,其他变量是固定的,于是得到子问题:

ε是常数,目标函数式省略了不含a1,a2的常数项.考虑不等式约束和等式约束,要求的是目标函数在一条平行于对角线的线段上        的最优值

                                 

问题变为单变量的最优化问题.假设初始可行解为aold,最优解为anew,考虑沿着约束方向未经剪辑的最优解anew,unc(即未考虑不等式约束).对该问题求偏导数,并令导数为0,代入

 原式,，

 经剪辑后a2的解是  L与H是a2new所在的对角线段端点的界.并解得：
变量的选择方法:在每个子问题中选择两个变量优化,其中至少一个变量是违反KKT条件的.第一个变量的选取标准是违反KKT条件最严重的样本点,第二个变量的选取标准是希望能使该变量有足够大的变化,一般可以选取使对应的|E1-E2|最大的点.在每次选取完点后,更新阈值b和差值Ei.

算法步骤：

1，循环所以训练样本，判断是否满足KKT条件，如不满足，则记录该点ai： KKT条件：

                                                                                                                                   

（2）不等式约束使得(αi,αj)在盒子[0, C]x[0, C]内，等式约束使得(αi, αj)在平行于盒子[0, C]x[0, C]的对角线的直线上。因此要求的是目标函数在一条平行于对角线的线段上的最优值。这使得两个变量的最优化问题成为实质的单变量的最优化问题。由图可以得到，αj的上下界可以通过下面的方法得到： 

                                                                   

我们优化的时候，αj必须要满足上面这个约束。也就是说上面是αj的可行域。然后我们开始寻找αj，使得目标函数最大化。通过推导得到αj的更新公式如下 ：

                                       

对应参数值：

                              

这里Ek可以看做对第k个样本，SVM的输出与期待输出，也就是样本标签的误差。 而η实际上是度量两个样本i和j的相似性的。在计算η的时候，我们需要使用核函数，那么就可以用核函数来取代上面的内积。

（3）进行不等式剪裁：

                                     

（4）计算αi：

                                        

4，计算阈值b

优化αi和αj后，我们就可以更新阈值b，使得对两个样本i和j都满足KKT条件。如果优化后αi不在边界上（也就是满足0<αi<C，这时候根据KKT条件，可以得到yigi(xi)=1，这样我们才可以计算b），那下面的阈值b1是有效的，因为当输入xi时它迫使SVM输出yi。

          

同样      

如果0<αi<C和0<αj<C都满足，那么b1和b2都有效，而且他们是相等的。如果他们两个都处于边界上（也就是αi=0或者αi=C，同时αj=0或者αj=C），那么在b1和b2之间的阈值都满足KKT条件，一般我们取他们的平均值b=(b1+b2)/2。所以，总的来说对b的更新如下：

                   

5，凸优化问题终止条件：

SMO算法的基本思路是：如果说有变量的解都满足此最优化问题的KKT条件，那么这个最优化问题的解就得到了。因为KKT条件是该最优化问题的充分必要条件（证明请参考文献）。所以我们可以监视原问题的KKT条件，所以所有的样本都满足KKT条件，那么就表示迭代结束了。但是由于KKT条件本身是比较苛刻的，所以也需要设定一个容忍值，即所有样本在容忍值范围内满足KKT条件则认为训练可以结束；当然了，对于对偶问题的凸优化还有其他终止条件，可以参考文献。

代码：

from numpy import *
import matplotlib.pyplot as plt
 
MOSHI = 'mul'
GAOSI_C = 3
def loadDataSet(filename): #读取数据
    dataMat=[]
    labelMat=[]
    fr=open(filename)
    for line in fr.readlines():
        lineArr=line.strip().split('###')
        dataMat.append([float(lineArr[0]),float(lineArr[1])])
        labelMat.append(float(lineArr[2]))
    return dataMat,labelMat #返回数据特征和数据类别
 
def selectJrand(i,m): #在0-m中随机选择一个不是i的整数
    j=i
    while (j==i):
        j=int(random.uniform(0,m))
    return j
 
def clipAlpha(aj,H,L):  #保证a在L和H范围内（L <= a <= H）
    if aj>H:
        aj=H
    if L>aj:
        aj=L
    return aj
 
def kernelTrans(X, A, kTup): #核函数，输入参数,X:支持向量的特征树；A：某一行特征数据；kTup：('lin',k1)核函数的类型和参数
    m,n = shape(X)
    K = mat(zeros((m,1)))
    if kTup[0]=='lin': #线性函数
        K = X * A.T
    elif kTup[0]=='rbf': # 径向基函数(radial bias function)gaosi
        for j in range(m):
            deltaRow = X[j,:] - A
            K[j] = deltaRow*deltaRow.T
        K = exp(K/(-1*kTup[1]**2)) #返回生成的结果
    elif kTup[0]=='mul': # duo xiang shi
        for j in range(m):
            deltaRow = X[j,:]*A.T
            K[j] = (deltaRow + 1)**kTup[1]
    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):  # 存储各类参数
        self.X = dataMatIn  #数据特征
        self.labelMat = classLabels #数据类别
        self.C = C #软间隔参数C，参数越大，非线性拟合能力越强
        self.tol = toler #停止阀值
        self.m = shape(dataMatIn)[0] #数据行数
        self.alphas = mat(zeros((self.m,1)))
        self.b = 0 #初始设为0
        self.eCache = mat(zeros((self.m,2))) #缓存
        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): #计算Ek（参考《统计学习方法》p127公式7.105）
    fXk = float(multiply(oS.alphas,oS.labelMat).T*oS.K[:,k] + oS.b)
    Ek = fXk - float(oS.labelMat[k])
    return Ek
 
#随机选取aj，并返回其E值
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:
            if k == i:
                continue
            Ek = calcEk(oS, k)
            deltaE = abs(Ei - Ek)
            if (deltaE > maxDeltaE): #返回步长最大的aj
                maxK = k
                maxDeltaE = deltaE
                Ej = Ek
        return maxK, Ej
    else:
        j = selectJrand(i, oS.m)
        Ej = calcEk(oS, j)
    return j, Ej
 
 
def updateEk(oS, k): #更新os数据
    Ek = calcEk(oS, k)
    oS.eCache[k] = [1,Ek]
 
#首先检验ai是否满足KKT条件，如果不满足，随机选择aj进行优化，更新ai,aj,b值
def innerL(i, oS): #输入参数i和所有参数数据
    Ei = calcEk(oS, i) #计算E值
    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)): #检验这行数据是否符合KKT条件 参考《统计学习方法》p128公式7.111-113
        j,Ej = selectJ(i, oS, Ei) #随机选取aj，并返回其E值
        alphaIold = oS.alphas[i].copy()
        alphaJold = oS.alphas[j].copy()
        if (oS.labelMat[i] != oS.labelMat[j]): #以下代码的公式参考《统计学习方法》p126
            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] #参考《统计学习方法》p127公式7.107 (x-y)**2
        if eta >= 0:
            print("eta>=0")
            return 0
        oS.alphas[j] -= oS.labelMat[j]*(Ei - Ej)/eta #参考《统计学习方法》p127公式7.106
        oS.alphas[j] = clipAlpha(oS.alphas[j],H,L) #参考《统计学习方法》p127公式7.108
        updateEk(oS, j)
        if (abs(oS.alphas[j] - alphaJold) < oS.tol): #alpha变化大小阀值（自己设定）
            print("j not moving enough")
            return 0
        oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJold - oS.alphas[j])#参考《统计学习方法》p127公式7.109
        updateEk(oS, i) #更新数据
        #以下求解b的过程，参考《统计学习方法》p129公式7.114-7.116
        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]<oS.C):
            oS.b = b1
        elif (0 < oS.alphas[j]<oS.C):
            oS.b = b2
        else:
            oS.b = (b1 + b2)/2.0
        return 1
    else:
        return 0
 
 
#SMO函数，用于快速求解出alpha
def smoP(dataMatIn, classLabels, C, toler, maxIter,kTup=('lin', 0)): #输入参数：数据特征，数据类别，参数C，阀值toler，最大迭代次数，核函数（默认线性核）
    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:
            for i in range(oS.m): #遍历所有数据
                alphaPairsChanged += innerL(i,oS)
                print("fullSet, iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged)) #显示第多少次迭代，那行特征数据使alpha发生了改变，这次改变了多少次alpha
        else:
            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
        elif (alphaPairsChanged == 0):
            entireSet = True
        print("iteration number: %d" % iter)
    return oS.b,oS.alphas
 
def showdata(dataArr,labelArr, ax):
    datMat = array(dataArr)
    labelMat = array(labelArr)
    labelMat_i = argwhere(labelMat == 1)
    labelMat_j = argwhere(labelMat == -1)
    datMat_i = datMat[labelMat_i].T
    datMat_j = datMat[labelMat_j].T
 
    plt.sca(ax)
    plt.scatter([datMat_i[1]], [datMat_i[0]], c='r', marker='x')
    plt.scatter([datMat_j[1]], [datMat_j[0]], c='g', marker='o')
 
def showdata_i(dataArr,nums, ax):
    datMat = array(dataArr)
    datMat_i = datMat[nums].T
    plt.sca(ax)
    plt.scatter([datMat_i[1]], [datMat_i[0]], c='b', marker='v')
 
def showdata_sv(dataArr,nums, ax):
    datMat = array(dataArr)
    datMat_i = datMat[nums].T
    plt.sca(ax)
    plt.scatter([datMat_i[1]], [datMat_i[0]], c='b', marker='P')
 
 
def testRbf(data_train,data_test):
    plt.figure(1, figsize=(16, 8))
    dataArr,labelArr = loadDataSet(data_train) #读取训练数据
    ax1 = plt.subplot(1, 2, 1)
    showdata(dataArr,labelArr, ax1)
    b,alphas = smoP(dataArr, labelArr, 200, 0.00001, 10000, (MOSHI, GAOSI_C)) #通过SMO算法得到b和alpha
    datMat=mat(dataArr)
    labelMat = mat(labelArr).transpose()
    svInd=nonzero(alphas)[0]  #选取不为0数据的行数（也就是支持向量）
    showdata_sv(dataArr, svInd, ax1)
    sVs=datMat[svInd] #支持向量的特征数据
    labelSV = labelMat[svInd] #支持向量的类别（1或-1）00001
    print("there are %d Support Vectors" % shape(sVs)[0]) #打印出共有多少的支持向量
    m,n = shape(datMat) #训练数据的行列数
    errorCount = 0
    wT = multiply(labelSV, alphas[svInd]).T
    nums = []
    for i in range(m):
        kernelEval = kernelTrans(sVs,datMat[i,:],(MOSHI, GAOSI_C)) #将支持向量转化为核函数
        predict=wT * kernelEval + b  #这一行的预测结果（代码来源于《统计学习方法》p133里面最后用于预测的公式）注意最后确定的分离平面只有那些支持向量决定。
        if sign(predict)!=sign(labelArr[i]): #sign函数 -1 if x < 0, 0 if x==0, 1 if x > 0
            errorCount += 1
            nums.append(i)
    showdata_i(dataArr, nums, ax1)
    print("the training error rate is: %f" % (float(errorCount)/m)) #打印出错误率
    dataArr_test,labelArr_test = loadDataSet(data_test) #读取测试数据
    ax2 = plt.subplot(1, 2, 2)
    showdata(dataArr_test, labelArr_test, ax2)
    errorCount_test = 0
    datMat_test=mat(dataArr_test)
    labelMat = mat(labelArr_test).transpose()
    m,n = shape(datMat_test)
    nums_test = []
    for i in range(m): #在测试数据上检验错误率
        kernelEval = kernelTrans(sVs,datMat_test[i,:],(MOSHI, GAOSI_C))
        predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + b
        if sign(predict)!=sign(labelArr_test[i]):
            errorCount_test += 1
            nums_test.append(i)
    print("the test error rate is: %f" % (float(errorCount_test)/m))
    showdata_i(dataArr_test, nums_test, ax2)
    plt.show()
 
#主程序
def main():
    filename_traindata='train_data2.TXT'
    filename_testdata='test_data1.TXT'
    testRbf(filename_traindata,filename_testdata)
 
if __name__=='__main__':
    main()